# T2-simulation — Finite-sample performance of three cluster-robust confidence intervals for a mean: a simulation study

# Scientific question

When observations are grouped in a small number of clusters, how reliable are the usual cluster-robust confidence intervals for a population mean? This is a complete simulation-methodology project: design → implementation → Monte Carlo evaluation with its own uncertainty → identification of failure regimes → figure and table → an executable technical note. The estimand of the inferential procedures is the population mean μ of y, estimated by the sample mean ȳ of N clustered observations; the estimands of the study are, for each of the nine design cells (G ∈ {5, 10, 30} clusters × intra-cluster correlation ICC ∈ {0.0, 0.2, 0.5}) and each of three interval procedures, (a) the actual coverage probability of the nominal 95% interval, (b) for the two sandwich estimators, the relative bias of the variance estimator, E[V̂] / Var(ȳ), where Var(ȳ) is the exact variance of the sample mean under the design, (c) the Monte Carlo standard error of each coverage estimate, and (d) the set of (procedure, cell) pairs whose coverage falls below 0.90 — the failure regimes.

# Inputs

data/design.csv — simulated-data design written by the task designer (nothing here is observed data): one row per cell with cell_id, G, icc, cluster_sizes (the G cluster sizes separated by |: the pattern 8, 16, 32, 64 repeated), N, R = 2000 replications, B = 399 bootstrap draws, seed, mu = 0, total_variance = 1, nominal_level = 0.95. See data/README.md.

# Information-availability rules

Data-generating process (per replication): y_ig = μ + u_g + e_ig with u_g ~ N(0, ICC) i.i.d. across clusters and e_ig ~ N(0, 1 − ICC) i.i.d. across observations, independent of u; cluster g has n_g observations as listed in cluster_sizes. Exact variance of the sample mean: Var(ȳ) = [ICC · Σ_g n_g² + (1 − ICC) · N] / N².
Interval procedures (all three required, exactly as specified, with ȳ as the point estimate and residuals ê_ig = y_ig − ȳ):
- CR0: V̂_CR0 = Σ_g (Σ_i ê_ig)² / N²; interval ȳ ± 1.96 · √V̂_CR0;
- CR1: V̂_CR1 = [G / (G − 1)] · V̂_CR0; interval ȳ ± t_{0.975, G−1} · √V̂_CR1 (Student t with G − 1 degrees of freedom);
- WCB (wild cluster bootstrap-t, unrestricted, Rademacher): for b = 1..B draw cluster weights w_g^b ∈ {−1, +1} with equal probability, form y*_ig = ȳ + w_g^b · ê_ig, compute ȳ* and its V̂_CR1* on the bootstrap sample, t*_b = (ȳ* − ȳ) / √V̂_CR1*; critical value q = the 0.95 empirical quantile of |t*_1|, …, |t*_B|; interval ȳ ± q · √V̂_CR1.
Use R = 2000 replications per cell, B = 399 bootstrap draws per replication and the cell's seed from design.csv; do not change the grid, R or B. A procedure's coverage in a cell is the fraction of the R intervals that contain μ. The Monte Carlo standard error of a coverage estimate p̂ is √(p̂(1 − p̂)/R).
A (procedure, cell) pair is a failure regime when its estimated coverage is below 0.90. A pair is consistent with nominal coverage when |p̂ − 0.95| ≤ 2 · MC SE.

# Required outputs

out/table1.csv — 27 rows (3 procedures × 9 cells) with columns estimator (cr0, cr1, wcb), cell_id, G, icc, coverage, mc_se_coverage, rel_bias_var (empty for wcb), mean_ci_halfwidth, R.
out/results.json with the keys:

key	type	unit	meaning
`task_id`	string	—	"T2-simulation"
`R`	integer	—	replications per cell actually used
`B`	integer	—	bootstrap draws per replication actually used
`n_cells`	integer	—	number of design cells simulated
`coverage`	object	—	{"_G_icc": number} for the 27 pairs, e.g. `"cr1_G5_icc0.2"` (icc written as 0.0, 0.2, 0.5)
`mc_se_coverage`	object	—	same keys: Monte Carlo SE of each coverage estimate
`rel_bias_var`	object	—	18 keys (cr0 and cr1 only): E[V̂]/Var(ȳ) estimated by the Monte Carlo mean
`failure_cells`	array of strings	—	the pairs with coverage < 0.90, same key format, any order
`n_failure_cells`	integer	—	length of `failure_cells`
`nominal_cells`	array of strings	—	the pairs whose coverage is consistent with nominal (
`best_procedure_small_G`	string	—	the procedure with the highest average coverage over the three G = 5 cells
`feasibility`	string	—	"complete", "partial" or "infeasible"
`feasibility_reason`	string or null	—	explanation if not complete
`data_issues`	array of {"issue","handling"}	—	problems detected in the design or during the simulation and how handled; empty list if none
`n_input_rows`	object	rows	{"design.csv": number of data rows}
`n_analysed`	integer or null	—	total number of replications simulated (cells × R)
`random_seed`	integer or null	—	the seed of cell C1 (the other seeds are in `design.csv`)

out/figure1.png — coverage against G for the three procedures, one panel per ICC (or an equivalent figure), with the nominal level and the 0.90 failure line marked.
out/report.md — an executable technical note with the sections: Question; Data and exclusions; Method; Results; Uncertainty; Limitations and open issues. "Uncertainty" must discuss the Monte Carlo error of the coverage and bias estimates; "Limitations and open issues" must state what the design does not cover (distributions, heterogeneity, regressors, other weights).
out/analysis.py — running ~/wp29/venv/bin/python out/analysis.py from the workspace root must regenerate out/results.json, out/table1.csv and out/figure1.png from data/design.csv only, deterministically (same seeds → same numbers).
out/claims.json — a list of the central claims of the report: [{"id": "C1", "text": "...", "results_keys": ["..."], "status": "supported"|"unsupported"|"not_evaluated", "caveat": "..."}].

# Constraints

No network access. Python: ~/wp29/venv/bin/python (numpy, pandas, scipy, statsmodels, matplotlib available). Do not modify data/.
Vectorise the bootstrap (a B × G weight matrix per replication is enough); the full grid runs in a few minutes in numpy. The analysis must complete within 10 minutes.
Budget class: transfer mini-project — 12.00 USD cap, 60 minutes wall, 150 turns per session. The project is designed to be completed in well under 45 minutes of work.

spb/wp29_uqo