1"""T2-simulation: finite-sample coverage of CR0, CR1 and wild-cluster-bootstrap-t CIs for a mean.23Run from the workspace root: ~/wp29/venv/bin/python out/analysis.py4Reads only data/design.csv; writes out/results.json, out/table1.csv, out/figure1.png.56RNG: for each cell, rng = numpy.random.default_rng(cell seed). Per replication the draw order is7 u = rng.standard_normal(G); e = rng.standard_normal(N); W = rng.integers(0, 2, (B, G))*2-1.8Baseline (P2): naive i.i.d. interval ybar +- 1.96*s/sqrt(N), ignoring clustering (baseline_* keys).9"""10import json11import os12import numpy as np13import pandas as pd14from scipy import stats15import matplotlib16matplotlib.use("Agg")17import matplotlib.pyplot as plt1819ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))20DATA = os.path.join(ROOT, "data", "design.csv")21OUT = os.path.join(ROOT, "out")2223design = pd.read_csv(DATA)24data_issues = []252627def key(est, G, icc):28 return f"{est}_G{int(G)}_icc{icc:.1f}"293031def simulate_cell(row):32 G, icc, N, R, B = int(row.G), float(row.icc), int(row.N), int(row.R), int(row.B)33 mu = float(row.mu)34 n = np.array([int(x) for x in str(row.cluster_sizes).split("|")])35 assert len(n) == G and n.sum() == N, f"inconsistent design row {row.cell_id}"36 gid = np.repeat(np.arange(G), n)37 var_exact = (icc * np.sum(n ** 2) + (1 - icc) * N) / N ** 238 z = 1.9639 tcrit = stats.t.ppf(0.975, G - 1)40 cfac = G / (G - 1)41 qidx = int(np.ceil((B + 1) * 0.95)) - 1 # 0-based index of the 380th order statistic for B=39942 rng = np.random.default_rng(int(row.seed))43 rec = {k: np.empty(R) for k in ["ybar", "v0", "v1", "q", "s2", "qalt"]}44 n_zero_vstar = 045 for r in range(R):46 u = rng.standard_normal(G) * np.sqrt(icc)47 e = rng.standard_normal(N) * np.sqrt(1 - icc)48 W = rng.integers(0, 2, size=(B, G)) * 2 - 149 y = mu + u[gid] + e50 ybar = y.mean()51 ehat = y - ybar52 S = np.bincount(gid, weights=ehat, minlength=G) # cluster residual sums53 v0 = np.sum(S ** 2) / N ** 254 v1 = cfac * v055 # wild cluster bootstrap: y* = ybar + w_g ehat_ig56 d = (W @ S) / N # ybar* - ybar, shape (B,)57 Sstar = W * S[None, :] - d[:, None] * n[None, :] # cluster sums of bootstrap residuals58 vstar = cfac * np.sum(Sstar ** 2, axis=1) / N ** 259 n_zero_vstar += int(np.sum(vstar <= 0))60 tabs = np.abs(d) / np.sqrt(vstar)61 q = np.sort(tabs)[qidx]62 rec["qalt"][r] = np.quantile(tabs, 0.95) # sensitivity: numpy linear-interpolation quantile63 rec["ybar"][r], rec["v0"][r], rec["v1"][r], rec["q"][r] = ybar, v0, v1, q64 rec["s2"][r] = ehat.var(ddof=1)65 ybar, v0, v1, q = rec["ybar"], rec["v0"], rec["v1"], rec["q"]66 hw = {"cr0": z * np.sqrt(v0), "cr1": tcrit * np.sqrt(v1), "wcb": q * np.sqrt(v1),67 "iid": z * np.sqrt(rec["s2"] / N), "wcb_alt": rec["qalt"] * np.sqrt(v1)}68 out = {}69 for est, h in hw.items():70 cov = float(np.mean(np.abs(ybar - mu) <= h))71 out[est] = {"coverage": cov, "mc_se": float(np.sqrt(cov * (1 - cov) / R)),72 "hw": float(np.mean(h))}73 out["cr0"]["rel_bias"] = float(np.mean(v0) / var_exact)74 out["cr1"]["rel_bias"] = float(np.mean(v1) / var_exact)75 out["cr0"]["rel_bias_se"] = float(np.std(v0, ddof=1) / np.sqrt(R) / var_exact)76 out["cr1"]["rel_bias_se"] = float(np.std(v1, ddof=1) / np.sqrt(R) / var_exact)77 out["meta"] = {"var_exact": float(var_exact), "mc_var_ybar": float(np.var(ybar, ddof=1)),78 "n_zero_vstar": n_zero_vstar, "n_distinct_q": int(len(np.unique(np.round(q, 12))))}79 return out808182ests = ["cr0", "cr1", "wcb"]83rows, coverage, mcse, relbias = [], {}, {}, {}84base_cov, base_se, meta, rb_se, wcb_alt = {}, {}, {}, {}, {}85for _, row in design.iterrows():86 res = simulate_cell(row)87 G, icc = int(row.G), float(row.icc)88 for est in ests:89 k = key(est, G, icc)90 coverage[k] = res[est]["coverage"]91 mcse[k] = res[est]["mc_se"]92 rb = res[est].get("rel_bias")93 if rb is not None:94 relbias[k] = rb95 rb_se[k] = res[est]["rel_bias_se"]96 rows.append({"estimator": est, "cell_id": row.cell_id, "G": G, "icc": icc,97 "coverage": res[est]["coverage"], "mc_se_coverage": res[est]["mc_se"],98 "rel_bias_var": rb, "mean_ci_halfwidth": res[est]["hw"], "R": int(row.R)})99 bk = f"G{G}_icc{icc:.1f}"100 base_cov[bk] = res["iid"]["coverage"]101 base_se[bk] = res["iid"]["mc_se"]102 wcb_alt[key("wcb", G, icc)] = res["wcb_alt"]["coverage"]103 meta[row.cell_id] = res["meta"]104 if res["meta"]["n_zero_vstar"] > 0:105 data_issues.append({"issue": f"{row.cell_id}: {res['meta']['n_zero_vstar']} bootstrap draws with zero V*",106 "handling": "t* set to inf/nan by numpy; kept as is"})107108table = pd.DataFrame(rows)109table["est_order"] = table.estimator.map({e: i for i, e in enumerate(ests)})110table = table.sort_values(["est_order", "cell_id"]).drop(columns="est_order")111table.to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.6f")112113failure = sorted([k for k, v in coverage.items() if v < 0.90])114nominal = sorted([k for k, v in coverage.items() if abs(v - 0.95) <= 2 * mcse[k]])115smallG = {e: float(np.mean([coverage[key(e, 5, i)] for i in [0.0, 0.2, 0.5]])) for e in ests}116best = max(smallG, key=smallG.get)117118if (design.G == 5).any():119 data_issues.append({"issue": "G=5 cells: Rademacher weights give only 2^5=32 distinct weight vectors (16 distinct |t*| values), "120 "so the WCB critical value takes few distinct values and B=399 draws mostly repeat them",121 "handling": "procedure kept exactly as specified (no Webb weights / enumeration); reported as a limitation"})122data_issues.append({"issue": "ICC=0.0 cells: cluster effect variance is zero (degenerate N(0,0))",123 "handling": "u_g drawn as sqrt(0)*N(0,1)=0 so the RNG stream structure is the same in every cell"})124125R = int(design.R.iloc[0]); B = int(design.B.iloc[0])126results = {127 "task_id": "T2-simulation",128 "R": R, "B": B, "n_cells": int(len(design)),129 "coverage": coverage, "mc_se_coverage": mcse, "rel_bias_var": relbias,130 "mc_se_rel_bias_var": rb_se,131 "failure_cells": failure, "n_failure_cells": len(failure), "nominal_cells": nominal,132 "best_procedure_small_G": best,133 "mean_coverage_G5": smallG,134 "feasibility": "complete", "feasibility_reason": None,135 "data_issues": data_issues,136 "n_input_rows": {"design.csv": int(len(design))},137 "n_analysed": int(design.R.sum()),138 "random_seed": int(design.loc[design.cell_id == "C1", "seed"].iloc[0]),139 "baseline_iid_coverage": base_cov, "baseline_iid_mc_se": base_se,140 "sensitivity_wcb_coverage_numpy_quantile": wcb_alt,141 "wcb_quantile_rule": "380th order statistic of B=399 |t*| values (ceil((B+1)*0.95))",142 "cell_diagnostics": meta,143}144with open(os.path.join(OUT, "results.json"), "w") as f:145 json.dump(results, f, indent=2)146147# Figure148fig, axes = plt.subplots(1, 3, figsize=(13, 4.2), sharey=True)149Gs = sorted(design.G.unique())150style = {"cr0": ("o-", "CR0 (z)"), "cr1": ("s-", "CR1 (t, G-1)"), "wcb": ("^-", "WCB-t (Rademacher)")}151for ax, icc in zip(axes, sorted(design.icc.unique())):152 for e in ests:153 c = np.array([coverage[key(e, g, icc)] for g in Gs])154 s = np.array([mcse[key(e, g, icc)] for g in Gs])155 ax.errorbar(Gs, c, yerr=1.96 * s, fmt=style[e][0], label=style[e][1], capsize=3)156 ax.axhline(0.95, color="k", ls="--", lw=1, label="nominal 0.95")157 ax.axhline(0.90, color="r", ls=":", lw=1, label="failure line 0.90")158 ax.set_xscale("log"); ax.minorticks_off(); ax.set_xticks(Gs); ax.set_xticklabels([str(g) for g in Gs])159 ax.set_title(f"ICC = {icc:.1f}"); ax.set_xlabel("number of clusters G")160axes[0].set_ylabel("coverage (±1.96 MC SE)")161axes[0].legend(fontsize=8, loc="lower right")162fig.suptitle(f"Coverage of nominal 95% cluster-robust CIs for the mean (R={R}, B={B})")163fig.tight_layout()164fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130)165166print(table.to_string(index=False))167print("failure cells:", failure)168print("baseline iid coverage:", base_cov)169print("WCB coverage, numpy-quantile sensitivity:", wcb_alt)170print("best procedure at G=5:", best, smallG)171