# Finite-sample coverage of cluster-robust confidence intervals for a mean (T2-simulation)

All numbers below come from out/analysis.py (run: ~/wp29/venv/bin/python out/analysis.py from the workspace root; ~2 s; deterministic given the seeds in data/design.csv). Full numbers: out/table1.csv, out/results.json; figure: out/figure1.png.

# Question

With few clusters, how close to 95% is the actual coverage of three cluster-robust intervals for a population mean μ — CR0 with a normal critical value, CR1 with a t(G−1) critical value, and the wild cluster bootstrap-t (WCB; Rademacher, unrestricted) — 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

There are no observed data. The only input is data/design.csv (9 rows = 9 cells: G ∈ {5,10,30} × ICC ∈ {0,0.2,0.5}). Cluster sizes repeat 8,16,32,64 (so cluster sizes are unbalanced, ranging from 8 to 64). The script checks that each cluster_sizes list has G entries and sums to N; both checks pass. Nothing is excluded. R = 2000 replications and B = 399 bootstrap draws per replication in each cell, so 18,000 replications in total.

# Method

• DGP: y_ig = μ + u_g + e_ig, u_g ~ N(0, ICC), e_ig ~ N(0, 1−ICC), μ = 0. Each cell uses its own numpy.random.default_rng(seed). Draws are made in a fixed order: U (R×G), then E (R×N), then the Rademacher weights (R×B×G).
• CR0: V̂ = Σ_g(Σ_i ê_ig)²/N², interval ȳ ± 1.96√V̂. CR1: V̂·G/(G−1), interval ȳ ± t_{0.975,G−1}√V̂_CR1.
• WCB: y* = ȳ + w_g ê_ig. The bootstrap is computed exactly through cluster residual sums S_g. The bootstrap estimate shifts by D = Σ w_g S_g / N, and the bootstrap cluster residual sums are w_g S_g − n_g D. From these, t* = D/√V̂*_CR1 and q = the 0.95 quantile of |t*| over B = 399 draws. The quantile uses numpy.quantile's default linear-interpolation rule; this is a documented choice, since the task does not fix a quantile definition. The interval is ȳ ± q√V̂_CR1.
• Coverage is the share of the R intervals that contain μ, with MC SE = √(p̂(1−p̂)/R). Relative bias is mean(V̂)/Var(ȳ), using the exact Var(ȳ) from TASK.md. A pair is a failure if p̂ < 0.90. It is "consistent with nominal" if |p̂ − 0.95| ≤ 2·MC SE.
• Implementation checks: the empirical variance of ȳ across replications matches the exact Var(ȳ) in every cell (e.g. C1: 0.00758 vs 0.00781). I also derived E[V̂_CR0] in closed form. With T_g = Σ_i(y_ig−μ), a_g = Var(T_g) = ICC·n_g² + (1−ICC)·n_g, and S_g = T_g − (n_g/N)ΣT, the expectation is E[S_g²] = a_g(1 − 2n_g/N) + (n_g/N)²Σa. The simulated relative biases agree with this formula to within 1.4 MC SE in every cell (e.g. C1 CR0: simulated 0.654 vs analytic 0.664; C3: 0.505 vs 0.488; C4: 0.853 vs 0.839). The analytic values are stored in results.json → supplementary.cell_diagnostics.

# Results

Coverage of the nominal 95% intervals (MC SE ≈ 0.005–0.010):

• Failure regimes (coverage < 0.90): 11 of 27 pairs. They are CR0 at G = 5 and G = 10 (all ICCs), CR1 at G = 5 with ICC 0.2 and 0.5, and WCB at G = 5 (all ICCs). No pair fails at G = 30.
• Consistent with nominal (|p̂−0.95| ≤ 2 MC SE): 5 pairs. These are CR1 at G=30 with ICC 0 and 0.2, and WCB at G=30 for all three ICCs. Every other pair, including CR0 everywhere, undercovers by a statistically clear margin.
• Variance bias: CR0 underestimates Var(ȳ) everywhere. It recovers only about 50% of the variance at G=5 with ICC>0, about 76% at G=10, and about 92% at G=30. Multiplying by G/(G−1) (CR1) helps but does not remove the bias at small G (0.63 at G=5, ICC>0). The downward bias is larger with ICC>0 than with ICC=0. That is because, under clustering, the variance is dominated by the largest clusters (size 64), whose own residual sums are shrunk the most when the overall mean is subtracted.
• Best procedure at G=5: CR1, with mean coverage 0.894 over the three G=5 cells, compared with WCB 0.856 and CR0 0.747. At G=10 with ICC>0, WCB is closest to nominal (0.936/0.935 vs CR1 0.917). At G=30, WCB and CR1 are both close to nominal, and CR0 still undercovers by 1.5–3.4 points.
• The WCB intervals at G=5 are on average wider than the CR1 intervals (e.g. mean half-width 0.88 vs 0.53 at ICC 0.2) yet cover less often. Its bootstrap critical value is larger on average (mean q ≈ 3.1–4.1 vs t_{0.975,4} = 2.78) but varies from replication to replication. The extra width therefore seems to go to samples that would have been covered anyway. I did not check this mechanism directly, e.g. via the joint distribution of q and V̂; it is a conjecture.

# Uncertainty

• Coverage: each coverage estimate is a binomial proportion over R = 2000 independent replications, with MC SE between 0.0047 (near 0.95) and 0.0100 (near 0.72). A 95% MC interval is therefore about ±0.01–0.02. Classifications near a threshold are fragile. The closest calls are CR1 G=5 ICC=0 (0.920, 3.3 SE above 0.90), CR1 G=5 ICC=0.5 (0.880, 2.8 SE below), and CR1 G=30 ICC=0.5 (0.929, about 3.7 SE below 0.95, so not nominal). WCB/CR1 differences within a cell are computed on the same simulated samples, so they are positively correlated; the marginal SEs therefore overstate the uncertainty of those differences.
• Relative bias: the MC SE of mean(V̂)/Var(ȳ) is about 0.012–0.016 at G=5/10 and about 0.007–0.009 at G=30 (stored in supplementary.rel_bias_var_mc_se). The closed-form check agrees within 1.4 MC SE in all 18 pairs.
• Bootstrap error: with B = 399, q itself is noisy. At G=5, Rademacher weights take only 2^5 = 32 distinct values, i.e. 16 distinct |t*| values, so q is computed from heavily repeated values. This is a property of the specified design, not an estimation error, and it is recorded in data_issues.
• Seeds: one seed per cell, as in the design. Results for other seeds would differ by roughly the MC SE.

# Limitations and open issues

• Distributions: the design uses Gaussian cluster effects and errors only. Skewed or heavy-tailed u_g/e_ig, which typically worsen small-G performance, are not covered.
• Heterogeneity: ICC and error variance are the same in every cluster. One fixed size pattern (8–64) is used, so size–effect correlation and heteroskedasticity across clusters are absent. Other degrees of size imbalance are not explored.
• Regressors: only an intercept-only mean is studied. Regression slopes, treated-cluster designs (e.g. few treated clusters in DiD), and fixed effects are not covered. Leverage-adjusted CR2/CR3 (Bell–McCaffrey) variants are not included.
• Other weights/variants: only unrestricted Rademacher WCB-t with a symmetric |t| quantile is studied. The restricted (null-imposed) WCB, Webb six-point or Mammen weights (relevant at G=5 because Rademacher gives only 32 distinct draws), percentile-t equal-tailed intervals, and alternative quantile definitions are not studied.
• The grid is coarse (G ∈ {5,10,30}, three ICCs), so the G at which each procedure becomes reliable is located only roughly.