1"""F2-conf-01: school-level programme effect with cluster-robust inference.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 np6import pandas as pd7import statsmodels.formula.api as smf8from scipy import stats910ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))11df = pd.read_csv(os.path.join(ROOT, "data", "schools.csv"))12n_raw = len(df)1314issues = []15excl = []16# Step 1: rows with missing values in required columns17n0 = len(df)18df = df.dropna(subset=["school_id", "student_id", "treated", "score"])19excl.append(("1", "drop rows with missing school_id/student_id/treated/score", n0, len(df), n0 - len(df),20 "required for estimation"))21# Step 2: duplicate student ids22n0 = len(df)23df = df.drop_duplicates(subset=["student_id"])24excl.append(("2", "drop duplicate student_id", n0, len(df), n0 - len(df), "student_id must be unique"))25for s, r, a, b, k, why in excl:26 if k > 0:27 issues.append({"issue": f"{k} rows removed: {r}", "handling": "excluded"})28# Design checks29within = df.groupby("school_id")["treated"].nunique()30if (within > 1).any():31 issues.append({"issue": f"treated varies within {int((within>1).sum())} schools", "handling": "reported; not altered"})32sizes = df.groupby("school_id").size()33if sizes.nunique() != 1:34 issues.append({"issue": "unbalanced school sizes", "handling": "kept; CR1 SE handles unequal clusters"})3536pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(37 os.path.join(ROOT, "out", "exclusions.csv"), index=False)3839G = df["school_id"].nunique()40G1 = int(df.groupby("school_id")["treated"].first().sum())41N = len(df)4243# Baseline: naive OLS with iid homoskedastic SE44ols = smf.ols("score ~ treated", data=df).fit()45tau = float(ols.params["treated"])46se_iid = float(ols.bse["treated"])4748# Main: cluster-robust (CR1, Stata-style small-sample correction) by school, t with G-1 df49cl = smf.ols("score ~ treated", data=df).fit(cov_type="cluster", cov_kwds={"groups": df["school_id"]},50 use_t=True)51se_cl = float(cl.bse["treated"])52dof = G - 153tcrit = float(stats.t.ppf(0.975, dof))54tstat = tau / se_cl55p = float(2 * stats.t.sf(abs(tstat), dof))56lo, hi = tau - tcrit * se_cl, tau + tcrit * se_cl5758# Cross-check: two-sample Welch t-test on school means (aggregated analysis)59sm = df.groupby("school_id").agg(m=("score", "mean"), t=("treated", "first"))60a, b = sm.m[sm.t == 1], sm.m[sm.t == 0]61agg = stats.ttest_ind(a, b, equal_var=False)62se_agg = float(np.sqrt(a.var(ddof=1) / len(a) + b.var(ddof=1) / len(b)))6364# ICC of residual (after removing arm means), one-way ANOVA estimator65df = df.assign(res=df["score"] - df.groupby("treated")["score"].transform("mean"))66grp = df.groupby("school_id")["res"]67ng = grp.size(); mg = grp.mean()68ssw = float(((df["res"] - grp.transform("mean")) ** 2).sum())69ssb = float((ng * (mg - df["res"].mean()) ** 2).sum())70msw = ssw / (N - G)71msb = ssb / (G - 2) # two df used by arm means72n0_eff = (N - (ng ** 2).sum() / N) / (G - 1)73sb2 = max((msb - msw) / n0_eff, 0.0)74icc = sb2 / (sb2 + msw)75deff = 1 + (float(ng.mean()) - 1) * icc7677res = {78 "task_id": "F2-conf-01",79 "tau_hat": tau,80 "se_iid": se_iid,81 "se_cluster": se_cl,82 "p_value": p,83 "tau_ci_low": float(lo),84 "tau_ci_high": float(hi),85 "effect_significant_5pct": bool(p < 0.05),86 "icc_hat": float(icc),87 "n_schools": int(G),88 "n_treated_schools": G1,89 "n_analysed": int(N),90 "feasibility": "complete",91 "feasibility_reason": None,92 "data_issues": issues,93 "n_input_rows": {"schools.csv": int(n_raw)},94 "random_seed": None,95 # supplementary96 "baseline_p_value_iid": float(ols.pvalues["treated"]),97 "baseline_ci_iid_low": float(ols.conf_int().loc["treated", 0]),98 "baseline_ci_iid_high": float(ols.conf_int().loc["treated", 1]),99 "t_stat_cluster": float(tstat),100 "df_cluster": int(dof),101 "se_ratio_cluster_to_iid": se_cl / se_iid,102 "design_effect_from_icc": float(deff),103 "school_means_welch_se": se_agg,104 "school_means_welch_p": float(agg.pvalue),105 "mean_score_treated": float(df.loc[df.treated == 1, "score"].mean()),106 "mean_score_control": float(df.loc[df.treated == 0, "score"].mean()),107}108with open(os.path.join(ROOT, "out", "results.json"), "w") as f:109 json.dump(res, f, indent=2)110print(json.dumps(res, indent=2))111