# Finite-sample coverage of cluster-robust 95% intervals for a mean (T2-simulation)

All numbers below come from out/results.json (key in brackets) or the identical values in out/table1.csv, produced by out/analysis.py (run: ~/wp29/venv/bin/python out/analysis.py, ≈2 s, deterministic).

# Question

With observations grouped in few, unequal-sized clusters, how close to 95% is the actual coverage of three cluster-robust intervals for the population mean μ: CR0 with a normal critical value, CR1 with a t(G−1) critical value, and the wild cluster bootstrap-t (WCB, Rademacher, unrestricted)? For each of 9 cells (G ∈ {5,10,30} × ICC ∈ {0, 0.2, 0.5}) we estimate coverage, its Monte Carlo (MC) SE, the relative bias E[V̂]/Var(ȳ) of the CR0/CR1 variance estimators, and the failure regimes (coverage < 0.90).

# Data and exclusions

No observed data. The only input is data/design.csv (9 rows [n_input_rows], SHA-256 f4ff004e…, audited in out/state/data_audit.md). Cluster sizes repeat 8,16,32,64 (N = 128, 264, 864); all sizes sum to N. Nothing was excluded (out/exclusions.csv: 9 → 9 cells). R = 2000 replications [R], B = 399 bootstrap draws [B], 9 cells [n_cells], 18 000 simulated samples in total [n_analysed]. Seeds are the per-cell seeds in the design (C1: 20261101 [random_seed]). Design notes recorded in [data_issues]: (i) at ICC = 0 the cluster effect has variance 0 (drawn as exact zeros); (ii) at G = 5 Rademacher weights take only 2⁵ = 32 distinct values (16 distinct |t*|), so the WCB critical value is coarse — kept as specified.

# Method

DGP per replication: y_ig = μ + u_g + e_ig, u_g ~ N(0, ICC), e_ig ~ N(0, 1−ICC), μ = 0. Per cell: numpy.random.default_rng(seed); draw u (R×G), then e (R×N), then Rademacher weights (B×G per replication, in chunks of 200 replications). With ê = y − ȳ and S_g = Σ_i ê_ig:

• CR0: V̂ = Σ S_g²/N², interval ȳ ± 1.96√V̂ (1.959964 used).
• CR1: V̂ = G/(G−1)·V̂_CR0, interval ȳ ± t_{0.975,G−1}√V̂.
• WCB: y* = ȳ + w_g ê_ig; t*_b = (ȳ* − ȳ)/√V̂*_CR1 computed on the bootstrap sample (residuals re-centred at ȳ*); q = 0.95 empirical quantile of |t*|, implemented as the ⌈0.95·B⌉ = 380th order statistic of the 399 sorted |t*| (the inverse-ECDF quantile; [wcb_quantile_rule]); interval ȳ ± q√V̂_CR1. Coverage = share of R intervals with |ȳ − μ| ≤ half-width; MC SE = √(p̂(1−p̂)/R). Relative bias = mean(V̂)/Var(ȳ), with exact Var(ȳ) = [ICC Σn_g² + (1−ICC)N]/N². Baseline (P2): naive i.i.d. interval ȳ ± 1.96·s/√N on the same draws [baseline_naive_iid_coverage].

# Results

Coverage (MC SE in parentheses) [coverage, mc_se_coverage]; relative variance bias [rel_bias_var]:

Figure: out/figure1.png (coverage vs G, one panel per ICC, ±1.96 MC SE bars, lines at 0.95 and 0.90).

Key findings:

1. Failure regimes [failure_cells, n_failure_cells = 11]: CR0 at G = 5 and G = 10 for every ICC (6 pairs); CR1 at G = 5 with ICC = 0.2 and 0.5; WCB at G = 5 for every ICC. No procedure fails at G = 30.
2. Best at G = 5 [best_procedure_small_G = cr1; avg_coverage_small_G: cr1 0.893, wcb 0.857, cr0 0.746]. Even CR1 under-covers at G = 5 when ICC > 0 (0.882, 0.880).
3. Sandwich variance is biased downward in every cell [rel_bias_var]: CR0 captures only about half of Var(ȳ) at G = 5 with ICC > 0 (0.507, 0.505) and 0.92–0.93 at G = 30. The CR1 G/(G−1) factor reduces but does not remove this (0.63 at G = 5, ICC > 0). The bias is worse when ICC > 0 because the unequal cluster sizes give the largest cluster (64 obs.) high leverage. The MC estimates agree with the exact analytic E[V̂_CR0]/Var(ȳ) in all cells (|z| ≤ 2.0; [diagnostics.*.exact_rel_bias_cr0], out/state/checks.md).
4. Consistent with nominal 95% [nominal_cells, 6 pairs]: CR1 at (G=10, ICC=0), (G=30, ICC=0), (G=30, ICC=0.2); WCB at G = 30 for all three ICCs. CR0 is never consistent with nominal.
5. WCB vs CR1: at G = 10 with ICC > 0, WCB (0.939, 0.937) beats CR1 (0.917, 0.917); at G = 5, WCB is worse than CR1 (e.g., 0.869 vs 0.919 at ICC = 0) even though its intervals are wider on average [table1.csv mean_ci_halfwidth: e.g. 0.244 vs 0.205 at G=5, ICC=0]. So at G = 5, WCB does not lose coverage because its intervals are too narrow on average. It loses because its interval width is poorly matched to the actual error from one sample to the next, given the small discrete set of 2⁵ weight vectors.
6. Baseline [baseline_naive_iid_coverage]: ignoring clustering gives 0.944–0.956 at ICC = 0 but only 0.30–0.48 at ICC > 0, so all three cluster-robust procedures are large improvements over the naive interval when clustering is present.

# Uncertainty

• Coverage MC SE ranges from 0.0048 to 0.0100 [mc_se_coverage]; a 95% MC interval for a coverage is about ±0.010–0.020. Classifications near a threshold are therefore uncertain. Close to 0.90: CR1 at G = 5 with ICC = 0.2/0.5 (0.882, 0.880; ≈2.5–2.8 SE below 0.90), CR0 at (G = 10, ICC = 0) (0.884, ≈2.3 SE below) and WCB at (G = 5, ICC = 0) (0.869, ≈4 SE below). Close to the nominal band: WCB at (G = 30, ICC = 0.5) (0.942, ~1.6 SE) and CR1 at (G = 30, ICC = 0.5) (0.929, ~3.7 SE, not nominal). A different seed could change the classification of pairs within about 2 SE of a threshold.
• The MC SE of the relative-bias estimates is 0.007–0.016 [rel_bias_var_mc_se] (SD of V̂ / √R / Var(ȳ)), so bias statements to two decimals are reliable. The exact analytic values agree with them (check in out/state/checks.md).
• The empirical variance of ȳ over replications is within 0.97–1.05 of the exact Var(ȳ) [diagnostics.*.emp_var_ybar / var_exact], which matches the R = 2000 sampling error.
• B = 399 adds bootstrap noise to each WCB critical value. That noise is already part of the WCB coverage estimate and its binomial MC SE. It is not separated out.

# Limitations and open issues

• Distributions: only Gaussian u and e. Skewed or heavy-tailed cluster effects or errors, which usually hurt few-cluster inference more, are not studied.
• Heterogeneity: one fixed unbalanced size pattern (8,16,32,64), homoskedastic errors, a common ICC across clusters, and no size-dependent cluster effects. Results with balanced clusters or different leverage patterns may differ.
• Regressors: the estimand is a mean (intercept-only model). Regression slopes, fixed effects and treatment designs with few treated clusters are not covered.
• Other weights/variants: only Rademacher weights with the unrestricted (non-null-imposed) WCB. Webb 6-point weights (relevant at G = 5, where Rademacher gives only 32 combinations), restricted WCB, CR2/CR3 (Bell–McCaffrey) corrections and degrees-of-freedom adjustments were not studied.
• The grid has only G ∈ {5,10,30} and R = 2000. Where coverage crosses 0.90 between these G values is not determined.
• The WCB quantile convention (380th order statistic of 399) is our reading of "0.95 empirical quantile". Other interpolation rules would shift q slightly.