1"""F2-conf-01: school-level programme effect with design-valid (cluster) uncertainty.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, hashlib4import numpy as np, pandas as pd5import statsmodels.formula.api as smf6from scipy import stats78DATA = "data/schools.csv"9raw = pd.read_csv(DATA)10n_raw = len(raw)11issues, excl = [], []1213# Audit-based filters (each logged even if it removes nothing)14d = raw.copy()15def step(name, rule, mask, why):16 global d17 nb = len(d); d = d[mask(d)].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", lambda x: x[["school_id", "treated", "score"]].notna().all(axis=1), "estimand requires all three")20step("2", "drop duplicated student_id", lambda x: ~x.student_id.duplicated(), "one row per student")21step("3", "drop schools with non-constant treated", lambda x: x.groupby("school_id").treated.transform("nunique") == 1, "treatment must be school-level")22pd.DataFrame(excl).to_csv("out/exclusions.csv", index=False)2324sizes = d.groupby("school_id").size()25if sizes.nunique() != 1:26 issues.append({"issue": "unbalanced school sizes", "handling": "kept; cluster-robust SE does not require balance"})2728G = d.school_id.nunique()29G1 = int(d.groupby("school_id").treated.first().sum())30N = len(d)3132# Naive baseline: conventional OLS SE33ols = smf.ols("score ~ treated", data=d).fit()34tau = float(ols.params["treated"]); se_iid = float(ols.bse["treated"])3536# Main: CR1 cluster-robust SE by school, t reference with G-1 df37cl = smf.ols("score ~ treated", data=d).fit(cov_type="cluster", cov_kwds={"groups": d.school_id.values}, use_t=True)38se_cl = float(cl.bse["treated"])39df_t = G - 140tstat = tau / se_cl41p = float(2 * stats.t.sf(abs(tstat), df_t))42q = stats.t.ppf(0.975, df_t)43lo, hi = tau - q * se_cl, tau + q * se_cl4445# Robustness: OLS on school means (exact in balanced design), t with G-2 df46m = d.groupby("school_id").agg(score=("score", "mean"), treated=("treated", "first")).reset_index()47agg = smf.ols("score ~ treated", data=m).fit()48se_agg = float(agg.bse["treated"]); p_agg = float(agg.pvalues["treated"])4950# ICC of residual (ANOVA estimator, treatment effect removed)51d["resid"] = ols.resid52msw = ((d.resid - d.groupby("school_id").resid.transform("mean")) ** 2).sum() / (N - G)53mbar = N / G54gm = d.groupby("school_id").resid.mean()55msb = (sizes * (gm - d.resid.mean()) ** 2).sum() / (G - 2) # 2 df used by intercept+treatment56s2b = max((msb - msw) / mbar, 0.0)57icc = float(s2b / (s2b + msw))58deff = 1 + (mbar - 1) * icc5960res = {61 "task_id": "F2-conf-01",62 "tau_hat": tau, "se_iid": se_iid, "se_cluster": se_cl, "p_value": p,63 "tau_ci_low": float(lo), "tau_ci_high": float(hi),64 "effect_significant_5pct": bool(p < 0.05),65 "icc_hat": icc, "n_schools": int(G), "n_treated_schools": G1, "n_analysed": int(N),66 "feasibility": "complete", "feasibility_reason": None,67 "data_issues": issues, "n_input_rows": {"schools.csv": int(n_raw)}, "random_seed": None,68 # supplementary69 "cluster_se_type": "CR1 (statsmodels), t reference with G-1 df",70 "t_df": int(df_t), "t_stat_cluster": float(tstat),71 "se_ratio_cluster_to_iid": se_cl / se_iid, "design_effect_from_icc": float(deff),72 "school_means_se": se_agg, "school_means_p_value": p_agg,73 "mean_control": float(d.loc[d.treated == 0, "score"].mean()),74 "mean_treated": float(d.loc[d.treated == 1, "score"].mean()),75 "students_per_school": float(mbar),76 "input_sha256": hashlib.sha256(open(DATA, "rb").read()).hexdigest(),77}78json.dump(res, open("out/results.json", "w"), indent=2)79print(json.dumps(res, indent=2))80