# Finite-sample coverage of CR0, CR1 and wild-cluster-bootstrap-t intervals for a mean with few clusters

All numbers below are produced by out/analysis.py (run: ~/wp29/venv/bin/python out/analysis.py, ≈2 s) and are cited by their key in out/results.json. Full per-cell values are in out/table1.csv; the figure is out/figure1.png.

# Question

When N observations are grouped in few (G = 5, 10, 30) unbalanced clusters with intra-cluster correlation ICC ∈ {0, 0.2, 0.5}, what is the actual coverage of three nominal-95% cluster-robust intervals for the mean μ (CR0 with z critical value; CR1 with t_{G−1}; unrestricted Rademacher wild cluster bootstrap-t, WCB), how biased are the CR0/CR1 variance estimators relative to the exact Var(ȳ), and in which (procedure, cell) pairs does coverage fall below 0.90?

# Data and exclusions

No observed data. The only input is data/design.csv (n_input_rows: 9 rows; SHA-256 recorded in out/state/data_audit.md): 9 cells, cluster sizes 8,16,32,64 repeated (N = 128, 264, 864), R = 2000 (R), B = 399 (B), one seed per cell (C1 seed = random_seed; all seeds in seeds_by_cell). The audit found no missing values, duplicates or inconsistencies (sizes count = G, sum = N). No exclusions: all n_cells = 9 cells × 2000 replications = n_analysed = 18000 replications were simulated and analysed (out/exclusions.csv, one no-op step). feasibility = "complete". One design feature is recorded in data_issues: at G = 5 Rademacher weights take only 32 distinct values (16 distinct |t*|), so the bootstrap critical value is coarse; kept as specified. No bootstrap draw had a non-positive variance (n_degenerate_bootstrap_draws = 0).

# Method

Per replication (numpy default_rng(cell seed), draws in fixed order u, e, W): y_ig = μ + u_g + e_ig, u_g ~ N(0, ICC), e_ig ~ N(0, 1−ICC), μ = 0. With ê = y − ȳ and S_g = Σ_i ê_ig: V̂_CR0 = Σ S_g²/N², CI ȳ ± 1.96√V̂_CR0; V̂_CR1 = G/(G−1)·V̂_CR0, CI ȳ ± t_{0.975,G−1}√V̂_CR1; WCB: B×G Rademacher matrix W, ȳ*−ȳ = Σ_g w_g S_g/N, bootstrap cluster residual sums w_g S_g − n_g(ȳ*−ȳ), V̂*_CR1 from these, t*_b = (ȳ*−ȳ)/√V̂*_CR1, q = 0.95 quantile of |t*| (numpy linear interpolation, wcb_quantile_method), CI ȳ ± q√V̂_CR1. Coverage = share of R intervals containing μ; MC SE = √(p̂(1−p̂)/R); relative bias = mean(V̂)/Var(ȳ) with Var(ȳ) = [ICC Σn_g² + (1−ICC)N]/N². Failure regime: coverage < 0.90; nominal-consistent: |p̂ − 0.95| ≤ 2 MC SE. Baseline (P2): the naive i.i.d. interval ȳ ± 1.96·s/√N on the same draws (baseline_iid_coverage).

Verification (out/state/checks.md): byte-identical re-run; vectorised WCB matches a brute-force construction of y* (diff 7e-15); MC relative bias of CR0 matches the closed-form E[V̂_CR0]/Var(ȳ) within 2 MC SE in 8/9 cells (z = +2.46 in G5_icc0.2); MC Var(ȳ)/exact in [0.93, 1.06] (check_mc_var_ybar_over_exact).

# Results

Coverage (coverage, MC SE in mc_se_coverage):

• Failure regimes (failure_cells, n_failure_cells = 11): CR0 at every G = 5 and G = 10 cell (6 pairs); WCB at all three G = 5 cells; CR1 at G = 5 with ICC = 0.2 and 0.5. No procedure fails at G = 30.
• Nominal-consistent pairs (nominal_cells, 5 pairs): cr1_G10_icc0.0, cr1_G30_icc0.0, wcb_G30_icc0.0, wcb_G30_icc0.2, wcb_G30_icc0.5. CR0 is never consistent with 0.95, even at G = 30 (0.917–0.928).
• Best at G = 5 (best_procedure_small_G = "cr1"): average G = 5 coverage 0.891 for CR1, 0.853 for WCB, 0.740 for CR0 (avg_coverage_G5). The CR1 − WCB difference is 0.038 with paired MC SE 0.0035 (avg_coverage_G5_paired_diff), so the ranking is not a Monte Carlo artefact. At G = 10 with ICC > 0, WCB covers better than CR1 (0.930/0.933 vs 0.916/0.907).
• Variance bias (rel_bias_var): both sandwich estimators are biased downward, most at small G and high ICC. CR0: 0.66, 0.53, 0.48 at G = 5 (ICC 0, 0.2, 0.5); 0.86, 0.76, 0.77 at G = 10; 0.94, 0.93, 0.93 at G = 30. CR1: 0.83, 0.66, 0.60 at G = 5; 0.95, 0.85, 0.86 at G = 10; 0.98, 0.96, 0.96 at G = 30. The unbalanced sizes (one 64-observation cluster dominates) make the bias much worse than the balanced-cluster factor (G−1)/G.
• Baseline: the naive i.i.d. interval is fine at ICC = 0 but covers only 0.29–0.49 when ICC > 0 (baseline_iid_coverage), whatever G is. All three cluster-robust procedures improve on it there.
• Mean CI half-widths are in table1.csv (mean_ci_halfwidth). At G = 5 the WCB interval is on average wider than CR1 but still covers less. One plausible explanation (not tested here) is that q varies a lot across replications, because the bootstrap distribution has only 16 distinct |t*| values.

# Uncertainty

Each coverage is a binomial proportion over R = 2000 independent replications. Its MC SE (mc_se_coverage) is 0.005–0.010, so a 95% MC interval is about ±0.010–0.020. The figure shows ±1.96 MC SE bars. Classifications near the 0.90 threshold are therefore fragile. cr1_G5_icc0.0 (0.912), cr0_G10_icc0.0 (0.897) and cr1_G10_icc0.5 (0.907) are within 2 MC SE of 0.90 (cells_within_2se_of_failure_line); a different seed could move them across the line. The "nominal-consistent" label is a 2-SE screen, not a test with multiplicity control: across 27 pairs, some misclassification is expected. With R = 2000 we cannot resolve coverage differences much smaller than ≈0.02. The relative-bias estimates have MC SE 0.007–0.016 (mc_se_rel_bias_var). Because V̂ is right-skewed at small G, the normal approximation behind that SE is rough, which matches the one 2.46-SE discrepancy from the exact value in G5_icc0.2. The paired SE for the G = 5 ranking uses the fact that all procedures are evaluated on the same replications. The bootstrap critical value adds its own simulation noise (B = 399), and that noise is included in the WCB coverage estimate.

# Limitations and open issues

• Distributions: only Gaussian random effects and errors are covered. Skewed or heavy-tailed u_g/e_ig, which are known to matter for few-cluster inference, are not.
• Heterogeneity: the design has a homoskedastic random intercept and a single fixed size pattern (8,16,32,64). It does not cover heteroskedasticity across clusters, cluster-size–outcome correlation, or other degrees of imbalance. Results may depend strongly on the dominant 64-observation cluster.
• Regressors: only the mean (intercept-only model) is studied. Regression slopes, leverage and fixed effects are not covered.
• Other weights / variants: only Rademacher weights, the unrestricted (non-null-imposed) WCB-t, and the symmetric percentile-t interval with numpy's linear-interpolation quantile are studied. Webb 6-point weights, restricted WCB (WCR), enumeration of all 2^G weight vectors at G = 5, CR2/CR3 (bias-reduced) estimators, and other degrees-of-freedom corrections are not evaluated. Results for the WCB at G = 5 are specific to these choices.
• G ∈ {5, 10, 30} only; the location of the transition between 10 and 30 is not resolved.
• The single 2.46-SE deviation between the MC and exact CR0 bias (G5_icc0.2) is reported, not removed. We read it as MC noise in a skewed estimator (see out/state/checks.md).