1"""F2-conf-01: school-level programme effect with cluster-robust inference.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 stats78d = pd.read_csv("data/schools.csv")9n_rows = len(d)10issues = []11# ---- data checks ----12assert d[["school_id", "student_id", "treated", "score"]].notna().all().all(), "missing values"13assert not d.student_id.duplicated().any(), "duplicate student ids"14assert set(d.treated.unique()) <= {0, 1}15assert (d.groupby("school_id").treated.nunique() == 1).all(), "treatment varies within school"16sizes = d.groupby("school_id").size()17assert sizes.nunique() == 1, "unbalanced"1819G = d.school_id.nunique()20G1 = int(d.groupby("school_id").treated.first().sum())21N = len(d)2223# ---- OLS ----24m_iid = smf.ols("score ~ treated", d).fit()25tau = float(m_iid.params["treated"])26se_iid = float(m_iid.bse["treated"])27# CR1 cluster-robust (Liang-Zeger with small-sample factor G/(G-1)*(N-1)/(N-K)),28# inference with t distribution on G-1 df29m_cl = smf.ols("score ~ treated", d).fit(cov_type="cluster",30 cov_kwds={"groups": d.school_id, "use_correction": True, "df_correction": True})31se_cl = float(m_cl.bse["treated"])32df = G - 133tstat = tau / se_cl34p = float(2 * stats.t.sf(abs(tstat), df))35q = stats.t.ppf(0.975, df)36lo, hi = tau - q * se_cl, tau + q * se_cl3738# ---- cross-check: two-sample t-test on school means (balanced design) ----39sm = d.groupby("school_id").agg(y=("score", "mean"), t=("treated", "first"))40tt = stats.ttest_ind(sm.y[sm.t == 1], sm.y[sm.t == 0], equal_var=False)41se_means = float(np.sqrt(sm.y[sm.t == 1].var() / G1 + sm.y[sm.t == 0].var() / (G - G1)))4243# ---- ICC of residual (one-way ANOVA estimator, balanced clusters) ----44d["resid"] = m_iid.resid45m = int(sizes.iloc[0])46grp = d.groupby("school_id").resid47msb = m * ((grp.mean() - d.resid.mean()) ** 2).sum() / (G - 1)48msw = ((d.resid - grp.transform("mean")) ** 2).sum() / (N - G)49s2b = max((msb - msw) / m, 0.0)50icc = float(s2b / (s2b + msw))51deff = 1 + (m - 1) * icc5253res = {54 "task_id": "F2-conf-01",55 "tau_hat": tau,56 "se_iid": se_iid,57 "se_cluster": se_cl,58 "p_value": p,59 "tau_ci_low": float(lo),60 "tau_ci_high": float(hi),61 "effect_significant_5pct": bool(p < 0.05),62 "icc_hat": icc,63 "n_schools": int(G),64 "n_treated_schools": G1,65 "n_analysed": int(N),66 "feasibility": "complete",67 "feasibility_reason": None,68 "data_issues": issues,69 "n_input_rows": {"schools.csv": int(n_rows)},70 "random_seed": None,71 "diagnostics": {72 "cluster_se_type": "CR1 (statsmodels cluster, G/(G-1)*(N-1)/(N-K) correction)",73 "inference_df": df,74 "t_stat_cluster": float(tstat),75 "t_stat_iid": float(tau / se_iid),76 "p_value_iid": float(m_iid.pvalues["treated"]),77 "school_means_welch_t": float(tt.statistic),78 "school_means_welch_p": float(tt.pvalue),79 "se_school_means": se_means,80 "students_per_school": m,81 "design_effect_1_plus_(m-1)icc": float(deff),82 "se_ratio_cluster_over_iid": se_cl / se_iid,83 },84}85os.makedirs("out", exist_ok=True)86with open("out/results.json", "w") as f:87 json.dump(res, f, indent=2)88print(json.dumps(res, indent=2))89