# Coverage of three cluster-robust 95% confidence intervals for a mean when there are few clusters: a simulation study (T2-simulation)

Reproduce with `~/wp29/venv/bin/python out/analysis.py`, run from the workspace root. The script reads only `data/design.csv` and writes `out/results.json`, `out/table1.csv` and `out/figure1.png`. It takes about 2 seconds and gives the same output every time it runs. Every number below is quoted from `out/results.json`, and the key it comes from is shown in brackets. A key such as `coverage[cr1_G5_icc0.2]` means the entry `cr1_G5_icc0.2` of the object `coverage`.

## Question
Observations come in G clusters. We estimate the population mean μ of y by the sample mean ȳ. How close to 95% is the actual coverage of three intervals?
- **CR0**: the cluster sandwich variance with a normal critical value of 1.96.
- **CR1**: CR0 multiplied by G/(G−1), with a Student t critical value on G−1 degrees of freedom.
- **WCB**: the unrestricted wild cluster bootstrap-t with Rademacher weights and B = 399 draws.

The design has 9 cells: G ∈ {5, 10, 30} crossed with ICC ∈ {0, 0.2, 0.5}. For each procedure and cell we estimate:
- the coverage;
- its Monte Carlo (MC) standard error;
- for CR0 and CR1, the relative bias of the variance estimator, E[V̂]/Var(ȳ), using the exact Var(ȳ);
- whether the pair is a failure regime (coverage < 0.90).

## Data and exclusions
- **Input.** There is no observed data. The only input is the design file `data/design.csv`, with 9 rows (`n_input_rows`). Its SHA-256 matches the bundle manifest. The audit is in `out/state/data_audit.md`.
- **Design checks.** Cluster sizes follow 8, 16, 32, 64 repeated. In every row the number of sizes equals G and their sum equals N.
- **Settings used.** R = 2000 replications per cell (`R`) and B = 399 bootstrap draws (`B`), for 9 cells (`n_cells`) and 18,000 simulated samples in total (`n_analysed`). Seeds are per cell, from `design.csv`; cell C1 uses 20261101 (`random_seed`).
- **Exclusions.** None. `out/exclusions.csv` records 0 removed, so 18,000 − 0 = 18,000.
- **Design issues found** (`data_issues`):
  1. When ICC = 0, u_g ~ N(0, 0) has zero variance. I draw u_g = √ICC·Z, which is exactly 0.
  2. With G = 5, Rademacher weights allow only 2⁵ = 32 weight vectors. Because flipping every sign leaves |t*| unchanged, there are at most 16 distinct |t*| values, so the bootstrap critical value can take only a few values.
  3. "The 0.95 empirical quantile" is not uniquely defined for finite B. I used numpy's default linear-interpolation quantile.

## Method
**Data-generating process.** In each replication, y_ig = μ + u_g + e_ig, with μ = 0, u_g ~ N(0, ICC) and e_ig ~ N(0, 1−ICC). The script creates one `numpy.random.default_rng(seed)` per cell. Within each replication it draws, in this order: G normals for u, then N normals for e, then a B×G matrix of Rademacher weights.

**Shared quantities.** The residuals are ê = y − ȳ, and S_g is the sum of the residuals in cluster g.
- V̂_CR0 = ΣS_g²/N².
- V̂_CR1 = G/(G−1)·V̂_CR0.

**The three intervals.**
- CR0: ȳ ± 1.96·√V̂_CR0.
- CR1: ȳ ± t_{0.975,G−1}·√V̂_CR1.
- WCB: draw y* = ȳ + w_g·ê. Then ȳ* − ȳ = Σ w_g S_g / N, and the cluster sums of the bootstrap residuals are w_g S_g − n_g(ȳ* − ȳ). From these, compute V̂*_CR1 and t*_b = (ȳ* − ȳ)/√V̂*_CR1. This is algebraically identical to recomputing on the full bootstrap sample, and it is vectorised as one B×G matrix product. The critical value q is the 0.95 quantile of |t*|, and the interval is ȳ ± q·√V̂_CR1.

**Outputs computed.**
- Coverage is the fraction of the R intervals that contain μ. Its MC SE is √(p̂(1−p̂)/R).
- The relative bias is mean(V̂) divided by the exact Var(ȳ) = [ICC·Σn_g² + (1−ICC)N]/N².
- A pair is a failure if p̂ < 0.90. It is consistent with nominal coverage if |p̂ − 0.95| ≤ 2·MC SE.

**Baseline (phase P2).** A naive interval that ignores clustering: ȳ ± 1.96·s/√N (`baseline_iid_coverage`, `baseline_iid_rel_bias_var`).

**Validation.** I derived E[V̂_CR0] analytically as Σ_g Var(T_g − n_g ȳ)/N², where T_g is the cluster total (`exact_rel_bias_var_analytic`), and compared it with the MC estimate. The checks are listed in `out/state/checks.md`.

## Results
Full results are in `out/table1.csv` and `out/figure1.png`. Coverage, with the MC SE in parentheses (`coverage`, `mc_se_coverage`):

| cell | CR0 | CR1 | WCB |
|---|---|---|---|
| G=5, ICC=0.0 | 0.7805 (0.0093) | 0.9115 (0.0064) | 0.8550 (0.0079) |
| G=5, ICC=0.2 | 0.7210 (0.0100) | 0.8845 (0.0071) | 0.8540 (0.0079) |
| G=5, ICC=0.5 | 0.7170 (0.0101) | 0.8755 (0.0074) | 0.8490 (0.0080) |
| G=10, ICC=0.0 | 0.8965 (0.0068) | 0.9435 (0.0052) | 0.9360 (0.0055) |
| G=10, ICC=0.2 | 0.8570 (0.0078) | 0.9155 (0.0062) | 0.9295 (0.0057) |
| G=10, ICC=0.5 | 0.8530 (0.0079) | 0.9065 (0.0065) | 0.9325 (0.0056) |
| G=30, ICC=0.0 | 0.9275 (0.0058) | 0.9400 (0.0053) | 0.9430 (0.0052) |
| G=30, ICC=0.2 | 0.9210 (0.0060) | 0.9350 (0.0055) | 0.9460 (0.0051) |
| G=30, ICC=0.5 | 0.9165 (0.0062) | 0.9300 (0.0057) | 0.9425 (0.0052) |

Relative bias of the variance estimator, E[V̂]/Var(ȳ) (`rel_bias_var`):

| cell | CR0 | CR1 |
|---|---|---|
| G=5, ICC=0.0 / 0.2 / 0.5 | 0.662 / 0.531 / 0.483 | 0.827 / 0.664 / 0.604 |
| G=10, ICC=0.0 / 0.2 / 0.5 | 0.856 / 0.762 / 0.774 | 0.951 / 0.847 / 0.860 |
| G=30, ICC=0.0 / 0.2 / 0.5 | 0.944 / 0.931 / 0.925 | 0.976 / 0.963 / 0.957 |

**Main findings.**
1. **Failure regimes.** 11 of the 27 pairs have coverage below 0.90 (`n_failure_cells`, `failure_cells`):
   - CR0 in all six cells with G ≤ 10;
   - CR1 at G = 5 with ICC = 0.2 and ICC = 0.5;
   - WCB in all three G = 5 cells.
   No procedure fails at G = 30.
2. **CR0 understates the variance substantially when there are few clusters.** At G = 5 its E[V̂]/Var(ȳ) is between 0.48 and 0.66 (`rel_bias_var[cr0_G5_icc*]`). Its coverage is 0.72–0.78 at G = 5 and 0.85–0.90 at G = 10 (`coverage[cr0_*]`).
3. **CR1 is the best of the three at G = 5, but still below nominal there.** Its average coverage over the G = 5 cells is 0.8905, against 0.8527 for WCB and 0.7395 for CR0 (`avg_coverage_small_G`), so `best_procedure_small_G` = "cr1". Even after the G/(G−1) correction, CR1's variance estimate is still biased downward (0.60–0.83 at G = 5).
4. **WCB is the best of the three at G = 10 with ICC > 0, and at G = 30.** Its coverage is 0.9295 and 0.9325 at G = 10 with ICC 0.2 and 0.5, and 0.943–0.946 at G = 30 (`coverage[wcb_*]`).
5. **Only five pairs are consistent with nominal coverage** (`nominal_cells`): CR1 at G = 10 and G = 30 with ICC = 0, and WCB in all three G = 30 cells.
6. **Bias worsens as ICC rises from 0 to 0.2 at small G.** At G = 5 and 10, all three procedures lose coverage when ICC goes from 0 to 0.2, and change little between 0.2 and 0.5. The likely reason is that unequal cluster sizes (8 to 64) let the large clusters dominate Var(ȳ) once there is cluster-level variance.
7. **The baseline shows why a clustered interval is needed.** The naive interval has about 0.94–0.96 coverage at ICC = 0, but only 0.29–0.49 when ICC > 0 (`baseline_iid_coverage`).

## Uncertainty
- **Coverage.** Each coverage estimate is a binomial proportion from R = 2000 replications. Its MC SE is 0.005–0.010 (`mc_se_coverage`); the largest is 0.0101, for CR0 at G = 5 and ICC = 0.5. A 95% MC interval is therefore about ±0.010–0.020 around each estimate.
- **Classification of failures.** Most failure classifications are far from the borderline. The exception is `cr0_G10_icc0.0` (0.8965, SE 0.0068), which is only 0.5 SE below 0.90, so whether it counts as a failure is uncertain.
- **CR1 at G = 5.** `cr1_G5_icc0.0` (0.9115) is above 0.90 by 1.8 SE, but it is still well below 0.95.
- **Relative bias.** The MC SE of the relative-bias estimates is 0.007–0.016 (`mc_se_rel_bias_var`). The MC estimates agree with the analytic E[V̂_CR0]/Var(ȳ) (`exact_rel_bias_var_analytic`) within 2.5 MC SE in all 9 cells. The largest deviation, z = 2.46, is in the G = 5, ICC = 0.2 cell; with 9 comparisons, one deviation of this size is not unusual. When citing CR0/CR1 bias, the analytic values should be preferred, because they have no MC error.
- **Bootstrap.** Using B = 399 adds bootstrap noise to each replication's critical value. Its effect is already included in the WCB coverage estimate. It is not reported separately.
- **Multiple comparisons.** The 27 "consistent with nominal" checks are not adjusted for multiplicity, so roughly one pair could be misclassified by chance.

## Limitations and open issues
- **Distributions.** Only Gaussian, homoskedastic u and e were studied. Skewed or heavy-tailed errors, which are known to hurt few-cluster inference more, are not covered.
- **Heterogeneity.** There is a single fixed pattern of cluster sizes (8/16/32/64) and the ICC is the same across clusters. The design does not vary the degree of size imbalance, and it has no cluster-specific variances or heterogeneous effects. Results may differ with balanced clusters or more extreme imbalance.
- **Regressors.** Only an intercept (the mean) is estimated. Inference on regression slopes, leverage from regressors, treatment indicators that vary only at the cluster level, and few treated clusters are not studied.
- **Other estimators and bootstrap variants.** The following are not evaluated:
  - CR2/CR3 (bias-reduced linearisation) and Bell–McCaffrey or Satterthwaite degrees of freedom;
  - the restricted (null-imposed) wild bootstrap and inverted tests;
  - Webb six-point or Mammen weights;
  - full enumeration of the 2^G weight vectors.
- **Why WCB does poorly at G = 5 is not settled.** With Rademacher weights, the G = 5 bootstrap distribution has at most 16 distinct |t*| values. WCB's poor coverage there is at least partly caused by this discreteness, and Webb weights would probably behave differently. That is not tested here.
- **Choice of quantile rule.** The WCB quantile uses numpy's linear interpolation. Another quantile rule, such as the 380th order statistic of 399, could shift WCB coverage slightly, especially at G = 5. This sensitivity was not run, to stay within the specified design.
- **Scope of the grid.** Only G ∈ {5, 10, 30} and ICC ∈ {0, 0.2, 0.5} were run. Between grid points the behaviour is interpolation, not measurement.