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.5Deterministic: each cell uses numpy.random.default_rng(seed of that cell). Within a replication6the draw order is: u (G normals), e (N normals), then a B x G Rademacher weight matrix.7"""8import json9import os1011import numpy as np12import pandas as pd13from scipy import stats14import matplotlib15matplotlib.use("Agg")16import matplotlib.pyplot as plt17import matplotlib.ticker1819ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))20DATA = os.path.join(ROOT, "data", "design.csv")21OUT = os.path.join(ROOT, "out")22EST = ["cr0", "cr1", "wcb"]232425def key(est, G, icc):26 return f"{est}_G{int(G)}_icc{float(icc):.1f}"272829def simulate_cell(row):30 G, icc, N, R, B = int(row.G), float(row.icc), int(row.N), int(row.R), int(row.B)31 mu, tv, level = float(row.mu), float(row.total_variance), float(row.nominal_level)32 sizes = np.array([int(s) for s in str(row.cluster_sizes).split("|")])33 assert len(sizes) == G and sizes.sum() == N34 cid = np.repeat(np.arange(G), sizes)35 var_u, var_e = icc * tv, (1.0 - icc) * tv36 var_ybar = (var_u * np.sum(sizes ** 2) + var_e * N) / N ** 2 # exact Var(ybar)37 z = stats.norm.ppf(1 - (1 - level) / 2)38 z = 1.96 if abs(z - 1.96) < 1e-3 else z # TASK specifies 1.9639 tcrit = stats.t.ppf(1 - (1 - level) / 2, G - 1)40 cr1f = G / (G - 1)41 rng = np.random.default_rng(int(row.seed))4243 out = {k: np.empty(R) for k in ["ybar", "v0", "v1", "q", "v_iid"]}44 n_degenerate = 045 for r in range(R):46 u = rng.normal(0.0, np.sqrt(var_u), G)47 e = rng.normal(0.0, np.sqrt(var_e), N)48 y = mu + u[cid] + e49 ybar = y.mean()50 res = y - ybar51 S = np.bincount(cid, weights=res, minlength=G) # cluster sums of residuals52 v0 = np.sum(S ** 2) / N ** 253 v1 = cr1f * v054 # wild cluster bootstrap-t (unrestricted, Rademacher)55 W = rng.choice(np.array([-1.0, 1.0]), size=(B, G))56 delta = (W @ S) / N # ybar* - ybar57 Sstar = W * S[None, :] - delta[:, None] * sizes[None, :] # cluster sums of bootstrap residuals58 v1s = cr1f * np.sum(Sstar ** 2, axis=1) / N ** 259 ok = v1s > 060 n_degenerate += int((~ok).sum())61 tstar = np.abs(delta[ok]) / np.sqrt(v1s[ok])62 q = np.quantile(tstar, level, method="inverted_cdf") # = ceil(0.95*B)-th order stat (380th of 399)63 out["ybar"][r], out["v0"][r], out["v1"][r], out["q"][r] = ybar, v0, v1, q64 out["v_iid"][r] = np.sum(res ** 2) / (N - 1) / N # naive iid variance (baseline only)6566 err = np.abs(out["ybar"] - mu)67 hw = {68 "cr0": z * np.sqrt(out["v0"]),69 "cr1": tcrit * np.sqrt(out["v1"]),70 "wcb": out["q"] * np.sqrt(out["v1"]),71 "iid": z * np.sqrt(out["v_iid"]),72 }73 res = {}74 for est, h in hw.items():75 cov = float(np.mean(err <= h))76 res[est] = dict(coverage=cov, mc_se=float(np.sqrt(cov * (1 - cov) / R)),77 mean_hw=float(np.mean(h)))78 res["cr0"]["rel_bias"] = float(np.mean(out["v0"]) / var_ybar)79 res["cr1"]["rel_bias"] = float(np.mean(out["v1"]) / var_ybar)80 res["cr0"]["rel_bias_mcse"] = float(np.std(out["v0"], ddof=1) / np.sqrt(R) / var_ybar)81 res["cr1"]["rel_bias_mcse"] = float(np.std(out["v1"], ddof=1) / np.sqrt(R) / var_ybar)82 res["iid"]["rel_bias"] = float(np.mean(out["v_iid"]) / var_ybar)83 # exact E[V_CR0]/Var(ybar) (analytic check of the Monte Carlo bias estimate)84 Mm = np.eye(N) - 1.0 / N85 Sig = var_u * (cid[:, None] == cid[None, :]) + var_e * np.eye(N)86 A = np.zeros((G, N)); A[cid, np.arange(N)] = 1.087 ev0 = float(np.trace(A @ Mm @ Sig @ Mm @ A.T) / N ** 2)88 res["cr0"]["rel_bias_exact"] = ev0 / var_ybar89 res["cr1"]["rel_bias_exact"] = cr1f * ev0 / var_ybar90 res["meta"] = dict(var_ybar=float(var_ybar), n_degenerate_boot=n_degenerate,91 mean_ybar=float(out["ybar"].mean()),92 var_ybar_mc=float(np.var(out["ybar"], ddof=1)),93 n_distinct_q_min=int(len(np.unique(np.round(out["q"], 12)))))94 return res959697def main():98 design = pd.read_csv(DATA)99 R_all, B_all = design.R.unique(), design.B.unique()100 assert len(R_all) == 1 and len(B_all) == 1101 R, B = int(R_all[0]), int(B_all[0])102103 rows, cov, mcse, rb, extra = [], {}, {}, {}, {}104 base_cov, base_mcse = {}, {}105 n_degenerate = 0106 for _, row in design.iterrows():107 r = simulate_cell(row)108 n_degenerate += r["meta"]["n_degenerate_boot"]109 for est in EST:110 k = key(est, row.G, row.icc)111 cov[k], mcse[k] = r[est]["coverage"], r[est]["mc_se"]112 if est != "wcb":113 rb[k] = r[est]["rel_bias"]114 extra[k + "_mcse"] = r[est]["rel_bias_mcse"]115 extra[k + "_exact"] = r[est]["rel_bias_exact"]116 rows.append(dict(estimator=est, cell_id=row.cell_id, G=int(row.G), icc=float(row.icc),117 coverage=r[est]["coverage"], mc_se_coverage=r[est]["mc_se"],118 rel_bias_var=r[est].get("rel_bias", np.nan) if est != "wcb" else np.nan,119 mean_ci_halfwidth=r[est]["mean_hw"], R=int(row.R)))120 kb = f"iid_G{int(row.G)}_icc{float(row.icc):.1f}"121 base_cov[kb], base_mcse[kb] = r["iid"]["coverage"], r["iid"]["mc_se"]122 extra[f"meta_{row.cell_id}"] = r["meta"]123124 table = pd.DataFrame(rows)125 table.to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.6f")126127 failure = [k for k, v in cov.items() if v < 0.90]128 nominal = [k for k, v in cov.items() if abs(v - 0.95) <= 2 * mcse[k]]129 g5 = {est: float(np.mean([cov[key(est, 5, i)] for i in sorted(design.icc.unique())])) for est in EST}130 best = max(g5, key=lambda e: g5[e])131132 data_issues = [{133 "issue": "With G=5 and Rademacher weights only 2^5=32 distinct bootstrap weight vectors exist "134 "(16 distinct |t*| values), so the B=399 draws resample a tiny support and the WCB "135 "critical value is coarse and typically large.",136 "handling": "Procedure kept exactly as specified (no Webb weights, no enumeration); flagged as a "137 "limitation in the report."}]138 if n_degenerate:139 data_issues.append({"issue": f"{n_degenerate} bootstrap draws had V*_CR1 = 0",140 "handling": "dropped from that replication's |t*| set before taking the quantile"})141142 results = {143 "task_id": "T2-simulation",144 "R": R, "B": B, "n_cells": int(len(design)),145 "coverage": cov, "mc_se_coverage": mcse, "rel_bias_var": rb,146 "failure_cells": failure, "n_failure_cells": len(failure),147 "nominal_cells": nominal,148 "best_procedure_small_G": best,149 "feasibility": "complete", "feasibility_reason": None,150 "data_issues": data_issues,151 "n_input_rows": {"design.csv": int(len(design))},152 "n_analysed": int(len(design) * R),153 "random_seed": int(design.loc[design.cell_id == "C1", "seed"].iloc[0]),154 "mean_coverage_G5": g5,155 "rel_bias_var_mcse": {k.replace("_mcse", ""): v for k, v in extra.items() if k.endswith("_mcse")},156 "rel_bias_var_exact": {k.replace("_exact", ""): v for k, v in extra.items() if k.endswith("_exact")},157 "baseline_iid_coverage": base_cov,158 "baseline_iid_mc_se_coverage": base_mcse,159 "cell_diagnostics": {k[5:]: v for k, v in extra.items() if k.startswith("meta_")},160 "seeds": {str(c): int(s) for c, s in zip(design.cell_id, design.seed)},161 "wcb_quantile_rule": "inverted_cdf 0.95 quantile of |t*| (380th order statistic of 399)",162 }163 with open(os.path.join(OUT, "results.json"), "w") as f:164 json.dump(results, f, indent=2)165166 # figure167 iccs = sorted(design.icc.unique())168 Gs = sorted(design.G.unique())169 fig, axes = plt.subplots(1, len(iccs), figsize=(13, 4.2), sharey=True)170 styles = {"cr0": ("o-", "CR0 (z)"), "cr1": ("s-", "CR1 (t, G-1)"), "wcb": ("^-", "WCB-t (Rademacher)")}171 for ax, icc in zip(axes, iccs):172 for est in EST:173 p = np.array([cov[key(est, g, icc)] for g in Gs])174 se = np.array([mcse[key(est, g, icc)] for g in Gs])175 ax.errorbar(Gs, p, yerr=1.96 * se, fmt=styles[est][0], label=styles[est][1], capsize=3)176 ax.axhline(0.95, color="k", ls="--", lw=1, label="nominal 0.95")177 ax.axhline(0.90, color="r", ls=":", lw=1, label="failure line 0.90")178 ax.set_xscale("log"); ax.set_xticks(Gs); ax.set_xticklabels([str(g) for g in Gs])179 ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())180 ax.set_title(f"ICC = {icc:.1f}"); ax.set_xlabel("number of clusters G")181 axes[0].set_ylabel("coverage (±1.96 MC SE)")182 axes[0].legend(fontsize=8, loc="lower right")183 fig.suptitle(f"Coverage of nominal 95% CIs for the mean (R={R}, B={B})")184 fig.tight_layout()185 fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)186 plt.close(fig)187188 print(table.to_string(index=False))189 print("failure cells:", failure)190 print("best procedure at G=5:", best, g5)191192193if __name__ == "__main__":194 main()195