1"""F2-conf-01: school-level programme effect with cluster-robust uncertainty.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np, pandas as pd5import statsmodels.formula.api as smf6from scipy import stats78ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))9OUT = os.path.join(ROOT, "out")10df_raw = pd.read_csv(os.path.join(ROOT, "data", "schools.csv"))11n_raw = len(df_raw)1213# ---- audit / exclusions ----14issues, excl = [], []15def step(name, rule, before, after, why):16 excl.append(dict(step=name, rule=rule, n_before=before, n_after=after, n_removed=before - after, rationale=why))17d = df_raw.dropna(subset=["school_id", "treated", "score"])18step("1", "drop rows with missing school_id/treated/score", n_raw, len(d), "required for estimation")19b = len(d); d = d.drop_duplicates(subset=["student_id"])20step("2", "drop duplicated student_id", b, len(d), "one row per student")21if (d.groupby("school_id")["treated"].nunique() > 1).any():22 issues.append({"issue": "treated varies within school", "handling": "flagged; not expected by design"})23sizes = d.groupby("school_id").size()24if sizes.nunique() != 1:25 issues.append({"issue": "unbalanced school sizes", "handling": "kept; student-weighted estimand unaffected"})26pd.DataFrame(excl).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)2728# ---- estimation ----29N = len(d)30G = d["school_id"].nunique()31G1 = int(d.groupby("school_id")["treated"].first().sum())32ols = smf.ols("score ~ treated", data=d).fit()33tau = float(ols.params["treated"]); se_iid = float(ols.bse["treated"])34# CR1 cluster-robust by school (small-sample factor G/(G-1)*(N-1)/(N-K)), t(G-1) reference35cl = smf.ols("score ~ treated", data=d).fit(cov_type="cluster", cov_kwds={"groups": d["school_id"]},36 use_t=True, df_correction=True)37se_cl = float(cl.bse["treated"])38dof = G - 139p = float(2 * stats.t.sf(abs(tau / se_cl), dof))40q = stats.t.ppf(0.975, dof)41lo, hi = tau - q * se_cl, tau + q * se_cl4243# ---- ICC of OLS residual (one-way ANOVA estimator; residual school means have G-2 df) ----44d = d.assign(resid=ols.resid)45grp = d.groupby("school_id")["resid"]46m = float(sizes.mean()) # balanced design -> common cluster size47msb = (grp.size() * grp.mean() ** 2).sum() / (G - 2)48msw = ((d["resid"] - grp.transform("mean")) ** 2).sum() / (N - G)49s2b = max((msb - msw) / m, 0.0)50icc = float(s2b / (s2b + msw))5152# ---- robustness: two-sample t-test on school means (G-2 df) ----53sm = d.groupby("school_id").agg(score=("score", "mean"), treated=("treated", "first"))54a, c = sm.score[sm.treated == 1], sm.score[sm.treated == 0]55tt = stats.ttest_ind(a, c, equal_var=True)56se_sm = float((a.mean() - c.mean()) / tt.statistic)5758res = {59 "task_id": "F2-conf-01",60 "tau_hat": tau, "se_iid": se_iid, "se_cluster": se_cl, "p_value": p,61 "tau_ci_low": float(lo), "tau_ci_high": float(hi),62 "effect_significant_5pct": bool(p < 0.05),63 "icc_hat": icc, "n_schools": int(G), "n_treated_schools": G1, "n_analysed": int(N),64 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,65 "n_input_rows": {"schools.csv": int(n_raw)}, "random_seed": None,66 "baseline_p_value_iid": float(ols.pvalues["treated"]),67 "design_effect_se_ratio": se_cl / se_iid,68 "cluster_t_df": int(dof),69 "robust_school_mean_ttest_se": se_sm,70 "robust_school_mean_ttest_p": float(tt.pvalue),71 "mean_score_treated": float(d.score[d.treated == 1].mean()),72 "mean_score_control": float(d.score[d.treated == 0].mean()),73 "students_per_school": int(m),74}75with open(os.path.join(OUT, "results.json"), "w") as f:76 json.dump(res, f, indent=2)77print(json.dumps(res, indent=2))78