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

Executable note: every number below is produced by ~/wp29/venv/bin/python out/analysis.py (≈2 s), which writes out/results.json, out/table1.csv and out/figure1.png from data/design.csv only. Keys of out/results.json are cited in brackets.

# Question

With N observations in G clusters (G ∈ {5, 10, 30}) and intra-cluster correlation ICC ∈ {0.0, 0.2, 0.5}, how close to 95% is the actual coverage of three cluster-robust intervals for the population mean μ (estimated by ȳ): CR0 with a normal critical value, CR1 with a t(G−1) critical value, and the unrestricted Rademacher wild cluster bootstrap-t (WCB)? We also estimate the relative bias E[V̂]/Var(ȳ) of the CR0/CR1 variance estimators, the Monte Carlo (MC) SE of each coverage, 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…; audit in out/state/data_audit.md): one row per cell with the cluster sizes (pattern 8, 16, 32, 64 repeated; N = 128, 264, 864), R = 2000 [R], B = 399 [B], and a per-cell seed (C1 seed 20261101 [random_seed]). All rows passed consistency checks (G sizes listed, sizes summing to N), so no cell was excluded (out/exclusions.csv: 9 → 9). 9 cells [n_cells] × 2000 replications = 18000 simulated samples [n_analysed]. The design issues found [data_issues]: (i) at G = 5 Rademacher weights have only 32 distinct vectors (16 distinct |t*| values), so the WCB critical value is coarse; (ii) at ICC = 0 the cluster effect is degenerate (zero variance). Neither changed the procedures. Feasibility: complete [feasibility].

# Method

DGP per replication: y_ig = μ + u_g + e_ig, u_g ~ N(0, ICC), e_ig ~ N(0, 1−ICC), μ = 0. Exact Var(ȳ) = [ICC·Σn_g² + (1−ICC)N]/N². For each cell, numpy.random.default_rng(seed); per replication draws in the order u (G), e (N), Rademacher weight matrix (B × G).

• CR0: V̂ = Σ_g(Σ_i ê_ig)²/N², interval ȳ ± 1.96√V̂.
• CR1: V̂_CR1 = G/(G−1)·V̂_CR0, interval ȳ ± t_{0.975,G−1}√V̂_CR1.
• WCB: y* = ȳ + w_g ê_ig, t*_b = (ȳ* − ȳ)/√V̂*_CR1 (V̂* recomputed on the bootstrap sample), q = 0.95 empirical quantile of |t*|, interval ȳ ± q√V̂_CR1. Vectorised via cluster residual sums; checked against a brute-force construction of y* (max difference 1e-14). Empirical-quantile convention: the ⌈(B+1)·0.95⌉ = 380th order statistic of the 399 |t*| values (the standard convention for B = 399); the numpy linear-interpolation quantile is reported as a sensitivity [sensitivity_wcb_coverage_numpy_quantile].
• Coverage = share of R intervals containing μ; MC SE = √(p̂(1−p̂)/R); relative bias = mean(V̂)/exact Var(ȳ), with MC SE sd(V̂)/(√R·Var(ȳ)) [mc_se_rel_bias_var].
• Baseline (P2): the naive i.i.d. interval ȳ ± 1.96·s/√N that ignores clustering [baseline_iid_coverage].

# Results

Full table: out/table1.csv; figure: out/figure1.png (coverage vs G, one panel per ICC, ±1.96 MC SE bars, lines at 0.95 and 0.90).

Coverage [coverage] (MC SE ≈ 0.005–0.010 [mc_se_coverage]):

1. Baseline. Ignoring clustering gives near-nominal coverage only at ICC = 0 and coverage of 0.29–0.49 when ICC > 0 [baseline_iid_coverage]; all three cluster-robust intervals do much better whenever ICC > 0.
2. Failure regimes. 11 of 27 pairs have coverage < 0.90 [n_failure_cells, failure_cells]: CR0 at G = 5 and G = 10 (all ICC; 0.717–0.897), CR1 at G = 5 with ICC = 0.2 and 0.5 (0.885, 0.876), and WCB at G = 5 (all ICC; 0.851–0.858). No pair fails at G = 30.
3. Variance bias. CR0 underestimates Var(ȳ) strongly with few clusters: E[V̂]/Var(ȳ) = 0.48–0.66 at G = 5, 0.76–0.86 at G = 10, 0.93–0.94 at G = 30 [rel_bias_var]. The G/(G−1) factor of CR1 reduces but does not remove the bias at G = 5 (0.60–0.83); at G = 30 CR1 is within 0.96–0.98. The downward bias is larger with ICC > 0, plausibly because the unequal cluster sizes make the large cluster (n = 64) carry much of the variance and its residual sum is shrunk most by centring on ȳ.
4. Comparison of procedures. At G = 5, CR1 has the highest average coverage (0.8905 vs WCB 0.854 and CR0 0.740) [best_procedure_small_G, mean_coverage_G5]. At G = 10 and 30 with ICC > 0, WCB is closest to nominal (0.933–0.947). WCB intervals at G = 5 are much wider on average than CR1 (mean_ci_halfwidth in table1.csv, e.g. 1.39 vs 0.79 at ICC = 0.5) yet cover less often, which points to an unstable, data-dependent critical value when only 16 distinct |t*| values exist.
5. Consistent with nominal (|p̂ − 0.95| ≤ 2 MC SE): only 5 of 27 pairs [nominal_cells]: CR1 at (G=10, ICC=0) and (G=30, ICC=0), and WCB at G = 30 for every ICC. CR0 is never consistent with nominal, not even at G = 30.

# Uncertainty

Each coverage estimate is a binomial proportion with R = 2000; its MC SE [mc_se_coverage] ranges from about 0.005 (coverage near 0.95) to 0.010 (coverage near 0.72), so a 95% MC interval is about ±0.01–0.02. Pairs near the 0.90 threshold are therefore classified with limited confidence: CR0 at (G=10, ICC=0) = 0.8965 (SE 0.0068), CR1 at (G=5, ICC=0) = 0.9115 (SE 0.0064) and CR1 at (G=5, ICC=0.2) = 0.8845 (SE 0.0071) all lie within 2 MC SE of 0.90, so their failure/non-failure label could flip with another seed. The "nominal-consistent" label is a 2-SE screen applied to 27 pairs without multiplicity adjustment, so it is not a formal test. Classifications far from the boundaries (CR0 at G = 5, WCB at G = 5, everything at G = 30 vs 0.90) are robust. The relative-bias estimates have MC SEs of 0.007–0.016 [mc_se_rel_bias_var]. As a check, the MC relative bias of CR0 was compared with its exact analytic expectation (verification script in out/state/verify/): |z| ≤ 1.4 in 8 of 9 cells and z = 2.46 in C2 (G=5, ICC=0.2: MC 0.531 vs exact 0.500), a deviation that fits the somewhat large realised variance of ȳ in that cell; the conclusions about bias do not depend on it. The WCB critical value carries its own bootstrap noise (B = 399), which is part of the procedure's coverage. Switching to numpy's interpolated quantile changes WCB coverage by at most 0.003 and leaves the failure set unchanged [sensitivity_wcb_coverage_numpy_quantile]. Seeds are fixed, so re-running gives identical numbers; a different seed would give MC variation of the size of these SEs.

# Limitations and open issues

• Distributions: only Gaussian cluster effects and errors. Skewed or heavy-tailed u_g/e_ig, which usually hurt few-cluster inference more, are not covered.
• Heterogeneity: only one fixed unequal size pattern (8/16/32/64) and homoskedastic, equal ICC across clusters. Cluster-size-dependent effects, heteroskedasticity across clusters and informative cluster sizes are not studied. Results may depend on the size pattern (e.g. the dominant 64-observation cluster at G = 5).
• Regressors: only the mean (intercept-only model) is studied. Regression slopes, leverage from treated clusters and few treated clusters are outside the design.
• Other weights / variants: only Rademacher weights with the unrestricted (non-null-imposed) WCB-t were studied. Webb six-point or Mammen weights, the restricted WCR bootstrap, full enumeration of the 32 sign vectors at G = 5, CR2/CR3 (bias-reduced linearisation) and Satterthwaite/Bell–McCaffrey degrees of freedom were not evaluated. The G = 5 WCB results should not be generalised to those variants.
• The grid covers only G ∈ {5, 10, 30} and ICC ∈ {0, 0.2, 0.5}; the transition between G = 10 and 30 is not resolved.
• Interval length was summarised only by the mean half-width; the WCB half-width distribution at G = 5 (which seems heavy-tailed) was not analysed further.