1"""F2-conf-01: school-level programme effect with cluster-robust uncertainty.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json4import numpy as np5import pandas as pd6import statsmodels.formula.api as smf7from scipy import stats89d = pd.read_csv("data/schools.csv")10n_input = len(d)11issues = []12# Data checks13if d.isna().any().any():14 issues.append({"issue": "missing values", "handling": "rows with missing values dropped"})15 d = d.dropna()16if d["student_id"].duplicated().any():17 issues.append({"issue": "duplicate student_id", "handling": "duplicates dropped (first kept)"})18 d = d.drop_duplicates("student_id")19tconst = d.groupby("school_id")["treated"].nunique()20if (tconst > 1).any():21 issues.append({"issue": "treated not constant within school", "handling": "reported; not altered"})22sizes = d.groupby("school_id").size()2324G = d["school_id"].nunique()25G1 = int(d.groupby("school_id")["treated"].first().sum())2627ols = smf.ols("score ~ treated", data=d).fit()28tau = float(ols.params["treated"])29se_iid = float(ols.bse["treated"])3031# Cluster-robust (CR1, Stata-style small-sample correction) by school,32# inference with t distribution on G-1 degrees of freedom.33cl = smf.ols("score ~ treated", data=d).fit(cov_type="cluster",34 cov_kwds={"groups": d["school_id"], "use_correction": True, "df_correction": True})35se_cl = float(cl.bse["treated"])36df = G - 137tstat = tau / se_cl38p = float(2 * stats.t.sf(abs(tstat), df))39q = stats.t.ppf(0.975, df)40lo, hi = tau - q * se_cl, tau + q * se_cl4142# Cross-check: two-sample t-test on school means (equivalent analysis at school level)43sm = d.groupby("school_id").agg(y=("score", "mean"), t=("treated", "first"))44tt = stats.ttest_ind(sm.y[sm.t == 1], sm.y[sm.t == 0], equal_var=True)45se_means = float((sm.y[sm.t == 1].mean() - sm.y[sm.t == 0].mean()) / tt.statistic)4647# ICC of residual: one-way ANOVA estimator on OLS residuals by school48d = d.assign(r=ols.resid)49m = sizes.iloc[0] if sizes.nunique() == 1 else (len(d) - (sizes**2).sum() / len(d)) / (G - 1)50gm = d.groupby("school_id")["r"].mean()51msb = (sizes * (gm - d["r"].mean())**2).sum() / (G - 1)52msw = ((d["r"] - d["school_id"].map(gm))**2).sum() / (len(d) - G)53s2b = max((msb - msw) / m, 0.0)54icc = float(s2b / (s2b + msw))5556res = {57 "task_id": "F2-conf-01",58 "tau_hat": tau,59 "se_iid": se_iid,60 "se_cluster": se_cl,61 "p_value": p,62 "tau_ci_low": float(lo),63 "tau_ci_high": float(hi),64 "effect_significant_5pct": bool(p < 0.05),65 "icc_hat": icc,66 "n_schools": int(G),67 "n_treated_schools": G1,68 "n_analysed": int(len(d)),69 "feasibility": "complete",70 "feasibility_reason": None,71 "data_issues": issues,72 "n_input_rows": {"schools.csv": int(n_input)},73 "random_seed": None,74 "_diagnostics": {75 "t_stat_cluster": float(tstat), "df_cluster": int(df),76 "p_value_iid": float(ols.pvalues["treated"]),77 "se_school_means_ttest": se_means, "p_school_means_ttest": float(tt.pvalue),78 "design_effect_implied": float((se_cl / se_iid)**2),79 "design_effect_icc": float(1 + (m - 1) * icc),80 "students_per_school": float(m),81 "msb": float(msb), "msw": float(msw),82 },83}84with open("out/results.json", "w") as f:85 json.dump(res, f, indent=2)86print(json.dumps(res, indent=2))87