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

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

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.