1"""F2-conf-01: school-level programme effect with cluster-robust uncertainty.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Deterministic; no randomness used."""4import json, os, csv5import numpy as np, pandas as pd6import statsmodels.api as sm7from scipy import stats89D = "data/schools.csv"10d = pd.read_csv(D)11n_raw = len(d)12issues = []13ledger = [("0", "raw input", n_raw, n_raw, 0, "all rows read")]1415# Audit-driven filters (none expected to remove rows; recorded for the ledger)16n0 = len(d); d = d.dropna(subset=["school_id", "treated", "score"])17ledger.append(("1", "drop rows with missing school_id/treated/score", n0, len(d), n0 - len(d), "cannot be used in OLS"))18n0 = len(d); d = d.drop_duplicates(subset=["student_id"])19ledger.append(("2", "drop duplicate student_id", n0, len(d), n0 - len(d), "student_id must be unique"))20if (d.groupby("school_id")["treated"].nunique() > 1).any():21 issues.append({"issue": "treated varies within school", "handling": "reported; not altered"})2223y = d["score"].to_numpy(float); t = d["treated"].to_numpy(float)24X = sm.add_constant(t)25ols = sm.OLS(y, X).fit()26tau = float(ols.params[1]); se_iid = float(ols.bse[1])2728# Cluster-robust (CR1, Stata-style small-sample correction), inference with t(G-1)29cl = sm.OLS(y, X).fit(cov_type="cluster", cov_kwds={"groups": d["school_id"].to_numpy(), "use_correction": True, "df_correction": True})30se_cl = float(cl.bse[1])31G = int(d["school_id"].nunique()); dfree = G - 132tstat = tau / se_cl33p = float(2 * stats.t.sf(abs(tstat), dfree))34q = stats.t.ppf(0.975, dfree)35lo, hi = tau - q * se_cl, tau + q * se_cl3637# Cross-check: two-sample t-test on school means (Welch & pooled, df=G-2)38sm_ = d.groupby("school_id").agg(m=("score", "mean"), tr=("treated", "first"))39a, b = sm_.loc[sm_.tr == 1, "m"], sm_.loc[sm_.tr == 0, "m"]40pooled = stats.ttest_ind(a, b, equal_var=True)41se_means = float(np.sqrt(a.var(ddof=1) / len(a) + b.var(ddof=1) / len(b)))4243# ICC of residual: one-way ANOVA estimator on OLS residuals by school44r = pd.Series(ols.resid, index=d.index); gid = d["school_id"]45m = d.groupby("school_id").size(); n_bar = float(m.mean())46gm = r.groupby(gid).transform("mean")47ssw = float(((r - gm) ** 2).sum()); ssb = float((m * r.groupby(gid).mean() ** 2).sum())48msw = ssw / (len(d) - G); msb = ssb / (G - 1)49# note: residual df within schools unaffected since treated is constant within school50s2b = max((msb - msw) / n_bar, 0.0)51icc = s2b / (s2b + msw)52deff = 1 + (n_bar - 1) * icc5354res = {55 "task_id": "F2-conf-01",56 "tau_hat": tau, "se_iid": se_iid, "se_cluster": se_cl,57 "p_value": p, "tau_ci_low": float(lo), "tau_ci_high": float(hi),58 "effect_significant_5pct": bool(p < 0.05),59 "icc_hat": float(icc),60 "n_schools": G, "n_treated_schools": int(sm_.tr.sum()), "n_analysed": int(len(d)),61 "feasibility": "complete", "feasibility_reason": None,62 "data_issues": issues, "n_input_rows": {"schools.csv": int(n_raw)},63 "random_seed": None,64 "cluster_df": dfree, "t_stat_cluster": float(tstat),65 "se_ratio_cluster_to_iid": se_cl / se_iid,66 "design_effect_from_icc": float(deff),67 "check_se_school_means": se_means,68 "check_p_school_means_pooled_t": float(pooled.pvalue),69 "p_value_iid": float(ols.pvalues[1]),70 "students_per_school": n_bar,71}72os.makedirs("out", exist_ok=True)73with open("out/results.json", "w") as f:74 json.dump(res, f, indent=2)75with open("out/exclusions.csv", "w", newline="") as f:76 w = csv.writer(f); w.writerow(["step", "rule", "n_before", "n_after", "n_removed", "rationale"]); w.writerows(ledger)77print(json.dumps(res, indent=2))78