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

## Question
When N observations come from a small number G of unequal-size clusters, how close to 95% is the actual coverage of three cluster-robust intervals for the population mean μ, using the sample mean ȳ as estimator? For each of nine cells (G ∈ {5, 10, 30} × ICC ∈ {0.0, 0.2, 0.5}) we estimate coverage (`coverage`), its Monte Carlo SE (`mc_se_coverage`), the relative bias E[V̂]/Var(ȳ) of the CR0 and CR1 variance estimators (`rel_bias_var`), and the failure regimes where coverage is below 0.90 (`failure_cells`).

## Data and exclusions
There is no observed data. The only input is `data/design.csv` (9 rows, `n_input_rows`; SHA-256 f4ff004e…c50, see `out/state/data_audit.md`). Every row passed the consistency checks (number of sizes = G, sum of sizes = N). No exclusions were applied (`out/exclusions.csv`: 9 → 9 cells). In each cell R = 2000 replications (`R`) and B = 399 bootstrap draws (`B`) were used, giving 18,000 simulated samples (`n_analysed`). Cell seeds are 20261101–20261109; C1's seed is 20261101 (`random_seed`, all seeds under `seeds`). Cluster sizes follow 8, 16, 32, 64 repeated, so with G = 5 the 64-observation cluster contains half the sample. One design-related issue is recorded in `data_issues`: with G = 5, Rademacher weights can produce only 32 distinct weight vectors. No bootstrap draw had a zero variance (`cell_diagnostics.*.n_degenerate_boot` = 0).

## Method
Each replication draws y_ig = μ + u_g + e_ig, with u_g ~ N(0, ICC) and e_ig ~ N(0, 1−ICC), using `numpy.random.default_rng(cell seed)`. Within each replication the draws are made in this order: u, e, then the B×G Rademacher matrix. The exact variance is Var(ȳ) = [ICC·Σn_g² + (1−ICC)N]/N².
- **CR0**: V̂ = Σ_g S_g²/N², where S_g is the sum of residuals y − ȳ in cluster g; interval ȳ ± 1.96√V̂.
- **CR1**: V̂_CR1 = G/(G−1)·V̂_CR0; interval ȳ ± t_{0.975,G−1}√V̂_CR1.
- **WCB**: unrestricted wild cluster bootstrap-t with Rademacher weights. This is vectorised: ȳ* − ȳ = Σ w_g S_g/N, and the bootstrap cluster residual sums are w_g S_g − n_g(ȳ* − ȳ). These give V̂*_CR1 and |t*_b|. The critical value q is the 0.95 empirical quantile of the 399 |t*| values. I took this to be the `inverted_cdf` quantile, which is the 380th order statistic (`wcb_quantile_rule`). The interval is ȳ ± q√V̂_CR1. A brute-force check that builds y* explicitly matched the vectorised version to 1e-12 (`out/state/checks.md`).
- Coverage is the share of replications with |ȳ − μ| ≤ half-width. MC SE = √(p̂(1−p̂)/R). Relative bias is the mean of V̂ over replications divided by the exact Var(ȳ).
- **Baseline (P2)**: a naive i.i.d.-SE z-interval computed on the same replications (`baseline_iid_coverage`).
- Running `out/analysis.py` regenerates `results.json`, `table1.csv` and `figure1.png` in about 2 s. A clean re-run gave byte-identical output.

## Results
Full numbers are in `out/table1.csv` and `out/figure1.png`.
- **Naive baseline**: ignoring clustering gives coverage of about 0.29–0.49 whenever ICC > 0 (e.g. `baseline_iid_coverage.iid_G30_icc0.5` = 0.2985). At ICC = 0 it gives about 0.94–0.96. Some cluster correction is therefore necessary.
- **CR0** under-covers in every cell:
  - G = 5: 0.7805 / 0.7210 / 0.7170 for ICC 0 / 0.2 / 0.5 (`coverage.cr0_G5_icc0.0`, `cr0_G5_icc0.2`, `cr0_G5_icc0.5`);
  - G = 10: 0.8965 / 0.8570 / 0.8530 (`coverage.cr0_G10_*`);
  - G = 30: still only 0.9275 / 0.9210 / 0.9165 (`coverage.cr0_G30_*`).
- **CR0 variance bias**: E[V̂_CR0]/Var(ȳ) is 0.66 / 0.53 / 0.48 at G = 5 (`rel_bias_var.cr0_G5_*`), 0.86 / 0.76 / 0.77 at G = 10 and 0.94 / 0.93 / 0.93 at G = 30. These Monte Carlo values agree with the analytic expectation (`rel_bias_var_exact`) within 2.5 MC SE in all 18 cases.
- **CR1** removes only part of the variance bias: 0.83 / 0.66 / 0.60 at G = 5 (`rel_bias_var.cr1_G5_*`) and 0.96–0.98 at G = 30. Its t critical value adds more width.
  - Coverage at G = 5: 0.9115 / 0.8845 / 0.8755 (`coverage.cr1_G5_*`).
  - Coverage at G = 10: 0.9435 / 0.9155 / 0.9065.
  - Coverage at G = 30: 0.9400 / 0.9350 / 0.9300.
- **WCB** has the best coverage at G = 10 and G = 30 when ICC > 0 (0.9325 and 0.9350 at G = 10, `coverage.wcb_G10_icc0.2`, `wcb_G10_icc0.5`; 0.9445–0.9470 at G = 30). At G = 5 it fails: 0.8575 / 0.8545 / 0.8510 (`coverage.wcb_G5_*`). This happens even though its mean half-width there is the largest of the three procedures (`table1.csv`, `mean_ci_halfwidth`).
- **Failure regimes** (coverage < 0.90): 11 of the 27 pairs (`n_failure_cells`, `failure_cells`):
  - all three CR0 cells at G = 5 and all three at G = 10;
  - all three WCB cells at G = 5;
  - CR1 at G = 5 with ICC 0.2 and 0.5.
- **Best procedure at G = 5**: CR1 has the highest average coverage over the G = 5 cells (0.8905 vs WCB 0.8543 and CR0 0.7395; `best_procedure_small_G`, `mean_coverage_G5`).
- **Consistent with nominal 95%** (|p̂ − 0.95| ≤ 2 MC SE): only 5 pairs (`nominal_cells`): CR1 at G = 10 and G = 30 with ICC = 0, and WCB at G = 30 for all three ICC values.

## Uncertainty
Each coverage estimate is a binomial proportion based on R = 2000 replications. Its MC SE (`mc_se_coverage`) runs from 0.0050 (near 0.95) to 0.0101 (near 0.72). A ±2 SE band near the nominal level is therefore about ±0.010–0.011, and differences between procedures smaller than about 0.015 within a cell should not be over-interpreted. However, all procedures in a cell share the same replications, so paired comparisons are more precise than independent SEs suggest. That pairing was not separately quantified.

The MC SE of the relative-bias estimates (`rel_bias_var_mcse`) is about 0.012–0.014 at G = 5–10 and 0.007–0.009 at G = 30. The largest MC-vs-analytic discrepancy is 2.46 SE (`out/state/checks.md`).

The "nominal" classification uses a 2-SE rule applied separately to 27 pairs, with no multiplicity adjustment. Pairs close to the boundary (e.g. WCB at G = 10, 0.9325–0.9375) could change class under another seed. Similarly, failure classifications near 0.90 (CR0 at G = 10/ICC 0, 0.8965; CR1 at G = 10/ICC 0.5, 0.9065) are within one MC SE of the threshold.

The WCB critical value has its own Monte Carlo noise from B = 399. That noise is included in the coverage estimate but not separated from it.

## Limitations and open issues
- **Distributions**: only Gaussian random effects and errors were simulated. Skewed or heavy-tailed u_g or e_ig, which are known to hurt few-cluster inference, are not covered.
- **Heterogeneity**: the ICC is common to all clusters, and errors are homoskedastic within and across clusters. Heterogeneity tied to cluster size, and informative cluster sizes, are not covered.
- **Cluster-size design**: only one unbalanced pattern (8, 16, 32, 64) was used. With G = 5 a single cluster holds 50% of N, which drives the CR0/CR1 downward bias (high leverage). Balanced designs, or other degrees of imbalance, could give different rankings.
- **Regressors**: the estimand is a plain mean (intercept-only model). Regression coefficients, especially for treatments that vary only at the cluster level or are given to few clusters, are not covered.
- **Other weights and variants**:
  - Only Rademacher weights were used, in the unrestricted (WCU) form. Webb six-point weights, Mammen weights, restricted (null-imposed, WCR) bootstrap and full enumeration of the 2^5 weight vectors at G = 5 were not studied.
  - With G = 5, Rademacher weights give only 16 distinct |t*| values (`data_issues`). Under-coverage of WCB at G = 5 is therefore at least partly a property of this weight choice.
  - CR2/CR3 (Bell–McCaffrey) corrections and other degrees-of-freedom rules were not studied.
- **Quantile convention**: the WCB quantile rule (380th of 399 order statistics) is my choice for the unspecified "empirical quantile". Other interpolation rules would change q slightly.
- **Grid**: only G ∈ {5, 10, 30} and three ICC values were simulated. The G at which each procedure reaches nominal coverage is not pinned down.