# 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(ȳ):

| G | ICC | CR0 cov | CR1 cov | WCB cov | CR0 rel. bias | CR1 rel. bias |
|---|---|---|---|---|---|---|
| 5 | 0.0 | 0.800 (0.009) | 0.920 (0.006) | 0.870 (0.008) | 0.654 | 0.818 |
| 5 | 0.2 | 0.720 (0.010) | 0.882 (0.007) | 0.847 (0.008) | 0.507 | 0.634 |
| 5 | 0.5 | 0.721 (0.010) | 0.880 (0.007) | 0.857 (0.008) | 0.505 | 0.632 |
| 10 | 0.0 | 0.877 (0.007) | 0.934 (0.006) | 0.929 (0.006) | 0.853 | 0.948 |
| 10 | 0.2 | 0.859 (0.008) | 0.917 (0.006) | 0.939 (0.005) | 0.760 | 0.844 |
| 10 | 0.5 | 0.853 (0.008) | 0.917 (0.006) | 0.937 (0.005) | 0.767 | 0.852 |
| 30 | 0.0 | 0.935 (0.006) | 0.951 (0.005) | 0.957 (0.005) | 0.940 | 0.972 |
| 30 | 0.2 | 0.924 (0.006) | 0.940 (0.005) | 0.952 (0.005) | 0.919 | 0.950 |
| 30 | 0.5 | 0.916 (0.006) | 0.929 (0.006) | 0.942 (0.005) | 0.915 | 0.947 |

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.