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, hashlib, os5import numpy as np, pandas as pd6import statsmodels.formula.api as smf7from scipy import stats89IN = "data/schools.csv"10os.makedirs("out", exist_ok=True)11d = pd.read_csv(IN)12n_raw = len(d)13issues, excl = [], []1415def step(name, rule, mask, why):16 global d17 nb = len(d); d = d[mask].copy(); excl.append([name, rule, nb, len(d), nb - len(d), why])1819# Audit-driven filters (each logged even if it removes nothing)20step("1", "drop rows with missing school_id/treated/score", d[["school_id", "treated", "score"]].notna().all(axis=1), "required for estimation")21step("2", "drop duplicate student_id (keep first)", ~d.student_id.duplicated(), "one row per student")22nonconst = d.groupby("school_id").treated.transform("nunique") > 123if nonconst.any():24 issues.append({"issue": "treated not constant within some schools", "handling": "schools dropped"})25step("3", "drop schools with non-constant treated", ~nonconst, "design: assignment at school level")26pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv("out/exclusions.csv", index=False)2728sizes = d.groupby("school_id").size()29if sizes.nunique() != 1:30 issues.append({"issue": "unequal school sizes", "handling": "kept; CR1 cluster SE is valid regardless"})3132G = d.school_id.nunique()33G1 = int(d.groupby("school_id").treated.first().sum())3435# Baseline: OLS with iid SE36ols = smf.ols("score ~ treated", data=d).fit()37tau = float(ols.params["treated"]); se_iid = float(ols.bse["treated"])3839# Main: CR1 cluster-robust SE by school, t reference with G-1 df40cl = smf.ols("score ~ treated", data=d).fit(cov_type="cluster", cov_kwds={"groups": d.school_id.values, "use_correction": True, "df_correction": True})41se_cl = float(cl.bse["treated"])42df_t = G - 143tstat = tau / se_cl44p = float(2 * stats.t.sf(abs(tstat), df_t))45q = stats.t.ppf(0.975, df_t)46lo, hi = tau - q * se_cl, tau + q * se_cl4748# Diagnostic: two-sample t-test on school means (equivalent aggregated analysis)49sm = d.groupby("school_id").agg(score=("score", "mean"), treated=("treated", "first"))50a, b = sm.score[sm.treated == 1], sm.score[sm.treated == 0]51tt = stats.ttest_ind(a, b, equal_var=True)52se_means = float(np.sqrt(a.var(ddof=1) / len(a) + b.var(ddof=1) / len(b)))5354# ICC of residual (after removing treatment-arm means): one-way ANOVA estimator55d["resid"] = ols.resid56m = sizes.mean() # balanced57grp = d.groupby("school_id").resid58msb = m * ((grp.mean() - d.resid.mean()) ** 2).sum() / (G - 2) # 2 df used by arm means59msw = ((d.resid - grp.transform("mean")) ** 2).sum() / (len(d) - G)60s2b = max((msb - msw) / m, 0.0)61icc = float(s2b / (s2b + msw))62deff = 1 + (m - 1) * icc6364res = {65 "task_id": "F2-conf-01",66 "tau_hat": tau, "se_iid": se_iid, "se_cluster": se_cl, "p_value": p,67 "tau_ci_low": float(lo), "tau_ci_high": float(hi),68 "effect_significant_5pct": bool(p < 0.05),69 "icc_hat": icc,70 "n_schools": int(G), "n_treated_schools": G1, "n_analysed": int(len(d)),71 "feasibility": "complete", "feasibility_reason": None,72 "data_issues": issues,73 "n_input_rows": {"schools.csv": int(n_raw)},74 "random_seed": None,75 "diagnostics": {76 "input_sha256": hashlib.sha256(open(IN, "rb").read()).hexdigest(),77 "cluster_se_type": "CR1 (Stata-style small-sample correction), t reference with G-1 df",78 "t_df": int(df_t), "t_stat_cluster": float(tstat),79 "p_value_iid": float(ols.pvalues["treated"]),80 "se_ratio_cluster_over_iid": se_cl / se_iid,81 "students_per_school": float(m),82 "design_effect_moulton": float(deff), "sqrt_design_effect": float(np.sqrt(deff)),83 "school_means_ttest_t": float(tt.statistic), "school_means_ttest_p_df38": float(tt.pvalue),84 "school_means_se": se_means,85 "msb": float(msb), "msw": float(msw), "sigma2_between": float(s2b),86 },87}88with open("out/results.json", "w") as f:89 json.dump(res, f, indent=2)90print(json.dumps({k: res[k] for k in ["tau_hat", "se_iid", "se_cluster", "p_value", "tau_ci_low", "tau_ci_high", "effect_significant_5pct", "icc_hat"]}, indent=2))91