1"""F2-conf-01: school-level programme effect with cluster-robust uncertainty."""2import json, os3import numpy as np, pandas as pd4import statsmodels.formula.api as smf5from scipy import stats67ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))8d = pd.read_csv(os.path.join(ROOT, "data", "schools.csv"))9n_rows = len(d)10issues = []11# Data checks12if d.isna().any().any():13 issues.append({"issue": "missing values", "handling": "rows dropped"}); d = d.dropna()14if d.student_id.duplicated().any():15 issues.append({"issue": "duplicate student_id", "handling": "first kept"}); d = d.drop_duplicates("student_id")16g = d.groupby("school_id")17assert (g.treated.nunique() == 1).all(), "treatment varies within school"18assert set(d.treated.unique()) <= {0, 1}19sizes = g.size()2021ols = smf.ols("score ~ treated", data=d).fit()22tau = float(ols.params["treated"]); se_iid = float(ols.bse["treated"])23# Cluster-robust (CR1, Stata-style small-sample correction), inference with t(G-1)24cl = smf.ols("score ~ treated", data=d).fit(cov_type="cluster",25 cov_kwds={"groups": d.school_id, "use_correction": True, "df_correction": True})26se_cl = float(cl.bse["treated"])27G = d.school_id.nunique(); dfree = G - 128tstat = tau / se_cl29p = float(2 * stats.t.sf(abs(tstat), dfree))30q = stats.t.ppf(0.975, dfree)31lo, hi = tau - q * se_cl, tau + q * se_cl3233# Cross-check: school-means two-sample t-test (exact analogue in balanced design)34m = g.agg(score=("score", "mean"), treated=("treated", "first"))35tt = stats.ttest_ind(m.score[m.treated == 1], m.score[m.treated == 0])3637# ICC of residuals: one-way ANOVA estimator on residuals, balanced clusters38d = d.assign(res=ols.resid)39k = int(sizes.iloc[0])40grp_mean = d.groupby("school_id").res.transform("mean")41msb = k * ((d.groupby("school_id").res.mean() - d.res.mean()) ** 2).sum() / (G - 2) # 2 params removed42msw = ((d.res - grp_mean) ** 2).sum() / (len(d) - G)43s2b = max((msb - msw) / k, 0.0)44icc = float(s2b / (s2b + msw))45deff = 1 + (k - 1) * icc4647res = {48 "task_id": "F2-conf-01",49 "tau_hat": tau, "se_iid": se_iid, "se_cluster": se_cl, "p_value": p,50 "tau_ci_low": float(lo), "tau_ci_high": float(hi),51 "effect_significant_5pct": bool(p < 0.05),52 "icc_hat": icc,53 "n_schools": int(G), "n_treated_schools": int(m.treated.sum()), "n_analysed": int(len(d)),54 "feasibility": "complete", "feasibility_reason": None,55 "data_issues": issues, "n_input_rows": {"schools.csv": int(n_rows)},56 "random_seed": None,57 "_diagnostics": {58 "cluster_se_type": "CR1 (G/(G-1)*(N-1)/(N-K)), t with G-1 df",59 "t_stat_cluster": float(tstat), "df": dfree,60 "school_means_ttest_t": float(tt.statistic), "school_means_ttest_p": float(tt.pvalue),61 "se_ratio_cluster_over_iid": se_cl / se_iid,62 "design_effect_from_icc": deff, "sqrt_design_effect": float(np.sqrt(deff)),63 "students_per_school": k, "p_value_iid": float(ols.pvalues["treated"]),64 },65}66os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)67with open(os.path.join(ROOT, "out", "results.json"), "w") as f:68 json.dump(res, f, indent=2)69print(json.dumps(res, indent=2))70