1"""F2-conf-01: school-level programme effect with cluster-robust uncertainty.2Reads data/schools.csv only; writes out/results.json (deterministic, no randomness)."""3import json, os4import numpy as np, pandas as pd5import statsmodels.formula.api as smf6from scipy import stats78ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))9d = pd.read_csv(os.path.join(ROOT, "data", "schools.csv"))10n_raw = len(d)11issues = []12# Audit checks (no exclusions needed if all pass)13assert d[["school_id", "student_id", "treated", "score"]].notna().all().all(), "missing values"14assert not d.student_id.duplicated().any(), "duplicate student ids"15assert (d.groupby("school_id").treated.nunique() == 1).all(), "treatment varies within school"16assert set(d.treated.unique()) <= {0, 1}17sizes = d.groupby("school_id").size()18if sizes.nunique() != 1:19 issues.append({"issue": "unbalanced school sizes", "handling": "kept; CR1 handles it"})2021# OLS of score on treated with intercept22m_iid = smf.ols("score ~ treated", data=d).fit()23G = d.school_id.nunique()24m_cl = smf.ols("score ~ treated", data=d).fit(cov_type="cluster",25 cov_kwds={"groups": d.school_id}, use_t=True)26tau = float(m_iid.params["treated"])27se_iid = float(m_iid.bse["treated"])28se_cl = float(m_cl.bse["treated"]) # CR1: G/(G-1)*(N-1)/(N-K) small-sample correction29df = G - 130tstat = tau / se_cl31p = float(2 * stats.t.sf(abs(tstat), df))32q = stats.t.ppf(0.975, df)33lo, hi = tau - q * se_cl, tau + q * se_cl3435# Sensitivity: two-sample t-test on school means (balanced design)36sm = d.groupby("school_id").agg(score=("score", "mean"), treated=("treated", "first"))37tt = stats.ttest_ind(sm.score[sm.treated == 1], sm.score[sm.treated == 0], equal_var=True)3839# ICC of residual: one-way ANOVA estimator on OLS residuals by school40d["resid"] = m_iid.resid41m = int(sizes.iloc[0]) if sizes.nunique() == 1 else float(sizes.mean())42gm = d.groupby("school_id").resid.mean()43msb = m * ((gm - d.resid.mean()) ** 2).sum() / (G - 1)44msw = ((d.resid - d.school_id.map(gm)) ** 2).sum() / (len(d) - G)45s2b = max((msb - msw) / m, 0.0)46icc = float(s2b / (s2b + msw))47deff = 1 + (m - 1) * icc4849res = {50 "task_id": "F2-conf-01",51 "tau_hat": tau,52 "se_iid": se_iid,53 "se_cluster": se_cl,54 "p_value": p,55 "tau_ci_low": float(lo),56 "tau_ci_high": float(hi),57 "effect_significant_5pct": bool(p < 0.05),58 "icc_hat": icc,59 "n_schools": int(G),60 "n_treated_schools": int(sm.treated.sum()),61 "n_analysed": int(len(d)),62 "feasibility": "complete",63 "feasibility_reason": None,64 "data_issues": issues,65 "n_input_rows": {"schools.csv": int(n_raw)},66 "random_seed": None,67 "t_df": int(df),68 "t_stat_cluster": float(tstat),69 "p_value_iid": float(m_iid.pvalues["treated"]),70 "se_ratio_cluster_to_iid": se_cl / se_iid,71 "design_effect_from_icc": float(deff),72 "school_mean_ttest_t": float(tt.statistic),73 "school_mean_ttest_p": float(tt.pvalue),74 "students_per_school": m,75}76os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)77with open(os.path.join(ROOT, "out", "results.json"), "w") as f:78 json.dump(res, f, indent=2)79print(json.dumps(res, indent=2))80