Clustering and Few-Cluster Inference

The Clustering Problem#

You run a regression of student test scores on a school-level policy. You have 500,000 students in 50 states. Your standard errors look tight. Your p-value is 0.001. You feel confident.

You should not.

The policy varies at the state level. Students within the same state share the same policy, the same funding formula, the same regulatory environment. Their residuals are correlated — not because of anything you did wrong, but because the world is structured in groups. Ignoring that correlation makes your standard errors too small, your t-statistics too large, and your confidence intervals too narrow.

Within-cluster correlation is the clustering problem, and getting it wrong is one of the most common inferential errors in applied research.

The Moulton Factor: Why Naive SEs Fail#

The intuition is simple. If observations within clusters are correlated, they contain less independent information than their raw count suggests. Your effective sample size is not the number of observations — it is closer to the number of clusters.

The Moulton factor quantifies exactly how much naive standard errors understate true uncertainty. For a cluster-level regressor, the inflation factor is approximately:

where $m$ is the average cluster size and $\rho$ is the intraclass correlation coefficient — the fraction of total residual variance that is between-cluster rather than within-cluster.

The numbers are startling. Even a modest intraclass correlation of $\rho = 0.05$ with an average cluster size of $m = 50$ yields:

Your naive standard errors are 86% too small. A t-statistic of 2.5 — seemingly significant at the 1% level — becomes 1.34 after correction. Not significant at any conventional level.

Why Does Correlation Arise?#

Within-cluster correlation of residuals arises from several sources:

1. Shared treatments. If a policy is implemented at the state level, all individuals in the state experience the same treatment shock. Any state-level unobservable that affects outcomes creates residual correlation.

2. Shared environments. Students in the same school share teachers, resources, and peer effects. Workers in the same firm share management quality and organizational culture.

3. Sampling design. If you sample clusters (schools, villages) and then sample individuals within clusters, the sampling itself induces correlation.

4. Serial correlation. In panel data, the same unit observed over time has correlated residuals across periods. Serial correlation is clustering over time within units.

When to Cluster: The Abadie et al. Framework#

Not all regressions require clustered standard errors. The decision of whether and how to cluster is not a mechanical one — it depends on the research design, the source of variation, and the sampling process.

Abadie, Athey, Imbens, and Wooldridge (2023) provide a principled framework. Their key insight is that clustering addresses two distinct problems, and you should cluster when either applies:

Reason 1: Clustered Treatment Assignment#

If treatment is assigned at the cluster level (e.g., a policy enacted at the state level), you must cluster at least at that level. The assignment mechanism induces correlation in treatment status within clusters, which in turn creates correlation in residuals.

Rule: Cluster at the level of treatment assignment or higher.

Reason 2: Clustered Sampling#

If you sample clusters from a population of clusters (e.g., you sample 50 schools from 5,000 schools in the country), you should cluster at the sampling level to account for the fact that your clusters are a random draw.

Rule: If clusters are sampled, cluster at the sampling level.

When Not to Cluster#

If treatment varies at the individual level and you observe the full population (or a simple random sample) of clusters, clustering may not be necessary. In this case, heteroscedasticity-robust standard errors suffice.

Cluster-Robust Standard Errors#

Cluster-robust standard errors are computed via a sandwich estimator that allows for arbitrary within-cluster correlation. The key formula for the variance of the OLS estimator $\hat{\beta}$ is:

where $g = 1, \ldots, G$ indexes clusters, $X_g$ is the matrix of regressors for cluster $g$, and $\hat{u}_g$ is the vector of residuals for cluster $g$.

CR0 vs. CR1 vs. CR2#

Different variants apply different finite-sample corrections:

CR0 (no correction): The formula above, used by Stata's default cluster() option (with a degrees-of-freedom adjustment of $G/(G-1) \cdot (N-1)/(N-K)$, making it effectively CR1).

CR1 (degrees-of-freedom correction): Multiplies the variance by $G/(G-1)$, analogous to the HC1 correction for heteroscedasticity-robust standard errors. CR1 is the most commonly used variant in practice.

CR2 (bias-corrected): Applies a cluster-level leverage correction analogous to HC2. The CR2 estimator uses:

where $A_g$ is a correction matrix based on the cluster-level hat matrix. CR2 is less biased than CR1 with few clusters and is particularly recommended when cluster sizes are unbalanced.

Degrees of Freedom#

With clustered standard errors, the relevant degrees of freedom for inference is approximately $G - 1$ (the number of clusters minus one), not $N - K$ (observations minus parameters). This matters enormously when $G$ is small. A t-statistic of 2.1 is significant at 5% with 100 degrees of freedom but not with 10 degrees of freedom (critical value $\approx$ 2.23).

Code: Cluster-Robust Standard Errors#

1library(fixest)
2library(clubSandwich)
3
4# --- CR1 (default in fixest) ---
5model <- feols(outcome ~ treatment + x1 + x2 | state_fe,
6             data = df,
7             cluster = ~state_id)
8summary(model)  # reports CR1 clustered SEs
9
10# --- CR2 (bias-corrected, via clubSandwich) ---
11model_lm <- lm(outcome ~ treatment + x1 + x2 + factor(state_fe),
12             data = df)
13cr2_vcov <- vcovCR(model_lm, cluster = df$state_id, type = "CR2")
14coef_test(model_lm, vcov = cr2_vcov, test = "Satterthwaite")
15
16# --- Two-way clustering (state and year) ---
17model_2way <- feols(outcome ~ treatment + x1 + x2 | state_fe + year_fe,
18                  data = df,
19                  cluster = ~state_id + year_id)
20summary(model_2way)

The Few-Cluster Problem#

Cluster-robust standard errors have an asymptotic justification: they are consistent as the number of clusters $G \to \infty$. But what happens when $G$ is small?

The answer: things break down. And "small" starts sooner than you might think.

What Counts as "Few"?#

There is no sharp cutoff, but the empirical evidence is clear:

Why CR1 Fails with Few Clusters#

The cluster-robust variance estimator involves a sum over $G$ cluster-level "scores" ($X_g' \hat{u}_g$). With few clusters, this sum has high variance — the estimate of the variance is itself highly variable. The result is that:

1. Standard errors are biased downward. With few clusters, CR1 standard errors tend to be too small, leading to over-rejection of the null.

2. The t-distribution approximation is poor. Even with the $G - 1$ degrees-of-freedom correction, the actual sampling distribution of the t-statistic is not well approximated by a t-distribution when $G$ is small.

3. Coverage of confidence intervals is too low. A 95% confidence interval based on CR1 may have actual coverage of 80% or worse with 10 clusters.

Cameron, Gelbach, and Miller (2008) showed through simulation that with $G = 10$ clusters, the actual rejection rate of a 5% test using CR1 standard errors can be 15–25% — three to five times the nominal rate.

Wild Cluster Bootstrap#

The wild cluster bootstrap is the workhorse solution to the few-cluster problem. Proposed by Cameron, Gelbach, and Miller (2008), it provides more reliable inference than conventional cluster-robust standard errors when the number of clusters is small.

How It Works#

The wild cluster bootstrap operates as follows:

1. Estimate the model under the null hypothesis $H_0: \beta_1 = 0$. Obtain restricted residuals $\tilde{u}_i$.

2. Generate bootstrap samples. For each bootstrap iteration $b = 1, \ldots, B$:

   • Draw a random weight $w_g^{(b)} \in \{-1, +1\}$ for each cluster $g$ (Rademacher weights). All observations within cluster $g$ receive the same weight.
   • Construct the bootstrap outcome: $y_i^{(b)} = \hat{y}_i^{R} + w_{g(i)}^{(b)} \cdot \tilde{u}_i$, where $\hat{y}_i^{R}$ is the restricted predicted value.
   • Re-estimate the model on $(y_i^{(b)}, X_i)$ and compute the t-statistic $t^{(b)}$.

3. Compute the bootstrap p-value. The p-value is the fraction of bootstrap t-statistics that exceed the observed t-statistic in absolute value:

Why It Works Better#

The wild cluster bootstrap resamples at the cluster level (the weight is constant within each cluster), preserving the within-cluster correlation structure. By imposing the null hypothesis and computing t-statistics (rather than raw coefficients), it provides an asymptotic refinement — the bootstrap distribution approximates the true finite-sample distribution of the test statistic more accurately than the $t_{G-1}$ approximation.

Webb Weights for Very Few Clusters#

With very few clusters ($G < 10$), even Rademacher weights ($\pm 1$) produce too few distinct bootstrap datasets. Webb (2023) proposed a six-point distribution $\{-\sqrt{3/2}, -\sqrt{2/2}, -\sqrt{1/2}, \sqrt{1/2}, \sqrt{2/2}, \sqrt{3/2}\}$ that provides more variation and better coverage with very small $G$.

Code: Wild Cluster Bootstrap#

1library(fwildclusterboot)
2library(fixest)
3
4# Step 1: Estimate the model with fixest
5model <- feols(outcome ~ treatment + x1 + x2 | state_fe,
6             data = df,
7             cluster = ~state_id)
8
9# Step 2: Wild cluster bootstrap
10boot_result <- boottest(
11model,
12param = "treatment",     # parameter to test
13clustid = ~state_id,     # cluster variable
14B = 9999,                # number of bootstrap iterations
15type = "rademacher"      # weight type ("rademacher", "webb")
16)
17
18# View results
19summary(boot_result)
20# Reports: bootstrap p-value and confidence interval
21
22# For very few clusters, use Webb weights
23boot_webb <- boottest(
24model,
25param = "treatment",
26clustid = ~state_id,
27B = 9999,
28type = "webb"
29)
30summary(boot_webb)

Alternative Solutions for Few Clusters#

The wild cluster bootstrap is the most commonly used approach, but other solutions exist for specific settings:

Randomization Inference#

When the treatment assignment mechanism is known — as in experiments or well-understood natural experiments — randomization inference provides exact p-values that do not depend on the number of clusters. Permute the treatment assignment across clusters and compute the test statistic for each permutation.

Effective Degrees of Freedom (Satterthwaite/BM)#

Rather than using $G - 1$ degrees of freedom, the Bell-McCaffrey (BM) correction estimates the effective degrees of freedom from the data, which can be much smaller than $G - 1$ when clusters are unbalanced. Combined with CR2 standard errors, this approach (CR2 + BM) performs well even with moderate numbers of clusters.

Aggregation to the Cluster Level#

A simple and sometimes underappreciated approach: if your treatment varies only at the cluster level, collapse the data to cluster-level means and run OLS on the collapsed data. With $G$ observations, you use standard (robust) standard errors. This approach is valid when cluster sizes are equal and treatment does not vary within clusters.

where $\bar{Y}_g$ is the cluster-level mean outcome. The drawback is that you lose the ability to control for individual-level covariates (which can improve precision) and that it is less efficient with unequal cluster sizes.

Code: Alternative Approaches#

1# --- CR2 + Satterthwaite degrees of freedom ---
2library(clubSandwich)
3
4model <- lm(outcome ~ treatment + x1 + x2 + factor(state_id),
5          data = df)
6cr2_test <- coef_test(model,
7                    vcov = vcovCR(model,
8                                  cluster = df$state_id,
9                                  type = "CR2"),
10                    test = "Satterthwaite")
11print(cr2_test)
12
13# --- Aggregation to cluster level ---
14df_collapsed <- df |>
15dplyr::group_by(state_id) |>
16dplyr::summarise(
17  outcome_mean = mean(outcome),
18  treatment = first(treatment),
19  x1_mean = mean(x1),
20  x2_mean = mean(x2)
21)
22
23model_collapsed <- lm(outcome_mean ~ treatment + x1_mean + x2_mean,
24                     data = df_collapsed)
25# Use HC2 robust SEs on collapsed data
26library(sandwich)
27library(lmtest)
28coeftest(model_collapsed, vcov = vcovHC(model_collapsed, type = "HC2"))

Two-Way Clustering#

Sometimes residuals are correlated along two dimensions simultaneously. In panel data, residuals may be correlated within firms (across time) and within years (across firms). Two-way clustering accounts for both sources of dependence.

The two-way cluster-robust variance estimator uses the Cameron-Gelbach-Miller (2011) formula:

where $\widehat{V}_{\text{cluster}_1}$ clusters by the first dimension, $\widehat{V}_{\text{cluster}_2}$ clusters by the second, and $\widehat{V}_{\text{intersection}}$ clusters by the intersection of both dimensions.

Failure Demo#

Error Detective#

Concept Checks#

Practical Decision Tree#

When you sit down to estimate standard errors, work through these questions:

1. Does the treatment vary at a higher level than the unit of observation?

• Yes: You must cluster. Go to question 2.
• No: Heteroscedasticity-robust SEs are typically sufficient. Consider clustering if there is a sampling design that induces correlation.

2. At what level should you cluster?

• Cluster at the level of treatment variation, or higher if there is a reason (e.g., sampling design).
• When in doubt, cluster at a broader level. It is conservative but valid.

3. How many clusters do you have?

• $G \geq 50$: CR1 standard errors with $t_{G-1}$ critical values are generally fine.
• $30 \leq G < 50$: Use CR2 + Satterthwaite degrees of freedom. Report wild cluster bootstrap as robustness.
• $G < 30$: Wild cluster bootstrap is strongly recommended. Report both CR1 and bootstrap results.
• $G < 10$: Bootstrap with Webb weights. Consider aggregation to cluster level or randomization inference.

4. Is there correlation along a second dimension?

• If residuals are correlated within both clusters (e.g., firms) and time periods (e.g., years), consider two-way clustering. Remember: the effective cluster count is the minimum of the two dimensions.

Common Mistakes#

How to Report#

A strong reporting practice includes:

Standard errors are clustered at the state level ($G$ = 42), the level at which the policy varies. With 42 clusters, we supplement conventional cluster-robust inference with wild cluster bootstrap p-values (9,999 replications, Rademacher weights) following Cameron, Gelbach, and Miller (2008). The bootstrap p-value for the main treatment effect is 0.041, compared to the conventional clustered p-value of 0.028, confirming that our results are robust to few-cluster concerns.

Key elements to include:

• The level of clustering and why
• The number of clusters ($G$)
• Wild cluster bootstrap results if $G < 50$
• Both conventional and bootstrap p-values for transparency