1"""F2-conf-01: school-level programme effect with cluster-robust uncertainty.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3No randomness is used (random_seed = null)."""4import json, os5import numpy as np, pandas as pd6import statsmodels.formula.api as smf7from scipy import stats89d = pd.read_csv("data/schools.csv")10n_raw = len(d)11issues = []12excl = []1314# Audit / filters (all checks; none remove rows on the supplied data)15def step(name, rule, mask, why):16 global d17 nb = len(d); d = d[mask].copy()18 excl.append(dict(step=name, rule=rule, n_before=nb, n_after=len(d), n_removed=nb - len(d), rationale=why))19step("1", "drop rows with missing school_id/treated/score", d[["school_id", "treated", "score"]].notna().all(axis=1), "cannot be used in OLS")20step("2", "drop duplicate student_id (keep first)", ~d.student_id.duplicated(), "one row per student expected")21tv = d.groupby("school_id").treated.transform("nunique")22step("3", "drop schools with non-constant treated", tv == 1, "treatment must be constant within school")23pd.DataFrame(excl).to_csv("out/exclusions.csv", index=False)24sizes = d.groupby("school_id").size()25if sizes.nunique() != 1:26 issues.append({"issue": "unbalanced school sizes", "handling": "kept; CR1 SE is valid regardless"})2728# Baseline: OLS with iid SE29m_iid = smf.ols("score ~ treated", data=d).fit()30tau = float(m_iid.params["treated"]); se_iid = float(m_iid.bse["treated"])3132# Main: cluster-robust (CR1, statsmodels small-sample correction) clustered by school, t(G-1) reference33m_cl = smf.ols("score ~ treated", data=d).fit(cov_type="cluster", cov_kwds={"groups": d.school_id, "use_correction": True, "df_correction": True})34se_cl = float(m_cl.bse["treated"])35G = int(d.school_id.nunique()); dof = G - 136tstat = tau / se_cl37p = float(2 * stats.t.sf(abs(tstat), dof))38q = stats.t.ppf(0.975, dof)39lo, hi = tau - q * se_cl, tau + q * se_cl4041# Cross-check: two-sample t on school means (balanced design -> equivalent analysis)42sm = d.groupby("school_id").agg(y=("score", "mean"), t=("treated", "first"))43tt = stats.ttest_ind(sm.y[sm.t == 1], sm.y[sm.t == 0], equal_var=False)4445# ICC of residuals: one-way ANOVA estimator46d["res"] = m_iid.resid47k = sizes.mean()48grp = d.groupby("school_id").res49msb = (sizes * (grp.mean() - d.res.mean()) ** 2).sum() / (G - 1)50msw = ((d.res - grp.transform("mean")) ** 2).sum() / (len(d) - G)51s2b = max((msb - msw) / k, 0.0)52icc = float(s2b / (s2b + msw))53deff = 1 + (k - 1) * icc5455res = {56 "task_id": "F2-conf-01",57 "tau_hat": tau, "se_iid": se_iid, "se_cluster": se_cl, "p_value": p,58 "tau_ci_low": float(lo), "tau_ci_high": float(hi),59 "effect_significant_5pct": bool(p < 0.05),60 "icc_hat": icc,61 "n_schools": G, "n_treated_schools": int(sm.t.sum()), "n_analysed": int(len(d)),62 "feasibility": "complete", "feasibility_reason": None,63 "data_issues": issues, "n_input_rows": {"schools.csv": n_raw}, "random_seed": None,64 "aux_t_stat_cluster": float(tstat), "aux_df_cluster": dof,65 "aux_p_value_iid": float(m_iid.pvalues["treated"]),66 "aux_se_ratio_cluster_over_iid": se_cl / se_iid,67 "aux_design_effect_moulton": float(deff), "aux_students_per_school": float(k),68 "aux_school_means_welch_t": float(tt.statistic), "aux_school_means_welch_p": float(tt.pvalue),69}70with open("out/results.json", "w") as f:71 json.dump(res, f, indent=2)72print(json.dumps(res, indent=2))73