Clustered standard errors

Clustered standard errors are often useful when treatment is assigned at the level of a cluster instead of at the individual level. For example, suppose that an educational researcher wants to discover whether a new teaching technique improves student test scores. She therefore assigns teachers in "treated" classrooms to try this new technique, while leaving "control" classrooms unaffected. When analyzing her results, she may want to keep the data at the student level (for example, to control for student-level observable characteristics). However, when estimating the standard error or confidence interval of her statistical model, she realizes that classical or even heteroscedasticity-robust standard errors are inappropriate because student test scores within each class are not independently distributed. Instead, students in classes with better teachers have especially high test scores (regardless of whether they receive the experimental treatment) while students in classes with worse teachers have especially low test scores. The researcher can cluster her standard errors at the level of a classroom to account for this aspect of her experiment.^[7]

While this example is very specific, similar issues arise in a wide variety of settings. For example, in many panel data settings (such as difference-in-differences) clustering often offers a simple and effective way to account for non-independence between periods within each unit (sometimes referred to as "autocorrelation in residuals").^[4] Another common and logically distinct justification for clustering arises when a full population cannot be randomly sampled, and so instead clusters are sampled and then units are randomized within cluster. In this case, clustered standard errors account for the uncertainty driven by the fact that the researcher does not observe large parts of the population of interest.^[8]

A useful mathematical illustration comes from the case of one-way clustering in an ordinary least squares (OLS) model. Consider a simple model with N observations that are subdivided in C clusters. Let ${\displaystyle Y}$ be an ${\displaystyle n\times 1}$ vector of outcomes, ${\displaystyle X}$ a ${\displaystyle n\times m}$ matrix of covariates, ${\displaystyle \beta }$ an ${\displaystyle m\times 1}$ vector of unknown parameters, and ${\displaystyle e}$ an ${\displaystyle n\times 1}$ vector of unexplained residuals:

${\displaystyle Y=X\beta +e}$

As is standard with OLS models, we minimize the sum of squared residuals ${\displaystyle e}$ to get an estimate ${\displaystyle {\hat {\beta }}}$:

${\displaystyle \min _{\beta }(Y-X\beta )^{2}}$

${\displaystyle \Rightarrow X'(Y-X{\hat {\beta }})=0}$

${\displaystyle \Rightarrow {\hat {\beta }}=(X'X)^{-1}X'Y}$

From there, we can derive the classic "sandwich" estimator:

${\displaystyle V({\hat {\beta }})=V((X'X)^{-1}X'Y)=V(\beta +(X'X)^{-1}X'e)=V((X'X)^{-1}X'e)=(X'X)^{-1}X'ee'X(X'X)^{-1}}$

Denoting ${\displaystyle \Omega \equiv ee'}$ yields a potentially more familiar form

${\displaystyle V({\hat {\beta }})=(X'X)^{-1}X'\Omega X(X'X)^{-1}}$

While one can develop a plug-in estimator by defining ${\displaystyle {\hat {e}}\equiv Y-X{\hat {\beta }}}$ and letting ${\displaystyle {\hat {\Omega }}\equiv {\hat {e}}{\hat {e}}'}$, this completely flexible estimator will not converge to ${\displaystyle V({\hat {\beta }})}$ as ${\displaystyle N\rightarrow \infty }$. Given the assumptions that a practitioner deems as reasonable, different types of standard errors solve this problem in different ways. For example, classic homoskedastic standard errors assume that ${\displaystyle \Omega }$ is diagonal with identical elements ${\displaystyle \sigma ^{2}}$, which simplifies the expression for ${\displaystyle V({\hat {\beta }})=\sigma ^{2}(X'X)^{-1}}$. Huber-White standard errors assume ${\displaystyle \Omega }$ is diagonal but that the diagonal value varies, while other types of standard errors (e.g. Newey–West, Moulton SEs, Conley spatial SEs) make other restrictions on the form of this matrix to reduce the number of parameters that the practitioner needs to estimate.

Clustered standard errors assume that ${\displaystyle \Omega }$ is block-diagonal according to the clusters in the sample, with unrestricted values in each block but zeros elsewhere. In this case, one can define ${\displaystyle X_{c}}$ and ${\displaystyle \Omega _{c}}$ as the within-block analogues of ${\displaystyle X}$ and ${\displaystyle \Omega }$ and derive the following mathematical fact:

${\displaystyle X'\Omega X=\sum _{c}X'_{c}\Omega _{c}X_{c}}$

By constructing plug-in matrices ${\displaystyle {\hat {\Omega }}_{c}}$, one can form an estimator for ${\displaystyle V({\hat {\beta }})}$ that is consistent as the number of clusters ${\displaystyle c}$ becomes large. While no specific number of clusters is statistically proven to be sufficient, practitioners often cite a number in the range of 30-50 and are comfortable using clustered standard errors when the number of clusters exceeds that threshold.

Alternatively, finite-sample modifications are also typically used, to reduce downwards bias ${\displaystyle V({\hat {\beta }})}$ due to finite C.^[9] Often practitioners use the following bias corrected estimator:

${\displaystyle {\hat {V}}({\hat {\beta }})=(X'X)^{-1}\sum _{c}X'_{c}{\hat {\Omega }}_{c}X_{c}(X'X)^{-1}{\frac {C}{C-1}}{\frac {n-1}{n-k}}.}$

However, more recent practice has shifted towards analogues of the heteroscedasticity-robust HC2 and HC3 estimators.^[9] Often called the CR2 and CR3 estimators, these estimators are unbiased under certain assumptions. They also have been shown, especially when combined with degrees of freedom corrections for use in building confidence intervals, to produce better coverage rates when the number of clusters is not large.^[10][11]