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

All numbers below are produced by out/analysis.py (run ~/wp29/venv/bin/python out/analysis.py from the workspace root; ~2 s runtime, deterministic given the seeds in data/design.csv). They are in out/results.json, out/table1.csv and out/cell_diagnostics.csv; the figure is out/figure1.png.

# Question

When only a few clusters are available, how close to 95% is the actual coverage of three cluster-robust intervals for the population mean μ, estimated by the sample mean ȳ: CR0 with a normal critical value, CR1 with a t(G−1) critical value, and the wild cluster bootstrap-t (WCB, Rademacher, unrestricted)? Also: how biased are the CR0/CR1 variance estimators relative to the exact Var(ȳ), and in which (procedure, cell) pairs does coverage drop below 0.90?

# Data and exclusions

No observed data are used. The only input is data/design.csv, which has 9 rows (cells C1–C9: G ∈ {5,10,30} × ICC ∈ {0,0.2,0.5}). Cluster sizes repeat the pattern 8,16,32,64, so N = 128, 264 and 864. The script checked that the listed cluster sizes match G and N in every row, and they all do. Nothing was excluded. Each cell uses R = 2000 replications and B = 399 bootstrap draws, so 18,000 replications were simulated in total. No bootstrap draw produced a zero bootstrap variance.

# Method

• Data-generating process (DGP): y_ig = μ + u_g + e_ig, with u_g ~ N(0, ICC) and e_ig ~ N(0, 1−ICC). Each cell uses numpy.random.default_rng(seed) with that cell's seed. All R replications are drawn in one go in this order: u (R×G; skipped when ICC = 0, where u ≡ 0), then e (R×N), then the Rademacher weights (R×B×G).
• CR0: V̂ = Σ_g(Σ_i ê_ig)²/N²; the interval is ȳ ± 1.96√V̂.
• CR1: V̂ = G/(G−1)·V̂_CR0; the interval is ȳ ± t_{0.975,G−1}√V̂.
• WCB: y* = ȳ + w_g ê_ig. Every bootstrap quantity depends only on the cluster residual sums S_g:
  • ȳ* − ȳ = Σ w_g S_g / N;
  • the bootstrap residual cluster sums are w_g S_g − n_g(ȳ* − ȳ);
  • t* = (ȳ* − ȳ)/√V̂*_CR1. This is fully vectorised and algebraically identical to recomputing everything on y*.
  • The critical value q is the "0.95 empirical quantile of |t*|". It was implemented as the ⌈0.95(B+1)⌉ = 380th order statistic of the 399 values of |t*|, the usual bootstrap convention (documented choice). The interval is ȳ ± q√V̂_CR1.
• Coverage, MC SE and relative bias: coverage is the share of R intervals that contain μ = 0, and its MC SE is √(p̂(1−p̂)/R). Relative bias is mean(V̂)/Var(ȳ), using the exact Var(ȳ) = [ICC·Σn_g² + (1−ICC)N]/N².
• Classification: a pair is a failure when p̂ < 0.90. It is "consistent with nominal" when |p̂ − 0.95| ≤ 2·MC SE.
• Sanity check: in every cell, the Monte Carlo variance of ȳ is within sampling error of the exact Var(ȳ). The ratios run from 0.93 to 1.04; with R = 2000 the sampling SD of this ratio is about 0.03.

# Results

Coverage (MC SE in parentheses) and relative variance bias E[V̂]/Var(ȳ):

Main findings:

1. Failure regimes (11 of 27 pairs below 0.90):
   • CR0 fails in all G = 5 and all G = 10 cells (coverage 0.72–0.88).
   • CR1 fails at G = 5 with ICC = 0.2 and 0.5 (about 0.88).
   • WCB fails in all three G = 5 cells (0.85–0.87).
   • No procedure fails at G = 30.
2. Consistent with nominal coverage: only 5 pairs, all at G = 30: CR1 at ICC 0 and 0.2, and WCB at ICC 0, 0.2 and 0.5. CR0 is never consistent with nominal, not even at G = 30, where it under-covers by 1.5–3.4 points.
3. Small G: at G = 5, CR1 has the highest average coverage (best_procedure_small_G = cr1; mean 0.894, compared with 0.858 for WCB and 0.747 for CR0). At G = 10, WCB is closest to nominal when ICC > 0 (0.937–0.939, compared with 0.917 for CR1).
4. Variance bias: the sandwich variance is biased downward in every cell.
   • At G = 5 with ICC > 0, CR0 captures only about half of Var(ȳ) and CR1 about 63%.
   • The bias shrinks with G. At G = 30, CR1 reaches 0.95–0.97, CR0 0.92–0.94.
   • The bias is worse with ICC > 0 than with ICC = 0. Because cluster sizes are unequal (8 to 64), the large clusters dominate Var(ȳ) when ICC > 0, which makes the effective number of clusters smaller than G.
5. WCB at G = 5: WCB under-covers even though its average half-width is larger than CR1's (for example 0.88 vs 0.53 at ICC = 0.2), and its average critical value (3.1–4.2) is above t_{0.975,4} = 2.78. So the bootstrap q is large on average but not where it is needed. A plausible explanation, which I did not test separately, involves two features:
   • With G = 5 there are only 2⁵ = 32 Rademacher weight vectors (16 distinct |t*| values), so the bootstrap distribution is very coarse.
   • q tends to be small in exactly the samples where one large cluster drives a large |t|.

# Uncertainty

• Coverage: each coverage estimate is a binomial proportion from R = 2000 independent replications. Its MC SE ranges from 0.0046 (coverage near 0.96) to 0.010 (coverage near 0.72), and is reported for every pair in table1.csv and results.json. A 95% MC interval is about ±0.01–0.02.
• Robust vs borderline conclusions: the failure classifications at G = 5 for CR0 and WCB, and for CR0 at G = 10, are many SEs below 0.90, so they are robust. Three are borderline:
  • CR1 at G = 5 with ICC 0.2 and 0.5 (0.882 and 0.880) is about 2.5–3 SE below 0.90.
  • CR1 at G = 5 with ICC 0 (0.920) is about 3 SE above it.
  • The 2·MC SE "nominal" band is a pointwise criterion that is not adjusted for multiplicity across 27 pairs. Pairs near its edge, such as CR1 at G = 30 with ICC 0.2 (0.940, SE 0.0053) or CR1 at G = 10 with ICC 0 (0.934), could be classified differently with another seed.
• Variance bias: the MC SE of the relative-bias estimates is sd(V̂)/(√R·Var(ȳ)). It is about 0.012–0.016 at G = 5 and 10, and about 0.007–0.009 at G = 30 (cell_diagnostics.csv). All the downward biases reported above are therefore far beyond Monte Carlo noise, except the difference between CR1 at G = 30 and 1, which is 3–6 SE.
• Bootstrap error: the WCB critical values carry extra simulation noise from using B = 399 rather than infinitely many draws. This noise is already included in the binomial MC SE of WCB coverage, but it means the results describe WCB with B = 399, not the ideal bootstrap.

# Limitations and open issues

• Distributions: only Gaussian cluster effects and errors were studied. Skewed or heavy-tailed u_g or e_ig, which typically hurt small-G inference more, are not covered.
• Heterogeneity: variances and ICC are the same in every cluster, with no correlation between cluster size and cluster effect. Only one pattern of cluster-size imbalance was used (8,16,32,64 repeated), and the G = 5 design has two clusters of size 8. Balanced clusters, or stronger imbalance, could change the ranking.
• Regressors: the estimand is a plain mean. Regression coefficients, few treated clusters and leverage effects (where CR2/CR3 and restricted WCB matter) are not studied.
• Other weights and variants: WCB used only Rademacher weights, unrestricted (no null imposed), with B = 399 random draws. Webb 6-point or Mammen weights, null-restricted WCB (WCR), full enumeration of the 32 weight vectors at G = 5, CR2/Bell–McCaffrey degrees of freedom, and CR3/jackknife were not evaluated.
• Other choices:
  • The bootstrap quantile was taken as the 380th order statistic of 399. Other quantile definitions differ slightly.
  • CR0 uses 1.96 as specified rather than the exact 1.95996.
• Unchecked explanation: the explanation offered for WCB's behaviour at G = 5 has not been checked by a separate analysis.
• Design limits: only G ∈ {5,10,30} and ICC ∈ {0,0.2,0.5} were studied. The transition between G = 10 and 30 is not resolved.