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, os5import numpy as np, pandas as pd6import statsmodels.formula.api as smf7from scipy import stats89ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))10d = pd.read_csv(os.path.join(ROOT, "data", "schools.csv"))11n_raw = len(d)12issues = []13excl = []1415def step(df, name, mask, rule, rationale):16 nb = len(df); df2 = df[mask].copy()17 excl.append(dict(step=name, rule=rule, n_before=nb, n_after=len(df2),18 n_removed=nb - len(df2), rationale=rationale))19 return df22021d = step(d, "1", d[["school_id", "student_id", "treated", "score"]].notna().all(axis=1),22 "drop rows with any missing value", "outcome/treatment/cluster id required")23d = step(d, "2", ~d.student_id.duplicated(), "drop duplicate student_id", "one row per student")24d = step(d, "3", d.treated.isin([0, 1]), "treated must be 0/1", "binary assignment")25tconst = d.groupby("school_id").treated.transform("nunique") == 126d = step(d, "4", tconst, "treated constant within school", "school-level assignment")27for e in excl:28 if e["n_removed"] > 0:29 issues.append({"issue": f"{e['rule']}: {e['n_removed']} rows", "handling": "excluded"})30sizes = d.groupby("school_id").size()31if sizes.nunique() != 1:32 issues.append({"issue": "unequal school sizes", "handling": "kept; CR1 handles unbalanced clusters"})3334# Baseline: OLS with conventional iid/homoskedastic SE35ols = smf.ols("score ~ treated", data=d).fit()36tau = float(ols.params["treated"]); se_iid = float(ols.bse["treated"])3738# Main: cluster-robust CR1 SE (small-sample factor G/(G-1)*(N-1)/(N-K)), clustered by school,39# inference with t distribution on G-1 degrees of freedom40G = int(d.school_id.nunique())41cl = smf.ols("score ~ treated", data=d).fit(42 cov_type="cluster", cov_kwds={"groups": d.school_id.values, "use_correction": True, "df_correction": True})43se_cl = float(cl.bse["treated"])44dfree = G - 145tstat = tau / se_cl46p = float(2 * stats.t.sf(abs(tstat), dfree))47crit = float(stats.t.ppf(0.975, dfree))48lo, hi = tau - crit * se_cl, tau + crit * se_cl4950# Cross-check: Welch t-test on school means (equal school sizes => identical point estimate)51sm = d.groupby("school_id").agg(score=("score", "mean"), treated=("treated", "first"))52a, b = sm.score[sm.treated == 1], sm.score[sm.treated == 0]53tt = stats.ttest_ind(a, b, equal_var=False)54se_means = float((a.mean() - b.mean()) / tt.statistic)5556# ICC of residual (score minus arm mean): one-way ANOVA estimator, schools as groups57d["resid"] = ols.resid58m = float(sizes.mean())59grp = d.groupby("school_id").resid60msb = float((sizes * grp.mean() ** 2).sum() / (G - 2)) # 2 arm means removed61msw = float(((d.resid - grp.transform("mean")) ** 2).sum() / (len(d) - G))62s2b = max((msb - msw) / m, 0.0)63icc = s2b / (s2b + msw)64deff = 1 + (m - 1) * icc6566res = {67 "task_id": "F2-conf-01",68 "tau_hat": tau, "se_iid": se_iid, "se_cluster": se_cl, "p_value": p,69 "tau_ci_low": float(lo), "tau_ci_high": float(hi),70 "effect_significant_5pct": bool(p < 0.05),71 "icc_hat": float(icc),72 "n_schools": G, "n_treated_schools": int(sm.treated.sum()), "n_analysed": int(len(d)),73 "feasibility": "complete", "feasibility_reason": None,74 "data_issues": issues, "n_input_rows": {"schools.csv": int(n_raw)}, "random_seed": None,75 # supplementary keys76 "t_stat_cluster": float(tstat), "df_cluster": dfree, "t_crit_975": crit,77 "se_ratio_cluster_over_iid": se_cl / se_iid,78 "p_value_iid": float(ols.pvalues["treated"]),79 "se_school_means_welch": se_means, "p_school_means_welch": float(tt.pvalue),80 "students_per_school": m, "design_effect_from_icc": float(deff),81 "sqrt_design_effect": float(np.sqrt(deff)),82 "sigma2_between_resid": float(s2b), "sigma2_within_resid": msw,83 "mean_score_treated": float(d.score[d.treated == 1].mean()),84 "mean_score_control": float(d.score[d.treated == 0].mean()),85}86with open(os.path.join(ROOT, "out", "results.json"), "w") as f:87 json.dump(res, f, indent=2)88pd.DataFrame(excl).to_csv(os.path.join(ROOT, "out", "exclusions.csv"), index=False)89print(json.dumps(res, indent=2))90