# Finite-sample coverage of CR0, CR1 and wild-cluster-bootstrap-t intervals for a mean (T2-simulation)

# Question

When N observations come from a small number G of unequal-size clusters, how close to 95% is the actual coverage of three cluster-robust intervals for the population mean μ, using the sample mean ȳ as estimator? For each of nine cells (G ∈ {5, 10, 30} × ICC ∈ {0.0, 0.2, 0.5}) we estimate coverage (coverage), its Monte Carlo SE (mc_se_coverage), the relative bias E[V̂]/Var(ȳ) of the CR0 and CR1 variance estimators (rel_bias_var), and the failure regimes where coverage is below 0.90 (failure_cells).

# Data and exclusions

There is no observed data. The only input is data/design.csv (9 rows, n_input_rows; SHA-256 f4ff004e…c50, see out/state/data_audit.md). Every row passed the consistency checks (number of sizes = G, sum of sizes = N). No exclusions were applied (out/exclusions.csv: 9 → 9 cells). In each cell R = 2000 replications (R) and B = 399 bootstrap draws (B) were used, giving 18,000 simulated samples (n_analysed). Cell seeds are 20261101–20261109; C1's seed is 20261101 (random_seed, all seeds under seeds). Cluster sizes follow 8, 16, 32, 64 repeated, so with G = 5 the 64-observation cluster contains half the sample. One design-related issue is recorded in data_issues: with G = 5, Rademacher weights can produce only 32 distinct weight vectors. No bootstrap draw had a zero variance (cell_diagnostics.*.n_degenerate_boot = 0).

# Method

Each replication draws y_ig = μ + u_g + e_ig, with u_g ~ N(0, ICC) and e_ig ~ N(0, 1−ICC), using numpy.random.default_rng(cell seed). Within each replication the draws are made in this order: u, e, then the B×G Rademacher matrix. The exact variance is Var(ȳ) = [ICC·Σn_g² + (1−ICC)N]/N².

• CR0: V̂ = Σ_g S_g²/N², where S_g is the sum of residuals y − ȳ in cluster g; interval ȳ ± 1.96√V̂.
• CR1: V̂_CR1 = G/(G−1)·V̂_CR0; interval ȳ ± t_{0.975,G−1}√V̂_CR1.
• WCB: unrestricted wild cluster bootstrap-t with Rademacher weights. This is vectorised: ȳ* − ȳ = Σ w_g S_g/N, and the bootstrap cluster residual sums are w_g S_g − n_g(ȳ* − ȳ). These give V̂*_CR1 and |t*_b|. The critical value q is the 0.95 empirical quantile of the 399 |t*| values. I took this to be the inverted_cdf quantile, which is the 380th order statistic (wcb_quantile_rule). The interval is ȳ ± q√V̂_CR1. A brute-force check that builds y* explicitly matched the vectorised version to 1e-12 (out/state/checks.md).
• Coverage is the share of replications with |ȳ − μ| ≤ half-width. MC SE = √(p̂(1−p̂)/R). Relative bias is the mean of V̂ over replications divided by the exact Var(ȳ).
• Baseline (P2): a naive i.i.d.-SE z-interval computed on the same replications (baseline_iid_coverage).
• Running out/analysis.py regenerates results.json, table1.csv and figure1.png in about 2 s. A clean re-run gave byte-identical output.

# Results

Full numbers are in out/table1.csv and out/figure1.png.

• Naive baseline: ignoring clustering gives coverage of about 0.29–0.49 whenever ICC > 0 (e.g. baseline_iid_coverage.iid_G30_icc0.5 = 0.2985). At ICC = 0 it gives about 0.94–0.96. Some cluster correction is therefore necessary.
• CR0 under-covers in every cell:
  • G = 5: 0.7805 / 0.7210 / 0.7170 for ICC 0 / 0.2 / 0.5 (coverage.cr0_G5_icc0.0, cr0_G5_icc0.2, cr0_G5_icc0.5);
  • G = 10: 0.8965 / 0.8570 / 0.8530 (coverage.cr0_G10_*);
  • G = 30: still only 0.9275 / 0.9210 / 0.9165 (coverage.cr0_G30_*).
• CR0 variance bias: E[V̂_CR0]/Var(ȳ) is 0.66 / 0.53 / 0.48 at G = 5 (rel_bias_var.cr0_G5_*), 0.86 / 0.76 / 0.77 at G = 10 and 0.94 / 0.93 / 0.93 at G = 30. These Monte Carlo values agree with the analytic expectation (rel_bias_var_exact) within 2.5 MC SE in all 18 cases.
• CR1 removes only part of the variance bias: 0.83 / 0.66 / 0.60 at G = 5 (rel_bias_var.cr1_G5_*) and 0.96–0.98 at G = 30. Its t critical value adds more width.
  • Coverage at G = 5: 0.9115 / 0.8845 / 0.8755 (coverage.cr1_G5_*).
  • Coverage at G = 10: 0.9435 / 0.9155 / 0.9065.
  • Coverage at G = 30: 0.9400 / 0.9350 / 0.9300.
• WCB has the best coverage at G = 10 and G = 30 when ICC > 0 (0.9325 and 0.9350 at G = 10, coverage.wcb_G10_icc0.2, wcb_G10_icc0.5; 0.9445–0.9470 at G = 30). At G = 5 it fails: 0.8575 / 0.8545 / 0.8510 (coverage.wcb_G5_*). This happens even though its mean half-width there is the largest of the three procedures (table1.csv, mean_ci_halfwidth).
• Failure regimes (coverage < 0.90): 11 of the 27 pairs (n_failure_cells, failure_cells):
  • all three CR0 cells at G = 5 and all three at G = 10;
  • all three WCB cells at G = 5;
  • CR1 at G = 5 with ICC 0.2 and 0.5.
• Best procedure at G = 5: CR1 has the highest average coverage over the G = 5 cells (0.8905 vs WCB 0.8543 and CR0 0.7395; best_procedure_small_G, mean_coverage_G5).
• Consistent with nominal 95% (|p̂ − 0.95| ≤ 2 MC SE): only 5 pairs (nominal_cells): CR1 at G = 10 and G = 30 with ICC = 0, and WCB at G = 30 for all three ICC values.

# Uncertainty

Each coverage estimate is a binomial proportion based on R = 2000 replications. Its MC SE (mc_se_coverage) runs from 0.0050 (near 0.95) to 0.0101 (near 0.72). A ±2 SE band near the nominal level is therefore about ±0.010–0.011, and differences between procedures smaller than about 0.015 within a cell should not be over-interpreted. However, all procedures in a cell share the same replications, so paired comparisons are more precise than independent SEs suggest. That pairing was not separately quantified.

The MC SE of the relative-bias estimates (rel_bias_var_mcse) is about 0.012–0.014 at G = 5–10 and 0.007–0.009 at G = 30. The largest MC-vs-analytic discrepancy is 2.46 SE (out/state/checks.md).

The "nominal" classification uses a 2-SE rule applied separately to 27 pairs, with no multiplicity adjustment. Pairs close to the boundary (e.g. WCB at G = 10, 0.9325–0.9375) could change class under another seed. Similarly, failure classifications near 0.90 (CR0 at G = 10/ICC 0, 0.8965; CR1 at G = 10/ICC 0.5, 0.9065) are within one MC SE of the threshold.

The WCB critical value has its own Monte Carlo noise from B = 399. That noise is included in the coverage estimate but not separated from it.

# Limitations and open issues

• Distributions: only Gaussian random effects and errors were simulated. Skewed or heavy-tailed u_g or e_ig, which are known to hurt few-cluster inference, are not covered.
• Heterogeneity: the ICC is common to all clusters, and errors are homoskedastic within and across clusters. Heterogeneity tied to cluster size, and informative cluster sizes, are not covered.
• Cluster-size design: only one unbalanced pattern (8, 16, 32, 64) was used. With G = 5 a single cluster holds 50% of N, which drives the CR0/CR1 downward bias (high leverage). Balanced designs, or other degrees of imbalance, could give different rankings.
• Regressors: the estimand is a plain mean (intercept-only model). Regression coefficients, especially for treatments that vary only at the cluster level or are given to few clusters, are not covered.
• Other weights and variants:
  • Only Rademacher weights were used, in the unrestricted (WCU) form. Webb six-point weights, Mammen weights, restricted (null-imposed, WCR) bootstrap and full enumeration of the 2^5 weight vectors at G = 5 were not studied.
  • With G = 5, Rademacher weights give only 16 distinct |t*| values (data_issues). Under-coverage of WCB at G = 5 is therefore at least partly a property of this weight choice.
  • CR2/CR3 (Bell–McCaffrey) corrections and other degrees-of-freedom rules were not studied.
• Quantile convention: the WCB quantile rule (380th of 399 order statistics) is my choice for the unspecified "empirical quantile". Other interpolation rules would change q slightly.
• Grid: only G ∈ {5, 10, 30} and three ICC values were simulated. The G at which each procedure reaches nominal coverage is not pinned down.