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

All numbers below are produced by `out/analysis.py` (run: `~/wp29/venv/bin/python out/analysis.py`, ≈2 s) and are cited by their key in `out/results.json`. Full per-cell values are in `out/table1.csv`; the figure is `out/figure1.png`.

## Question
When N observations are grouped in few (G = 5, 10, 30) unbalanced clusters with intra-cluster correlation ICC ∈ {0, 0.2, 0.5}, what is the actual coverage of three nominal-95% cluster-robust intervals for the mean μ (CR0 with z critical value; CR1 with t_{G−1}; unrestricted Rademacher wild cluster bootstrap-t, WCB), how biased are the CR0/CR1 variance estimators relative to the exact Var(ȳ), and in which (procedure, cell) pairs does coverage fall below 0.90?

## Data and exclusions
No observed data. The only input is `data/design.csv` (`n_input_rows`: 9 rows; SHA-256 recorded in `out/state/data_audit.md`): 9 cells, cluster sizes 8,16,32,64 repeated (N = 128, 264, 864), R = 2000 (`R`), B = 399 (`B`), one seed per cell (C1 seed = `random_seed`; all seeds in `seeds_by_cell`). The audit found no missing values, duplicates or inconsistencies (sizes count = G, sum = N). No exclusions: all `n_cells` = 9 cells × 2000 replications = `n_analysed` = 18000 replications were simulated and analysed (`out/exclusions.csv`, one no-op step). `feasibility` = "complete". One design feature is recorded in `data_issues`: at G = 5 Rademacher weights take only 32 distinct values (16 distinct |t*|), so the bootstrap critical value is coarse; kept as specified. No bootstrap draw had a non-positive variance (`n_degenerate_bootstrap_draws` = 0).

## Method
Per replication (numpy `default_rng(cell seed)`, draws in fixed order u, e, W): y_ig = μ + u_g + e_ig, u_g ~ N(0, ICC), e_ig ~ N(0, 1−ICC), μ = 0. With ê = y − ȳ and S_g = Σ_i ê_ig: V̂_CR0 = Σ S_g²/N², CI ȳ ± 1.96√V̂_CR0; V̂_CR1 = G/(G−1)·V̂_CR0, CI ȳ ± t_{0.975,G−1}√V̂_CR1; WCB: B×G Rademacher matrix W, ȳ*−ȳ = Σ_g w_g S_g/N, bootstrap cluster residual sums w_g S_g − n_g(ȳ*−ȳ), V̂*_CR1 from these, t*_b = (ȳ*−ȳ)/√V̂*_CR1, q = 0.95 quantile of |t*| (numpy linear interpolation, `wcb_quantile_method`), CI ȳ ± q√V̂_CR1. Coverage = share of R intervals containing μ; MC SE = √(p̂(1−p̂)/R); relative bias = mean(V̂)/Var(ȳ) with Var(ȳ) = [ICC Σn_g² + (1−ICC)N]/N². Failure regime: coverage < 0.90; nominal-consistent: |p̂ − 0.95| ≤ 2 MC SE. Baseline (P2): the naive i.i.d. interval ȳ ± 1.96·s/√N on the same draws (`baseline_iid_coverage`).

Verification (`out/state/checks.md`): byte-identical re-run; vectorised WCB matches a brute-force construction of y* (diff 7e-15); MC relative bias of CR0 matches the closed-form E[V̂_CR0]/Var(ȳ) within 2 MC SE in 8/9 cells (z = +2.46 in G5_icc0.2); MC Var(ȳ)/exact in [0.93, 1.06] (`check_mc_var_ybar_over_exact`).

## Results
Coverage (`coverage`, MC SE in `mc_se_coverage`):

| G | ICC | CR0 | CR1 | WCB | naive i.i.d. baseline |
|---|---|---|---|---|---|
| 5 | 0.0 | 0.781 | 0.912 | 0.855 | 0.941 |
| 5 | 0.2 | 0.721 | 0.885 | 0.854 | 0.460 |
| 5 | 0.5 | 0.717 | 0.876 | 0.849 | 0.287 |
| 10 | 0.0 | 0.897 | 0.944 | 0.936 | 0.957 |
| 10 | 0.2 | 0.857 | 0.916 | 0.930 | 0.488 |
| 10 | 0.5 | 0.853 | 0.907 | 0.933 | 0.313 |
| 30 | 0.0 | 0.928 | 0.940 | 0.943 | 0.952 |
| 30 | 0.2 | 0.921 | 0.935 | 0.946 | 0.485 |
| 30 | 0.5 | 0.917 | 0.930 | 0.943 | 0.299 |

- **Failure regimes** (`failure_cells`, `n_failure_cells` = 11): CR0 at every G = 5 and G = 10 cell (6 pairs); WCB at all three G = 5 cells; CR1 at G = 5 with ICC = 0.2 and 0.5. No procedure fails at G = 30.
- **Nominal-consistent pairs** (`nominal_cells`, 5 pairs): cr1_G10_icc0.0, cr1_G30_icc0.0, wcb_G30_icc0.0, wcb_G30_icc0.2, wcb_G30_icc0.5. CR0 is never consistent with 0.95, even at G = 30 (0.917–0.928).
- **Best at G = 5** (`best_procedure_small_G` = "cr1"): average G = 5 coverage 0.891 for CR1, 0.853 for WCB, 0.740 for CR0 (`avg_coverage_G5`). The CR1 − WCB difference is 0.038 with paired MC SE 0.0035 (`avg_coverage_G5_paired_diff`), so the ranking is not a Monte Carlo artefact. At G = 10 with ICC > 0, WCB covers better than CR1 (0.930/0.933 vs 0.916/0.907).
- **Variance bias** (`rel_bias_var`): both sandwich estimators are biased downward, most at small G and high ICC. CR0: 0.66, 0.53, 0.48 at G = 5 (ICC 0, 0.2, 0.5); 0.86, 0.76, 0.77 at G = 10; 0.94, 0.93, 0.93 at G = 30. CR1: 0.83, 0.66, 0.60 at G = 5; 0.95, 0.85, 0.86 at G = 10; 0.98, 0.96, 0.96 at G = 30. The unbalanced sizes (one 64-observation cluster dominates) make the bias much worse than the balanced-cluster factor (G−1)/G.
- **Baseline**: the naive i.i.d. interval is fine at ICC = 0 but covers only 0.29–0.49 when ICC > 0 (`baseline_iid_coverage`), whatever G is. All three cluster-robust procedures improve on it there.
- Mean CI half-widths are in `table1.csv` (`mean_ci_halfwidth`). At G = 5 the WCB interval is on average *wider* than CR1 but still covers less. One plausible explanation (not tested here) is that q varies a lot across replications, because the bootstrap distribution has only 16 distinct |t*| values.

## Uncertainty
Each coverage is a binomial proportion over R = 2000 independent replications. Its MC SE (`mc_se_coverage`) is 0.005–0.010, so a 95% MC interval is about ±0.010–0.020. The figure shows ±1.96 MC SE bars. Classifications near the 0.90 threshold are therefore fragile. cr1_G5_icc0.0 (0.912), cr0_G10_icc0.0 (0.897) and cr1_G10_icc0.5 (0.907) are within 2 MC SE of 0.90 (`cells_within_2se_of_failure_line`); a different seed could move them across the line. The "nominal-consistent" label is a 2-SE screen, not a test with multiplicity control: across 27 pairs, some misclassification is expected. With R = 2000 we cannot resolve coverage differences much smaller than ≈0.02. The relative-bias estimates have MC SE 0.007–0.016 (`mc_se_rel_bias_var`). Because V̂ is right-skewed at small G, the normal approximation behind that SE is rough, which matches the one 2.46-SE discrepancy from the exact value in G5_icc0.2. The paired SE for the G = 5 ranking uses the fact that all procedures are evaluated on the same replications. The bootstrap critical value adds its own simulation noise (B = 399), and that noise is included in the WCB coverage estimate.

## Limitations and open issues
- **Distributions**: only Gaussian random effects and errors are covered. Skewed or heavy-tailed u_g/e_ig, which are known to matter for few-cluster inference, are not.
- **Heterogeneity**: the design has a homoskedastic random intercept and a single fixed size pattern (8,16,32,64). It does not cover heteroskedasticity across clusters, cluster-size–outcome correlation, or other degrees of imbalance. Results may depend strongly on the dominant 64-observation cluster.
- **Regressors**: only the mean (intercept-only model) is studied. Regression slopes, leverage and fixed effects are not covered.
- **Other weights / variants**: only Rademacher weights, the unrestricted (non-null-imposed) WCB-t, and the symmetric percentile-t interval with numpy's linear-interpolation quantile are studied. Webb 6-point weights, restricted WCB (WCR), enumeration of all 2^G weight vectors at G = 5, CR2/CR3 (bias-reduced) estimators, and other degrees-of-freedom corrections are not evaluated. Results for the WCB at G = 5 are specific to these choices.
- G ∈ {5, 10, 30} only; the location of the transition between 10 and 30 is not resolved.
- The single 2.46-SE deviation between the MC and exact CR0 bias (G5_icc0.2) is reported, not removed. We read it as MC noise in a skewed estimator (see `out/state/checks.md`).