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.4Randomness: one numpy default_rng(seed) per cell (seed from design.csv). Per replication, draws in5fixed order: u (G normals), e (N normals), W (B x G Rademacher via rng.integers(0,2)).6"""7import json, os8import numpy as np9import pandas as pd10from scipy import stats11import matplotlib12matplotlib.use("Agg")13import matplotlib.pyplot as plt1415ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))16DATA = os.path.join(ROOT, "data", "design.csv")17OUT = os.path.join(ROOT, "out")18EST = ["cr0", "cr1", "wcb"]192021def simulate_cell(row):22 G, icc, N, R, B, mu = int(row.G), float(row.icc), int(row.N), int(row.R), int(row.B), float(row.mu)23 sizes = np.array([int(x) for x in row.cluster_sizes.split("|")])24 assert len(sizes) == G and sizes.sum() == N25 gid = np.repeat(np.arange(G), sizes)26 var_exact = (icc * np.sum(sizes.astype(float) ** 2) + (1 - icc) * N) / N ** 227 tcrit = stats.t.ppf(0.975, G - 1)28 zcrit = 1.9629 adj = G / (G - 1)30 rng = np.random.default_rng(int(row.seed))31 res = {k: np.empty(R) for k in ["ybar", "v0", "v1", "q", "v_iid"]}32 n_degenerate = 033 for r in range(R):34 u = rng.normal(0.0, np.sqrt(icc), G) # scale 0 when icc=0 -> u = 035 e = rng.normal(0.0, np.sqrt(1 - icc), N)36 y = mu + u[gid] + e37 ybar = y.mean()38 ehat = y - ybar39 S = np.bincount(gid, weights=ehat, minlength=G) # cluster residual sums40 v0 = np.sum(S ** 2) / N ** 241 v1 = adj * v042 # wild cluster bootstrap-t (unrestricted, Rademacher), vectorised over B43 W = rng.integers(0, 2, size=(B, G)) * 2 - 144 d = (W * S).sum(axis=1) / N # ybar* - ybar45 Sstar = W * S - np.outer(d, sizes) # cluster sums of bootstrap residuals46 v1s = adj * np.sum(Sstar ** 2, axis=1) / N ** 247 if np.any(v1s <= 0):48 n_degenerate += int(np.sum(v1s <= 0))49 v1s = np.where(v1s <= 0, np.nan, v1s)50 tabs = np.abs(d) / np.sqrt(v1s)51 q = np.nanquantile(tabs, 0.95) # numpy default (linear interpolation)52 res["ybar"][r], res["v0"][r], res["v1"][r], res["q"][r] = ybar, v0, v1, q53 res["v_iid"][r] = np.var(y, ddof=1) / N # naive i.i.d. baseline54 hw = {"cr0": zcrit * np.sqrt(res["v0"]), "cr1": tcrit * np.sqrt(res["v1"]),55 "wcb": res["q"] * np.sqrt(res["v1"]), "iid": zcrit * np.sqrt(res["v_iid"])}56 out = {"var_exact": var_exact, "n_degenerate": n_degenerate, "R": R}57 for k, h in hw.items():58 cov = np.abs(res["ybar"] - mu) <= h59 out[k] = {"coverage": float(cov.mean()), "hw": float(h.mean()), "ind": cov.astype(float)}60 for k, v in [("cr0", res["v0"]), ("cr1", res["v1"])]:61 out[k]["rel_bias"] = float(v.mean() / var_exact)62 out[k]["rel_bias_mcse"] = float(v.std(ddof=1) / np.sqrt(R) / var_exact)63 out["mc_var_ybar_ratio"] = float(np.var(res["ybar"], ddof=1) / var_exact)64 out["mean_ybar"] = float(res["ybar"].mean())65 return out666768def main():69 d = pd.read_csv(DATA)70 R_all, B_all = d.R.unique(), d.B.unique()71 assert len(R_all) == 1 and len(B_all) == 172 R, B = int(R_all[0]), int(B_all[0])73 cov, se, rb, rb_se, base, chk = {}, {}, {}, {}, {}, {}74 rows, issues = [], []75 n_deg = 076 ind_g5 = {k: [] for k in EST}77 for row in d.itertuples(index=False):78 o = simulate_cell(row)79 n_deg += o["n_degenerate"]80 tag = f"G{int(row.G)}_icc{float(row.icc):.1f}"81 for k in EST:82 key = f"{k}_{tag}"83 p = o[k]["coverage"]84 cov[key] = p85 se[key] = float(np.sqrt(p * (1 - p) / o["R"]))86 if k != "wcb":87 rb[key] = o[k]["rel_bias"]88 rb_se[key] = o[k]["rel_bias_mcse"]89 rows.append({"estimator": k, "cell_id": row.cell_id, "G": int(row.G), "icc": float(row.icc),90 "coverage": p, "mc_se_coverage": se[key],91 "rel_bias_var": o[k].get("rel_bias", None), "mean_ci_halfwidth": o[k]["hw"], "R": o["R"]})92 if int(row.G) == 5:93 for k in EST:94 ind_g5[k].append(o[k]["ind"])95 base[f"iid_{tag}"] = o["iid"]["coverage"]96 chk[tag] = {"mc_var_ybar_over_exact": o["mc_var_ybar_ratio"], "mean_ybar": o["mean_ybar"]}97 print(f"{row.cell_id} {tag}: " + ", ".join(f"{k}={cov[f'{k}_{tag}']:.4f}" for k in EST))98 table = pd.DataFrame(rows)99 table.to_csv(os.path.join(OUT, "table1.csv"), index=False)100101 failure = [k for k, p in cov.items() if p < 0.90]102 nominal = [k for k, p in cov.items() if abs(p - 0.95) <= 2 * se[k]]103 g5 = {k: np.mean([cov[f"{k}_G5_icc{i:.1f}"] for i in (0.0, 0.2, 0.5)]) for k in EST}104 best = max(g5, key=g5.get)105 # paired MC SE of the difference in average G=5 coverage between best and each other procedure106 # (same replications -> paired; cells independent -> variance of mean over 3 cells = sum/9)107 g5_diff = {}108 for k in EST:109 if k == best:110 continue111 cell_vars = [np.var(a - b, ddof=1) / len(a) for a, b in zip(ind_g5[best], ind_g5[k])]112 g5_diff[f"{best}_minus_{k}"] = {"diff": float(g5[best] - g5[k]), "mc_se": float(np.sqrt(sum(cell_vars)) / len(cell_vars))}113 near = [k for k, p in cov.items() if abs(p - 0.90) <= 2 * se[k]]114 if 5 in d.G.values:115 issues.append({"issue": "With G=5 and Rademacher weights only 2^5=32 distinct weight vectors exist (16 distinct |t*| values, as |t*| is invariant to w -> -w); B=399 draws therefore repeat these few values and the 0.95 quantile of |t*| is coarse.",116 "handling": "Kept as specified (B=399, Rademacher, random draws, numpy linear-interpolation quantile); discussed in report as a known feature of WCB with very few clusters."})117 if n_deg:118 issues.append({"issue": f"{n_deg} bootstrap draws had non-positive V_CR1*", "handling": "Excluded those draws from the |t*| quantile (nan-quantile)."})119 results = {120 "task_id": "T2-simulation", "R": R, "B": B, "n_cells": int(len(d)),121 "coverage": cov, "mc_se_coverage": se, "rel_bias_var": rb,122 "failure_cells": failure, "n_failure_cells": len(failure), "nominal_cells": nominal,123 "best_procedure_small_G": best,124 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,125 "n_input_rows": {"design.csv": int(len(d))},126 "n_analysed": int(len(d) * R), "random_seed": int(d.loc[d.cell_id == "C1", "seed"].iloc[0]),127 # supplementary keys128 "avg_coverage_G5": {k: float(v) for k, v in g5.items()},129 "mc_se_rel_bias_var": rb_se,130 "avg_coverage_G5_paired_diff": g5_diff,131 "cells_within_2se_of_failure_line": near,132 "baseline_iid_coverage": base,133 "check_mc_var_ybar_over_exact": chk,134 "n_degenerate_bootstrap_draws": int(n_deg),135 "wcb_quantile_method": "numpy.quantile linear interpolation (type 7) of |t*_b|, b=1..B",136 "seeds_by_cell": {r.cell_id: int(r.seed) for r in d.itertuples(index=False)},137 }138 with open(os.path.join(OUT, "results.json"), "w") as f:139 json.dump(results, f, indent=2)140141 # figure142 fig, axes = plt.subplots(1, 3, figsize=(13, 4.2), sharey=True)143 Gs = sorted(d.G.unique())144 style = {"cr0": ("o-", "CR0 (z)"), "cr1": ("s-", "CR1 (t, G-1)"), "wcb": ("^-", "WCB-t Rademacher")}145 for ax, icc in zip(axes, sorted(d.icc.unique())):146 for k in EST:147 p = np.array([cov[f"{k}_G{g}_icc{icc:.1f}"] for g in Gs])148 s = np.array([se[f"{k}_G{g}_icc{icc:.1f}"] for g in Gs])149 ax.errorbar(Gs, p, yerr=1.96 * s, fmt=style[k][0], capsize=3, label=style[k][1])150 ax.axhline(0.95, color="k", ls="--", lw=1, label="nominal 0.95")151 ax.axhline(0.90, color="r", ls=":", lw=1, label="failure line 0.90")152 ax.set_xscale("log"); ax.minorticks_off(); ax.set_xticks(Gs); ax.set_xticklabels(Gs)153 ax.set_title(f"ICC = {icc:.1f}"); ax.set_xlabel("number of clusters G")154 axes[0].set_ylabel("coverage (R=2000; bars ±1.96 MC SE)")155 axes[0].legend(loc="lower right", fontsize=8)156 fig.tight_layout()157 fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130)158 print("failure cells:", failure)159 print("best procedure at G=5:", best)160161162if __name__ == "__main__":163 main()164