1"""T2-simulation: coverage of CR0, CR1 and wild-cluster-bootstrap-t intervals for a mean.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np, pandas as pd6from scipy import stats7import matplotlib8matplotlib.use("Agg")9import matplotlib.pyplot as plt1011OUT = "out"12design = pd.read_csv("data/design.csv")13issues = []1415def simulate_cell(row):16 sizes = np.array([int(s) for s in str(row.cluster_sizes).split("|")])17 G, icc, R, B, mu = int(row.G), float(row.icc), int(row.R), int(row.B), float(row.mu)18 N = int(sizes.sum())19 if N != int(row.N) or len(sizes) != G:20 issues.append({"issue": f"{row.cell_id}: cluster_sizes inconsistent with G/N",21 "handling": "used cluster_sizes as given"})22 rng = np.random.default_rng(int(row.seed))23 # data: all R replications at once (order: u then e, then bootstrap weights)24 u = rng.normal(0.0, np.sqrt(icc), size=(R, G)) if icc > 0 else np.zeros((R, G))25 e = rng.normal(0.0, np.sqrt(1.0 - icc), size=(R, N))26 gid = np.repeat(np.arange(G), sizes)27 y = mu + u[:, gid] + e # R x N28 ybar = y.mean(axis=1) # R29 resid = y - ybar[:, None]30 S = np.zeros((R, G))31 for g in range(G):32 S[:, g] = resid[:, gid == g].sum(axis=1) # cluster residual sums33 V0 = (S ** 2).sum(axis=1) / N ** 234 V1 = G / (G - 1) * V035 # wild cluster bootstrap-t (Rademacher, unrestricted); statistics depend only on cluster sums:36 # ybar* - ybar = sum_g w_g S_g / N ; bootstrap residual cluster sums = w_g S_g - n_g (ybar*-ybar)37 W = rng.choice(np.array([-1.0, 1.0]), size=(R, B, G))38 wS = W * S[:, None, :] # R x B x G39 d = wS.sum(axis=2) / N # R x B40 Sstar = wS - sizes[None, None, :] * d[:, :, None]41 V1star = G / (G - 1) * (Sstar ** 2).sum(axis=2) / N ** 242 with np.errstate(divide="ignore", invalid="ignore"):43 tstar = np.abs(d) / np.sqrt(V1star)44 n_degen = int(np.sum(V1star <= 0))45 tstar = np.where(V1star > 0, tstar, 0.0)46 # 0.95 empirical quantile of |t*|: order statistic ceil(0.95*(B+1)) = 380th of 39947 k = int(np.ceil(0.95 * (B + 1))) - 148 q = np.sort(tstar, axis=1)[:, k]49 z = stats.norm.ppf(0.975) # 1.95996...; task says 1.9650 z = 1.9651 tc = stats.t.ppf(0.975, G - 1)52 hw = {"cr0": z * np.sqrt(V0), "cr1": tc * np.sqrt(V1), "wcb": q * np.sqrt(V1)}53 var_exact = (icc * (sizes ** 2).sum() + (1 - icc) * N) / N ** 254 out = []55 for est, h in hw.items():56 cov = float(np.mean(np.abs(ybar - mu) <= h))57 se = float(np.sqrt(cov * (1 - cov) / R))58 rb = None59 if est == "cr0": rb = float(V0.mean() / var_exact)60 if est == "cr1": rb = float(V1.mean() / var_exact)61 out.append(dict(estimator=est, cell_id=row.cell_id, G=G, icc=icc, coverage=cov,62 mc_se_coverage=se, rel_bias_var=rb, mean_ci_halfwidth=float(h.mean()), R=R))63 extra = dict(cell_id=row.cell_id, n_degenerate_boot=n_degen,64 var_exact=var_exact, emp_var_ybar=float(ybar.var(ddof=1)),65 sd_V0=float(V0.std(ddof=1)), sd_V1=float(V1.std(ddof=1)),66 mean_q=float(q.mean()), t_crit=float(tc),67 mcse_relbias_cr0=float(V0.std(ddof=1) / np.sqrt(R) / var_exact),68 mcse_relbias_cr1=float(V1.std(ddof=1) / np.sqrt(R) / var_exact),69 emp_var_ratio=float(ybar.var(ddof=1) / var_exact))70 return out, extra7172rows, extras = [], []73for _, r in design.iterrows():74 o, x = simulate_cell(r)75 rows += o; extras.append(x)76 print(r.cell_id, [(d["estimator"], round(d["coverage"], 4)) for d in o])7778tab = pd.DataFrame(rows)79tab.to_csv(os.path.join(OUT, "table1.csv"), index=False)80pd.DataFrame(extras).to_csv(os.path.join(OUT, "cell_diagnostics.csv"), index=False)8182key = lambda d: f"{d['estimator']}_G{d['G']}_icc{d['icc']:.1f}"83coverage = {key(d): d["coverage"] for d in rows}84mcse = {key(d): d["mc_se_coverage"] for d in rows}85relb = {key(d): d["rel_bias_var"] for d in rows if d["estimator"] != "wcb"}86fail = [key(d) for d in rows if d["coverage"] < 0.90]87nominal = [key(d) for d in rows if abs(d["coverage"] - 0.95) <= 2 * d["mc_se_coverage"]]88small = tab[tab.G == 5].groupby("estimator").coverage.mean()89best = str(small.idxmax())9091for x in extras:92 if x["n_degenerate_boot"] > 0:93 issues.append({"issue": f"{x['cell_id']}: {x['n_degenerate_boot']} bootstrap draws with zero V*_CR1",94 "handling": "|t*| set to 0 for those draws"})95issues.append({"issue": "With G=5 Rademacher weights there are only 2^5=32 distinct weight vectors "96 "(16 distinct |t*| values up to sign), so the WCB critical value from B=399 draws is "97 "a quantile of a highly discrete distribution",98 "handling": "kept the specified procedure (B=399 random draws, no enumeration); noted in report"})99100res = {101 "task_id": "T2-simulation", "R": int(design.R.iloc[0]), "B": int(design.B.iloc[0]),102 "n_cells": int(len(design)), "coverage": coverage, "mc_se_coverage": mcse,103 "rel_bias_var": relb, "failure_cells": fail, "n_failure_cells": len(fail),104 "nominal_cells": nominal, "best_procedure_small_G": best,105 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,106 "n_input_rows": {"design.csv": int(len(design))},107 "n_analysed": int(design.R.sum()), "random_seed": int(design.seed.iloc[0]),108}109with open(os.path.join(OUT, "results.json"), "w") as f:110 json.dump(res, f, indent=2)111112fig, axes = plt.subplots(1, 3, figsize=(13, 4.2), sharey=True)113for ax, icc in zip(axes, sorted(tab.icc.unique())):114 for est, m in zip(["cr0", "cr1", "wcb"], ["o", "s", "^"]):115 s = tab[(tab.icc == icc) & (tab.estimator == est)].sort_values("G")116 ax.errorbar(s.G, s.coverage, yerr=1.96 * s.mc_se_coverage, marker=m, capsize=3, label=est.upper())117 ax.axhline(0.95, color="k", ls="--", lw=1, label="nominal 0.95")118 ax.axhline(0.90, color="r", ls=":", lw=1, label="failure 0.90")119 ax.set_xscale("log"); ax.set_xticks([5, 10, 30]); ax.set_xticklabels(["5", "10", "30"]); ax.minorticks_off()120 ax.set_title(f"ICC = {icc:.1f}"); ax.set_xlabel("number of clusters G")121axes[0].set_ylabel("coverage (±1.96 MC SE)"); axes[0].legend(loc="lower right", fontsize=8)122fig.suptitle("Coverage of nominal 95% cluster-robust intervals for the mean (R=2000 per cell)")123fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130)124print(json.dumps({k: res[k] for k in ["failure_cells", "nominal_cells", "best_procedure_small_G"]}, indent=1))125