1"""T2-simulation: finite-sample coverage of CR0, CR1 and wild-cluster-bootstrap-t CIs for a mean.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/design.csv; writes out/results.json, out/table1.csv, out/figure1.png.4RNG: per cell numpy.random.default_rng(seed from design.csv). Within each replication, draws are5taken in the order: Z_u (G std normals), Z_e (N std normals), Rademacher weights (B x G, via rng.integers).6"""7import json, os8import numpy as np, pandas as pd9from scipy import stats10import matplotlib11matplotlib.use("Agg")12import matplotlib.pyplot as plt1314ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))15OUT = os.path.join(ROOT, "out")16design = pd.read_csv(os.path.join(ROOT, "data", "design.csv"))171819def key(est, G, icc):20 return f"{est}_G{int(G)}_icc{float(icc):.1f}"212223def run_cell(row):24 G, icc, N, R, B = int(row.G), float(row.icc), int(row.N), int(row.R), int(row.B)25 mu, level = float(row.mu), float(row.nominal_level)26 n = np.array([int(x) for x in row.cluster_sizes.split("|")])27 assert len(n) == G and n.sum() == N28 gid = np.repeat(np.arange(G), n)29 var_exact = (icc * np.sum(n.astype(float) ** 2) + (1 - icc) * N) / N ** 230 z = 1.96 # as specified in TASK.md31 tcrit = stats.t.ppf(0.5 + level / 2, G - 1)32 adj = G / (G - 1)33 rng = np.random.default_rng(int(row.seed))34 res = {k: np.empty(R) for k in ["ybar", "v0", "v1", "q", "s2"]}35 for r in range(R):36 u = np.sqrt(icc) * rng.standard_normal(G)37 e = np.sqrt(1 - icc) * rng.standard_normal(N)38 y = mu + u[gid] + e39 ybar = y.mean()40 ehat = y - ybar41 S = np.bincount(gid, weights=ehat, minlength=G) # cluster residual sums42 v0 = np.sum(S ** 2) / N ** 243 v1 = adj * v044 # wild cluster bootstrap-t (unrestricted, Rademacher), vectorised over B45 W = 2.0 * rng.integers(0, 2, size=(B, G)) - 1.0 # B x G46 d = (W @ S) / N # ybar* - ybar, length B47 # bootstrap residual cluster sums: sum_i (w_g ehat_ig - d) = w_g S_g - n_g d48 Sstar = W * S[None, :] - d[:, None] * n[None, :]49 vstar = adj * np.sum(Sstar ** 2, axis=1) / N ** 250 tstar = np.abs(d) / np.sqrt(vstar)51 q = np.quantile(tstar, 0.95) # numpy default (linear interpolation)52 res["ybar"][r], res["v0"][r], res["v1"][r], res["q"][r] = ybar, v0, v1, q53 res["s2"][r] = ehat @ ehat / (N - 1)54 hw = {"cr0": z * np.sqrt(res["v0"]),55 "cr1": tcrit * np.sqrt(res["v1"]),56 "wcb": res["q"] * np.sqrt(res["v1"]),57 "iid": z * np.sqrt(res["s2"] / N)} # baseline: naive i.i.d. interval58 out = {}59 for est, h in hw.items():60 cov = float(np.mean(np.abs(res["ybar"] - mu) <= h))61 out[est] = dict(coverage=cov, mc_se=float(np.sqrt(cov * (1 - cov) / R)),62 mean_hw=float(np.mean(h)))63 out["cr0"]["rel_bias"] = float(np.mean(res["v0"]) / var_exact)64 out["cr1"]["rel_bias"] = float(np.mean(res["v1"]) / var_exact)65 out["iid"]["rel_bias"] = float(np.mean(res["s2"] / N) / var_exact)66 # MC SE of the relative-bias estimate (sd of V-hat / sqrt(R) / Var exact)67 out["cr0"]["rel_bias_se"] = float(np.std(res["v0"], ddof=1) / np.sqrt(R) / var_exact)68 out["cr1"]["rel_bias_se"] = float(np.std(res["v1"], ddof=1) / np.sqrt(R) / var_exact)69 # analytic E[V_CR0] (for validation): Var(T_g - n_g*ybar) summed over g, / N^270 VT = icc * n.astype(float) ** 2 + (1 - icc) * n71 ev0 = np.sum(VT * (1 - 2 * n / N) + n.astype(float) ** 2 * var_exact) / N ** 272 out["cr0"]["exact_rel_bias"] = float(ev0 / var_exact)73 out["cr1"]["exact_rel_bias"] = float(adj * ev0 / var_exact)74 out["meta"] = dict(var_exact=var_exact, emp_var_ybar=float(np.var(res["ybar"], ddof=1)),75 mean_ybar=float(np.mean(res["ybar"])), tcrit=float(tcrit),76 mean_q=float(np.mean(res["q"])), n_distinct_wcb_q=int(len(np.unique(np.round(res["q"], 10)))))77 return out787980results = {"task_id": "T2-simulation"}81Rs, Bs = design.R.unique(), design.B.unique()82assert len(Rs) == 1 and len(Bs) == 183R, B = int(Rs[0]), int(Bs[0])84cov, mcse, rb, rb_se, rb_exact = {}, {}, {}, {}, {}85base_cov, base_rb, cellmeta = {}, {}, {}86rows = []87for _, row in design.iterrows():88 o = run_cell(row)89 for est in ["cr0", "cr1", "wcb"]:90 k = key(est, row.G, row.icc)91 cov[k], mcse[k] = o[est]["coverage"], o[est]["mc_se"]92 if est != "wcb":93 rb[k], rb_se[k] = o[est]["rel_bias"], o[est]["rel_bias_se"]94 rb_exact[k] = o[est]["exact_rel_bias"]95 rows.append(dict(estimator=est, cell_id=row.cell_id, G=int(row.G), icc=float(row.icc),96 coverage=o[est]["coverage"], mc_se_coverage=o[est]["mc_se"],97 rel_bias_var=(o[est]["rel_bias"] if est != "wcb" else None),98 mean_ci_halfwidth=o[est]["mean_hw"], R=int(row.R)))99 kb = key("iid", row.G, row.icc)100 base_cov[kb], base_rb[kb] = o["iid"]["coverage"], o["iid"]["rel_bias"]101 cellmeta[row.cell_id] = o["meta"]102 print(f"{row.cell_id} G={row.G} icc={row.icc}: cr0={o['cr0']['coverage']:.4f} "103 f"cr1={o['cr1']['coverage']:.4f} wcb={o['wcb']['coverage']:.4f} iid={o['iid']['coverage']:.4f}")104105table = pd.DataFrame(rows)106table = table.sort_values(["estimator", "cell_id"]).reset_index(drop=True)107table.to_csv(os.path.join(OUT, "table1.csv"), index=False)108109failure = sorted(k for k, v in cov.items() if v < 0.90)110nominal = sorted(k for k, v in cov.items() if abs(v - 0.95) <= 2 * mcse[k])111avg_small = {est: float(np.mean([cov[key(est, 5, i)] for i in [0.0, 0.2, 0.5]])) for est in ["cr0", "cr1", "wcb"]}112best = max(avg_small, key=avg_small.get)113114data_issues = [115 {"issue": "ICC = 0.0 cells imply u_g ~ N(0, 0), a degenerate distribution.",116 "handling": "u_g drawn as sqrt(icc)*Z, which is identically 0 for icc = 0; DGP otherwise unchanged."},117 {"issue": "With G = 5 and Rademacher weights there are only 2^5 = 32 distinct weight vectors and |t*| is invariant to a global sign flip, so at most 16 distinct |t*| values exist; the B = 399 bootstrap quantile is highly discrete (draws repeat the same few values).",118 "handling": "Kept exactly as specified (no Webb weights or enumeration); flagged as a limitation in the report."},119 {"issue": "'0.95 empirical quantile' of |t*| is not uniquely defined for finite B.",120 "handling": "numpy.quantile default (linear interpolation, position 0.95*(B-1)); documented in report."},121]122results.update({123 "R": R, "B": B, "n_cells": int(len(design)),124 "coverage": cov, "mc_se_coverage": mcse, "rel_bias_var": rb,125 "failure_cells": failure, "n_failure_cells": len(failure), "nominal_cells": nominal,126 "best_procedure_small_G": best,127 "feasibility": "complete", "feasibility_reason": None,128 "data_issues": data_issues,129 "n_input_rows": {"design.csv": int(len(design))},130 "n_analysed": int(len(design) * R),131 "random_seed": int(design.loc[design.cell_id == "C1", "seed"].iloc[0]),132 # supplementary keys133 "avg_coverage_small_G": avg_small,134 "mc_se_rel_bias_var": rb_se,135 "exact_rel_bias_var_analytic": rb_exact,136 "baseline_iid_coverage": base_cov,137 "baseline_iid_rel_bias_var": base_rb,138 "cell_diagnostics": cellmeta,139 "seeds": {r.cell_id: int(r.seed) for _, r in design.iterrows()},140})141with open(os.path.join(OUT, "results.json"), "w") as f:142 json.dump(results, f, indent=2)143144# Figure145fig, axes = plt.subplots(1, 3, figsize=(13, 4.2), sharey=True)146Gs = sorted(design.G.unique())147for ax, icc in zip(axes, [0.0, 0.2, 0.5]):148 for est, mk in [("cr0", "o"), ("cr1", "s"), ("wcb", "^")]:149 p = np.array([cov[key(est, G, icc)] for G in Gs])150 s = np.array([mcse[key(est, G, icc)] for G in Gs])151 ax.errorbar(Gs, p, yerr=1.96 * s, marker=mk, capsize=3, label=est.upper())152 ax.axhline(0.95, color="k", ls="--", lw=1, label="nominal 0.95")153 ax.axhline(0.90, color="r", ls=":", lw=1, label="failure line 0.90")154 ax.set_xscale("log"); ax.set_xticks(Gs); ax.set_xticklabels([str(g) for g in Gs]); ax.minorticks_off()155 ax.set_title(f"ICC = {icc:.1f}"); ax.set_xlabel("number of clusters G")156axes[0].set_ylabel("coverage of nominal 95% CI (±1.96 MC SE)")157axes[0].legend(fontsize=8, loc="lower right")158fig.suptitle(f"Coverage of cluster-robust CIs for a mean (R={R}, B={B})")159fig.tight_layout()160fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130)161print("failure cells:", failure)162print("best procedure (G=5):", best, avg_small)163