1"""F2-conf-01: school-level programme effect with cluster-robust uncertainty.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/schools.csv; writes out/results.json and out/exclusions.csv. No randomness."""4import json, os5import numpy as np, pandas as pd6import statsmodels.formula.api as smf7from scipy import stats89os.makedirs("out", exist_ok=True)10raw = pd.read_csv("data/schools.csv")11n_raw = len(raw)12issues, excl = [], []13d = raw.copy()1415def step(name, rule, mask, why):16 global d17 nb = len(d); d = d[mask(d)].copy(); na = len(d)18 excl.append(dict(step=name, rule=rule, n_before=nb, n_after=na, n_removed=nb - na, rationale=why))1920step("1", "drop rows with missing school_id/treated/score", lambda x: x[["school_id", "treated", "score"]].notna().all(axis=1), "cannot be used in OLS")21step("2", "drop exact duplicate rows", lambda x: ~x.duplicated(), "duplicate records")22step("3", "drop duplicated student_id (keep first)", lambda x: ~x.student_id.duplicated(), "student_id must be unique")23step("4", "treated must be 0/1", lambda x: x.treated.isin([0, 1]), "binary treatment")24# within-school treatment consistency25bad = d.groupby("school_id").treated.nunique()26bad = bad[bad > 1].index27step("5", "drop schools with non-constant treated", lambda x: ~x.school_id.isin(bad), "treatment must be constant within school")28for e in excl:29 if e["n_removed"] > 0:30 issues.append({"issue": f"{e['rule']}: {e['n_removed']} rows", "handling": "excluded"})31pd.DataFrame(excl).to_csv("out/exclusions.csv", index=False)3233sizes = d.groupby("school_id").size()34if sizes.nunique() != 1:35 issues.append({"issue": "unequal school sizes", "handling": "kept; student-weighted OLS as defined"})3637G = d.school_id.nunique()38G1 = int(d.groupby("school_id").treated.first().sum())39N = len(d)4041# Baseline: conventional iid OLS42m_iid = smf.ols("score ~ treated", data=d).fit()43tau = float(m_iid.params["treated"]); se_iid = float(m_iid.bse["treated"])4445# Main: cluster-robust (CR1, Liang-Zeger with small-sample factor G/(G-1)*(N-1)/(N-K)), t with G-1 df46m_cl = smf.ols("score ~ treated", data=d).fit(cov_type="cluster", cov_kwds={"groups": d.school_id}, use_t=True)47se_cl = float(m_cl.bse["treated"])48df_t = G - 149tstat = tau / se_cl50p = float(2 * stats.t.sf(abs(tstat), df_t))51q = stats.t.ppf(0.975, df_t)52lo, hi = tau - q * se_cl, tau + q * se_cl5354# Sensitivity: school-means Welch/pooled t-test (one obs per school)55sm = d.groupby("school_id").agg(y=("score", "mean"), t=("treated", "first"))56tt = stats.ttest_ind(sm.y[sm.t == 1], sm.y[sm.t == 0], equal_var=True)57se_means = float((sm.y[sm.t == 1].mean() - sm.y[sm.t == 0].mean()) / tt.statistic)5859# ICC of residual (one-way ANOVA estimator, school as random factor)60d["res"] = m_iid.resid61k = N / G if sizes.nunique() == 1 else (N - (sizes**2).sum() / N) / (G - 1)62gm = d.groupby("school_id").res.agg(["mean", "size"])63ssb = float((gm["size"] * (gm["mean"] - d.res.mean())**2).sum())64ssw = float(((d.res - d.school_id.map(gm["mean"]))**2).sum())65msb = ssb / (G - 2) # one df lost to treatment effect (between-school)66msw = ssw / (N - G)67s2b = max((msb - msw) / k, 0.0)68icc = s2b / (s2b + msw)69deff = 1 + (k - 1) * icc7071res = {72 "task_id": "F2-conf-01",73 "tau_hat": tau,74 "se_iid": se_iid,75 "se_cluster": se_cl,76 "p_value": p,77 "tau_ci_low": float(lo),78 "tau_ci_high": float(hi),79 "effect_significant_5pct": bool(p < 0.05),80 "icc_hat": float(icc),81 "n_schools": int(G),82 "n_treated_schools": G1,83 "n_analysed": int(N),84 "feasibility": "complete",85 "feasibility_reason": None,86 "data_issues": issues,87 "n_input_rows": {"schools.csv": int(n_raw)},88 "random_seed": None,89 # supplementary90 "t_stat_cluster": float(tstat),91 "df_cluster": int(df_t),92 "p_value_iid": float(m_iid.pvalues["treated"]),93 "se_ratio_cluster_over_iid": se_cl / se_iid,94 "design_effect_from_icc": float(deff),95 "school_size": float(k),96 "sens_school_means_se": se_means,97 "sens_school_means_p": float(tt.pvalue),98 "mean_treated": float(d.score[d.treated == 1].mean()),99 "mean_control": float(d.score[d.treated == 0].mean()),100}101with open("out/results.json", "w") as f:102 json.dump(res, f, indent=2)103print(json.dumps(res, indent=2))104