# 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`]):

| G | ICC | CR0 | CR1 | WCB | naive iid baseline |
|---|---|---|---|---|---|
| 5 | 0.0 | 0.781 | 0.912 | 0.858 | 0.941 |
| 5 | 0.2 | 0.721 | 0.885 | 0.855 | 0.460 |
| 5 | 0.5 | 0.717 | 0.876 | 0.851 | 0.287 |
| 10 | 0.0 | 0.897 | 0.944 | 0.938 | 0.957 |
| 10 | 0.2 | 0.857 | 0.916 | 0.933 | 0.488 |
| 10 | 0.5 | 0.853 | 0.907 | 0.935 | 0.313 |
| 30 | 0.0 | 0.928 | 0.940 | 0.945 | 0.952 |
| 30 | 0.2 | 0.921 | 0.935 | 0.947 | 0.485 |
| 30 | 0.5 | 0.917 | 0.930 | 0.944 | 0.299 |

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.