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

Technical note for T2-simulation. To reproduce everything: ~/wp29/venv/bin/python out/analysis.py (about 2 s). It reads only data/design.csv and writes out/results.json, out/table1.csv, out/figure1.png and out/exclusions.csv. Below, keys in code refer to out/results.json. Pairs are named <estimator>_G<G>_icc<icc>.

# Question

When the data are grouped in only a few, unequal-sized clusters, how close to the nominal 95% do three standard cluster-robust intervals for a population mean μ actually cover? The three intervals are CR0 with a normal critical value, CR1 with a t(G−1) critical value, and the wild cluster bootstrap-t with Rademacher weights (WCB). Which (procedure, design) pairs fall below 90% coverage (the "failure regimes")? How biased are the CR0/CR1 variance estimators relative to the exact Var(ȳ)?

# Data and exclusions

No observed data are used. The only input is the design file data/design.csv: 9 rows (n_input_rows), one per cell of G ∈ {5, 10, 30} × ICC ∈ {0, 0.2, 0.5}. Each row gives unequal cluster sizes (8, 16, 32, 64 repeated, so N = 128, 264 and 864), R = 2000, B = 399 and one seed per cell. The SHA-256 hashes match the bundle manifest. In every row the number of listed cluster sizes equals G and their sum equals N; the script asserts both. Nothing was excluded (out/exclusions.csv: 9 → 9 cells). We simulated 9 × 2000 = 18 000 replications (n_analysed, R, B, n_cells). C1 used seed 20261101 (random_seed); the other cells use their own seeds from the design file. One issue is recorded in data_issues: with G = 5 the Rademacher weights have only 32 sign patterns, so the WCB reference distribution is coarse.

# Method

• DGP. y_ig = μ + u_g + e_ig, with u_g ~ N(0, ICC) and e_ig ~ N(0, 1−ICC), and μ = 0. Each cell uses numpy.random.default_rng(seed). Draws are made in a fixed order: u (R×G), then e (R×N), then the Rademacher weights (R×B×G).
• Exact variance. Var(ȳ) = [ICC Σ n_g² + (1−ICC) N] / N² (exact_var_ybar).
• CR0. V̂ = Σ_g (Σ_i ê_ig)² / N², with interval ȳ ± 1.96 √V̂ (1.96 used literally, as specified).
• CR1. V̂ = G/(G−1) · V̂_CR0, with interval ȳ ± t_{0.975, G−1} √V̂.
• WCB. The bootstrap is vectorised. Write S_g for the cluster sums of the residuals. Then ȳ* − ȳ = D = Σ w_g S_g / N, and the bootstrap residuals have cluster sums w_g S_g − n_g D. These give V̂*_CR1 and |t*_b| = |D| / √V̂*_CR1. The critical value q is the 0.95 empirical quantile of the 399 values of |t*|, using the inverted-CDF definition: the ⌈0.95·399⌉ = 380th order statistic, the usual (B+1)(1−α) choice. The interval is ȳ ± q √V̂_CR1. A plain loop implementation gave the same q (max difference 4e−15; see out/state/checks.md).
• Evaluation. Coverage is the share of the R intervals that contain μ (coverage), with MC SE √(p̂(1−p̂)/R) (mc_se_coverage). The variance relative bias is mean(V̂)/Var(ȳ) (rel_bias_var). A pair is a failure when coverage < 0.90 (failure_cells). It is "consistent with nominal" when |p̂ − 0.95| ≤ 2·SE (nominal_cells).
• Baseline (P2). A naive i.i.d. interval, ȳ ± 1.96·s/√N, which ignores clustering, evaluated on the same replications (baseline_coverage_iid).
• Verification. E[V̂_CR0] is a quadratic form, so it was computed analytically: tr(AΣA′)/N² with A = C − n1′/N (exact_rel_bias_var). This was compared with the Monte Carlo means (z_rel_bias_mc_vs_exact).

# Results

Table 1 (out/table1.csv) and Figure 1 (out/figure1.png) give all 27 pairs. Coverage (coverage):

1. Failure regimes (failure_cells, n_failure_cells = 11):

   • CR0 fails in all six cells with G ≤ 10.
   • CR1 fails at G = 5 when ICC = 0.2 or 0.5; at ICC = 0 its coverage is 0.919.
   • WCB fails in all three G = 5 cells.
   • No procedure fails at G = 30.

2. Consistent with nominal 0.95 (nominal_cells): only 6 of the 27 pairs. These are cr1_G10_icc0.0, cr1_G30_icc0.0, cr1_G30_icc0.2 and WCB in all three G = 30 cells. Every other pair under-covers by more than 2 MC SE. Even at G = 30, CR0 (0.916–0.931) and CR1 at ICC = 0.5 (0.929) are significantly below nominal.

3. Best procedure at G = 5 (best_procedure_small_G = cr1). Mean coverage over the three G = 5 cells is 0.893 for CR1, 0.857 for WCB and 0.746 for CR0 (mean_coverage_small_G). At G = 10 with ICC > 0, WCB covers better than CR1 (0.939/0.937 vs 0.917/0.917). At G = 30, WCB is the closest to nominal.

4. WCB at G = 5: wider but covers less. WCB intervals are on average much wider than CR1 (mean half-width 0.88 vs 0.53 at ICC = 0.2; mean_ci_halfwidth), yet they cover less. The bootstrap critical value is very variable: when one large cluster dominates, |t*| takes few distinct values. The result is very wide intervals in some replications and too-narrow ones in others.

5. Variance bias (rel_bias_var). Both sandwich estimators are biased downward in every cell, and the bias is worst with few clusters:

   • CR0: E[V̂]/Var(ȳ) ≈ 0.51–0.65 at G = 5, 0.76–0.85 at G = 10 and 0.92–0.94 at G = 30.
   • CR1: 0.63–0.82, 0.84–0.94 and 0.95–0.97 at the same G.

   The downward bias comes from estimating μ by ȳ, which removes variability that a single dominant cluster carries. It is larger when ICC > 0 than when ICC = 0, because the large clusters then contribute more to Var(ȳ). The analytic values (exact_rel_bias_var) agree.

6. Baseline. Ignoring clustering (baseline_coverage_iid) is fine at ICC = 0 (0.94–0.96) but disastrous at ICC > 0 (0.30–0.48), at every G. All three cluster-robust procedures are a large improvement, but none reaches nominal coverage reliably with 5 clusters.

# Uncertainty

• Coverage. With R = 2000, the MC SE of a coverage estimate is about 0.0049 near 0.95 and at most 0.010 (CR0 at G = 5) (mc_se_coverage). So ±2 SE ≈ ±0.01–0.02.

• Classifications. Classifying a pair as a failure is robust when its coverage is far from 0.90. The cases nearest the threshold are:

  • cr1_G5_icc0.2 (0.8815, SE 0.0072) and cr1_G5_icc0.5 (0.8795, SE 0.0073): about 2.6–2.8 SE below 0.90.
  • cr0_G10_icc0.0 (0.8835, SE 0.0072): about 2.3 SE below.
  • cr1_G5_icc0.0 (0.919, SE 0.0061): about 3.1 SE above.

  Other seeds could plausibly move these by a few hundredths. The "nominal" set rests on a 2-SE rule applied to 27 pairs, so a borderline pair could be misclassified by chance (multiple comparisons; no correction was applied, as the task specifies the rule).

• Variance bias. The MC SE of the relative bias is 0.007–0.016 (mc_se_rel_bias_var). The Monte Carlo means agree with the analytic expectations (exact_rel_bias_var): z-scores range from −2.0 to 1.36 over the 9 independent cells (z_rel_bias_mc_vs_exact). This is consistent with pure MC noise and confirms the implementation.

• Bootstrap error. The WCB coverage also includes bootstrap randomness (B = 399), which is part of the procedure as specified. It is absorbed into the binomial MC SE.

• Comparisons between procedures. All procedures in a cell use the same replications, so their differences are positively correlated. The marginal SEs are therefore conservative for those differences. We did not compute paired SEs separately. The G = 5 gap between CR1 and WCB (0.036 on average) is several times the marginal SE of each mean (≈0.004).

# Limitations and open issues

• Distribution. Only Gaussian random intercepts and Gaussian errors. Skewed or heavy-tailed u_g or e_ig, which usually hurt few-cluster inference more, are not covered.

• Heterogeneity. Within-cluster variance and ICC are the same in every cluster. Cluster sizes follow one fixed unbalanced pattern (8/16/32/64) and are not random or informative about the outcome. Balanced designs, or designs even more dominated by one cluster, may behave differently. At G = 5, half of the sample is in one cluster, which drives much of the small-G behaviour.

• Regressors. The estimand is a simple mean, so there are no regressors and no leverage from covariates. Results for regression coefficients, fixed effects, or treatment indicators that vary only at the cluster level can differ substantially.

• Other methods and weights. We did not study:

  • CR2/CR3 (bias-reduced linearisation) or Bell–McCaffrey/Satterthwaite degrees of freedom;
  • the restricted (null-imposed) WCB or the WCB p-value inversion;
  • Webb six-point or Mammen weights;
  • full enumeration of the 32 Rademacher patterns at G = 5;
  • the equal-tailed rather than symmetric bootstrap-t.

  The WCB results apply only to the unrestricted, symmetric, Rademacher version with B = 399.

• Specification choices.

  • The "0.95 empirical quantile" was implemented as the inverted-CDF (380th of 399) order statistic. Interpolated quantile definitions would change q slightly.
  • CR0 uses the literal 1.96 as specified.

• Design coverage. Only the 9 cells and R = 2000 required by the task were run. No other G values or ICC values were explored.

• Feasibility. The task was completed as specified (feasibility = "complete").