A Practitioner’s Guide to Cluster-Robust Inference

A. A Simple Example

For simplicity, we begin with OLS with a single regressor that is nonstochastic, and assume no intercept in the model. The results extend to multiple regression with stochastic regressors.

Let y_i = βx_i + u_i, i = 1, …, N, where x_i is nonstochastic and E[u_i] = 0. The OLS estimator can be reexpressed as , so in general

(1)

If errors are uncorrelated over i, then . In the simplest case of homoskedastic errors, V[u_i] = σ² and Equation 1 simplifies to .

If instead errors are heteroskedastic, then Equation 1 becomes

using since E[u_i] = 0. Implementation seemingly requires consistent estimates of each of the N error variances . In a very influential paper, one that extends naturally to the clustered setting, White (1980) noted that instead all that is needed is an estimate of the scalar , and that one can simply use , where is the OLS residual, provided N → ∞. This leads to estimated variance

The resulting standard error for is often called a robust standard error, though a better, more precise, term is heteroskedastic-robust standard error.

What if errors are correlated over i? In the most general case where all errors are correlated with each other,

so

The obvious extension of White (1980) is to use but this equals zero since . Instead, one needs to first set a large fraction of the error correlations E[u_iu_j] to zero. For time series data with errors assumed to be correlated only up to, say, m periods apart as well as heteroskedastic, White’s result can be extended to yield a heteroskedastic-and autocorrelation-consistent (HAC) variance estimate; see Newey and West (1987).

In this paper, we consider clustered errors, with E[u_iu_j] = 0 unless observations i and j are in the same cluster (such as same region). Then

(2)

where the indicator function 1[A] equals 1 if event A happens and equals 0 if event A does not happen. Provided the number of clusters goes to infinity, we can use the variance estimate

(3)

This estimate is called a cluster-robust estimate though, more precisely, it is heteroskedastic-and cluster-robust. This estimate reduces to in the special case that there is only one observation in each cluster.

Typically exceeds due to the addition of terms when i ≠ j. The amount of increase is larger (1) the more positively associated are the regressors across observations in the same cluster (via x_ix_j in Equation 3), (2) the more correlated are the errors (via E[u_iu_j] in Equation 2), and (3) the more observations are in the same cluster (via 1[i, j in same cluster] in Equation 3).

There are several take-away messages. First there can be great loss of efficiency in OLS estimation if errors are correlated within cluster rather than completely uncorrelated. Intuitively, if errors are positively correlated within cluster, then an additional observation in the cluster no longer provides a completely independent piece of new information. Second, failure to control for this within-cluster error correlation can lead to using standard errors that are too small, with consequent overly narrow confidence intervals, overly large t-statistics, and overrejection of true null hypotheses. Third, it is straightforward to obtain cluster-robust standard errors, though they do rely on the assumption that the number of clusters goes to infinity (see Section VI for the few clusters case).

B. Clustered Errors and Two Leading Examples

Let i denote the i^th of N individuals in the sample and g denote the g^th of G clusters. Then for individual i in cluster g the linear model with (one-way) clustering is

(4)

where x_ig is a K × 1 vector. As usual, it is assumed that E[u_ig|x_ig] = 0. The key assumption is that errors are uncorrelated across clusters while errors for individuals belonging to the same cluster may be correlated. Thus

(5)

C. The Cluster-Robust Variance Matrix Estimate

Stacking all observations in the g^th cluster, the Model 4 can be written as

where y_g and u_g are N_g × 1 vectors, X_g is an N_g × K matrix, and there are N_g observations in cluster g. Further stacking y_g, X_g and u_g over the G clusters then yields the model

The OLS estimator is

In general, the variance matrix conditional on X is

(7)

with

(8)

Given error independence across clusters, V[u|X] has a block-diagonal structure, and Equation 8 simplifies to

(9)

The matrix B_clu, the middle part of the “sandwich matrix” Equation 7, corresponds to the numerator of Equation 2. B_clu can be written as:

where ω_{ig, jg} = E[u_igu_jg|Xg] is the error covariance for the ig^th and jg^th observations. We can gain a few insights from inspecting this equation. The term B (and hence ) will be bigger when: (1) regressors within cluster are correlated, (2) errors within cluster are correlated so ω_{ig, jg} is nonzero, (3) N_g is large, and (4) the within-cluster regressor and error correlations are of the same sign (the usual situation). These conditions mirror the more precise Moulton result for the special case of equicorrelated errors given in Equation 6. Both examples in Section II had high within-cluster correlation of the regressor, the DiD example additionally had high within-cluster (serial) correlation of the error and the Moulton (1990) example additionally had N_g large.

Implementation requires an estimate of B_clu given in Equation 9. The cluster-robust estimate of the variance matrix (CRVE) of the OLS estimator is the sandwich estimate

(10)

where

(11)

and is the vector of OLS residuals for the g^th cluster. Formally, Equation 10 and Equation 11 provides a consistent estimate of the variance matrix if as G → ∞. Initial derivations of this estimator, by White (1984, pp. 134–42) for balanced clusters, and by Liang and Zeger (1986) for unbalanced, assumed a finite number of observations per cluster. Hansen (2007a) showed that the CRVE can also be used if N_g → ∞, the case for long panels, in addition to G → ∞. Carter, Schnepel, and Steigerwald (2013) considers unbalanced panels with either N_g fixed or N_g → ∞. The sandwich formula for the CRVE extends to many estimators other than OLS; see Section VII.

Algebraically, the estimator Equation 10 and Equation11 equals Equations 7 and 9 with E[u_gu^′_g] replaced with . What is striking about this is that for each cluster g, the N_g × N_g matrix is bound to be a very poor estimate of the N_g × N_g matrix E[u_gu^′_g]—there is no averaging going on to enable use of a Law of Large Numbers. The “magic” of the CRVE is that despite this, by averaging across all G clusters in Equation 11, we are able to get a consistent variance estimate. This fact helps us to understand one of the limitations of this method in practice—the averaging that makes accurate for is an average based on the number of clusters G. In applications with few clusters, this can lead to problems that we discuss below in Section VI.

Finite-sample modifications of Equation 11 are typically used to reduce downward bias in due to finite G. Stata uses in Equation 11 rather than , with

(12)

In general c ≃ G / (G – 1), though Section IIIB addresses an important exception when fixed effects are directly estimated. Some other packages such as SAS use c = G / (G – 1), a simpler correction that is also used by Stata for extensions to nonlinear models. Either choice of c usually lessens, but does not fully eliminate, the usual downward bias in the CRVE. Other finite-cluster corrections are discussed in Section VI but there is no clear best correction.

D. Feasible GLS

If errors are correlated within cluster, then in general OLS is inefficient and feasible GLS may be more efficient.

Suppose we specify a model for Ω_g = E[u_gu′_g |X_g] in Equation 9, such as within-cluster equicorrelation. Then the GLS estimator is (X′Ω⁻¹X)⁻¹ X′Ω⁻¹y, where Ω = Diag[Ω_g]. Given a consistent estimate of Ω, the feasible GLS estimator of β is

(13)

The FGLS estimator is second-moment efficient, with variance matrix

(14)

under the strong assumption that the error variance Ω is correctly specified.

Remarkably, the cluster-robust method of the previous section can be extended to FGLS. Essentially, OLS is the special case where Ω_g = σ²I_Ng. The cluster-robust estimate of the asymptotic variance matrix of the FGLS estimator is

(15)

where . This estimator requires that u_g and u_h are uncorrelated when g ≠ h, and that G → ∞, but permits E[u_gu′_g | X_g ≠ Ω_g]. The approach of specifying a model for the error variances and then doing inference that guards against misspecification of this model is especially popular in the biostatistics literature that calls Ω_g a “working” variance matrix (see, for example, Liang and Zeger 1986).

There are many possible candidate models for Ω_g, depending on the type of data being analyzed.

For individual-level data clustered by region, the example in Section IIB1, a common starting point model is the random effects (RE) model. The error in Model 4 is specified to have two components:

(16)

where α_g is a cluster-specific error or common shock that is assumed to be independent and identically distributed (iid) , and ε_ig is an idiosyncratic error that is assumed to be iid . Then and . It follows that the intraclass correlation of the error so this model implies equicorrelated errors within cluster. Richer models that introduce heteroskedasticity include random coefficients models and hierarchical linear models.

For panel data, the example in Section IIB2, a range of time series models for u_it may be used, including autoregressive and moving average error models. Analysis of within-cluster residual correlation patterns after OLS estimation can be helpful in selecting a model for Ω_g.

Note that in all cases if cluster-specific fixed effects are included as regressors and N_g is small then bias-corrected FGLS should be used; see Section IIIC.

The efficiency gains of FGLS need not be great. As an extreme example, with equicorrelated errors, balanced clusters, and all regressors invariant within cluster (x_ig = x_g) the FGLS estimator equals the OLS estimator—and so there is no efficiency gain to FGLS. With equicorrelated errors and general X, Scott and Holt (1982) provides an upper bound to the maximum proportionate efficiency loss of OLS, compared to the variance of the FGLS estimator, of 1 / [1 + (4(1 − ρ_u)[1 + (N_max − 1)ρ_u) / (N_max × ρ_u)²], N_max = max{N₁, …, N_G}. This upper bound is increasing in the error correlation ρ_u and the maximum cluster size N_max. For low ρ_u, the maximal efficiency gain can be low. For example, Scott and Holt (1982) notes that for ρ_u = 0.05 and N_max = 20 there is at most a 12 percent efficiency loss of OLS compared to FGLS. With ρ_u = 0.2 and N_max = 100, the efficiency loss could be as much as 86 percent, though this depends on the nature of X.

There is no clear guide to when FGLS may lead to considerable improvement in efficiency, and the efficiency gains can be modest. However, especially in models without cluster-specific fixed effects, implementation of FGLS and use of Equation 15 to guard against misspecification of Ω_g is straightforward. And even modest efficiency gains can be beneficial. It is remarkable that current econometric practice with clustered errors ignores the potential efficiency gains of FGLS.

E. Implementation for OLS and FGLS

For regression software that provides a cluster-robust option, implementation of the CRVE for OLS simply requires defining for each observation a cluster identifier variable that takes one of G distinct values according to the observation’s cluster, and then passing this cluster identifier to the estimation command’s cluster-robust option. For example, if the cluster identifier is id_clu, then Stata OLS command regress y x becomes regress y x, vce(cluster id_clu).

Wald hypothesis tests and confidence intervals are then implemented in the usual way. In some cases, however, joint tests of several hypotheses and of overall statistical significance may not be possible. The CRVE is guaranteed to be positive semi-definite, so the estimated variance of the individual components of are necessarily nonnegative. But is not necessarily positive definite, so it is possible that the variance matrix of linear combinations of the components of is singular. The rank of equals the rank of defined in Equation 11. Because , where is a K × G matrix, it follows that the rank of , and hence that of , is at most the rank of C. Because , the rank of equals min(K, G – 1), though it can be less than this in some examples such as perfect collinearity of regressors and cluster-specific dummy regressors (see Section IIIB for the latter).

A common setting is to have a richly specified model with thousands of observations in far fewer clusters, leading to more regressors than clusters. Then rank-deficient, so it will not be possible to perform an overall F test of the joint statistical significance of all regressors. And in a log-wage regression with occupation dummies and clustering on state, we cannot test the joint statistical significance of the occupation dummies if there are more occupations than states. But it is still okay to perform statistical inference on individual regression coefficients and to do joint tests on a limited number of restrictions (potentially as many as min(K, G – 1)).

Regression software usually also includes a panel data component. Panel commands may enable not only OLS with cluster-robust standard errors but also FGLS for some models of within-cluster error correlation with default (and possibly cluster-robust) standard errors. It is important to note that those panel data commands that do not explicitly use time-series methods—an example is FGLS with equicorrelation of errors within-cluster—can be applied more generally to other forms of clustered data, such as individual-level data with clustering on geographic region.

For example, in Stata first give the command xtset id_clu to let Stata know that the cluster-identifier is variable id_clu. Then the Stata command xtreg y x, pa corr(ind) vce(robust) yields OLS estimates with cluster-robust standard errors. Note that for Stata xt commands, option vce(robust) is generally interpreted as meaning cluster-robust; this is always the case from Version 12.1 on. The xt commands use standard normal critical values whereas command regress uses Student’s T(G – 1) critical values; see Sections VI and VIIA for further discussion.

For FGLS estimation, the commands vary with the model for Ω_g. For equicorrelated errors, a starting point for the example in Section IIB1, the FGLS estimator can be obtained using command xtreg y x, pa corr(exch) or command xtreg y x, re; slightly different estimates are obtained due to slightly different estimates of the equicorrelation. For FGLS estimation of hierarchical models that are richer than a random effects model, use Stata command mixed (Version 13) or xtmixed (earlier versions). For FGLS with panel data and time variable time, first give the command xtset id_clu time to let Stata know both the cluster-identifier and time variable. A starting point for the example in Section IIB2 is an autoregressive error of order 1, estimated using command xtreg y x, pa corr(ar 1). Stata permits a wide range of possible models for serially correlated errors.

In all of these FGLS examples, the reported standard errors are the default ones that assume correct specification of Ω_g. Better practice is to add option vce(robust) for xtreg commands, or option vce(cluster id_clu) for mixed commands, as this yields standard errors that are based on the cluster-robust variance defined in Equation 15.

F. Cluster-Bootstrap Variance Matrix Estimate

Not all econometrics packages compute cluster-robust variance estimates, and even those that do may not do so for all estimators. In that case, one can use a pairs cluster bootstrap that, like the CRVE, gives a consistent estimate of when errors are clustered.

To implement this bootstrap, do the following steps B times: (1) form G clusters by resampling with replacement G times from the original sample of clusters, and (2) compute the estimate of β, denoted in the b^th bootstrap sample. Then, given the B estimates , compute the variance of these

where , and B = 400 should be more than adequate in most settings. It is important that the resampling be done over entire clusters rather than over individual observations. Each bootstrap resample will have exactly G clusters, with some of the original clusters not appearing at all while others of the original clusters may be repeated in the resample two or more times. The terms “pairs” is used as (y_g, X_g) are resampled as a pair. The term “nonparametric” is also used for this bootstrap. Some alternative bootstraps hold X_g fixed while resampling. For finite clusters, if uses in Equation 11 then for comparability should be multiplied by the constant c defined in Equation 12. The pairs cluster bootstrap leads to a cluster-robust variance matrix for OLS with rank K even if K > G.

An alternative resampling method that can be used is the leave-one-cluster-out jackknife. Then, letting denote the estimator of β when the g^th cluster is deleted,

where . This older method can be viewed as an approximation to the bootstrap that does not work as well for nonlinear estimators. It is used less often than the bootstrap and has the same rank as the CRVE.

Unlike a percentile-t cluster bootstrap, presented later, the pairs cluster bootstrap and cluster jackknife variance matrix estimates are no better asymptotically than the CRVE so it is best and quickest to use the CRVE if it is available. But the CRVE is not always available, especially for estimators more complicated than OLS. In that case, one can instead use the pairs cluster bootstrap, though the end of Section VIC indicates potential pitfalls if there are few clusters, or even the cluster jackknife.

In Stata, the pairs cluster bootstrap for OLS without fixed effects can be implemented in several equivalent ways including: regress y x, vce(boot, cluster(id_clu) reps(400) seed(10101)); xtreg y x, pa corr(ind) vce(boot, reps(400) seed(10101)); and bootstrap, cluster(id_clu) reps(400) seed(10101) : regress y x. The last variant can be used for estimation commands and user-written programs that do not have a vce(boot) option. We recommend 400 bootstrap iterations for published results and for replicability one should always set the seed.

For the jackknife the commands are instead, respectively, regress y x, vce(jack, cluster(id_clu)); xtreg y x, pa corr(ind) vce(jack); and jackknife, cluster(id_clu): regress y x. For Stata xt commands, options vce(boot) and vce(jack) are generally interpreted as meaning cluster bootstrap and cluster jackknife; always so from Stata 12.1 on.

A. Do Fixed Effects Fully Control for Within-Cluster Correlation?

While cluster-specific effects will control for part of the within-cluster correlation of the error, in general they will not completely control for within-cluster error correlation (not to mention heteroskedasticity). So the CRVE should still be used. There are several ways to make this important point.

Suppose we have data on students in classrooms in schools. A natural model, a special case of a hierarchical model, is to suppose that there is both an unobserved school effect and, on top of that, an unobserved classroom effect. Letting i denote individual, s school, and c classroom, we have . A regression with school-level fixed effects (or random effects) will control for within-school correlation but not the additional within-classroom correlation.

Suppose we have a short panel (T fixed, N → ∞) of uncorrelated individuals and estimate . Then the error u_it may be correlated over time (that is, within-cluster) due to omitted factors that evolve progressively over time. Inoue and Solon (2006) provides a test for this serial correlation. Cameron and Trivedi (2005, p. 710) presents an FE individual-level panel data log-earnings regressed on log-hours example with cluster-robust standard errors four times the default. Serial correlation in the error may be due to omitting lagged y as a regressor. When y_i,t–1 is included as an additional regressor in the FE model, the Arellano-Bond estimator is used and, even with y_i,t–1 included, the Arellano-Bond methodology requires that we first test whether the remaining error u_it is serially correlated.

Finally, suppose we have a single cross-section (or a single time series). This can be viewed as regression on a single cluster. Then in the model (or ), the intercept is the cluster-specific fixed effect. There are many reasons for why the error u_i (or u_t) may be correlated in this regression.

B. Cluster-Robust Variance Matrix with Fixed Effects

Since inclusion of cluster-specific fixed effects may not fully control for cluster correlation (and/or heteroskedasticity), default standard errors that assume u_ig to be iid may be invalid. So one should use cluster-robust standard errors.

Arellano (1987) showed that defined in Equations 10 and 11 remains valid for the within estimator that controls for inclusion of G cluster-specific fixed effects, provided G → ∞ and N_g is small. If instead one obtains the LSDV estimator, the CRVE formula gives the same CRVE for as that for the within estimator, with the important proviso that the same degrees-of-freedom adjustment must be used—see below. The fixed effects estimates obtained for the LSDV estimator are essentially based only on N_g observations, so is inconsistent for , just as is inconsistent for α_g.

Hansen (2007a, p. 600) shows that this CRVE can also be used if additionally N_g → ∞, for both the case where within-cluster correlation is always present (for example, for many individuals in each village) and for the case where within-cluster correlation eventually disappears (for example, for panel data where time series correlation disappears for observations far apart in time). The rates of convergence are in the first case and in the second case but the same asymptotic variance matrix is obtained in either case. Kézdi (2004) analyzed the CRVE in FE models for a range of values of G and N_g.

It is important to note that, while LSDV and within estimation lead to identical estimates of β, they can yield different standard errors due to different finite sample degrees-of-freedom correction.

It is well known that if default standard errors are used, that is, it is assumed that u_ig in Equation 17 is iid, then one can safely use standard errors after LSDV estimation as this correctly views the number of parameters as G + K rather than K. If instead the within estimator is used, however, manual OLS estimation of Equation 18 will mistakenly view the number of parameters to equal K rather than G + K. (Built-in panel estimation commands for the within estimator—that is, a fixed effects command—should remain okay to use because they should be programmed to use G + K in calculating the standard errors.)

It is not well known that if cluster-robust standard errors are used, and cluster sizes are small, then inference should be based on the within estimator standard errors. We thank Arindrajit Dube and Jason Lindo for bringing this issue to our attention. Within and LSDV estimation lead to the same cluster-robust standard errors if we apply Formula 11 to the respective regressions, or if we multiply this estimate by c = G / (G – 1). Differences arise, however, if we multiply by the small-sample correction c given in Equation12. Within estimation sets c = [G / (G −1)] / [(N −1) / (N − K + 1)] since there are only (K – 1) regressors—the within model is estimated without an intercept. LSDV estimation uses c = [G / (G −1)] / [(N −1) / (N − G – K + 1)] since the G cluster dummies are also included as regressors. For balanced clusters with N_g = N_* and G large relative to K, c ≃ 1 for within-estimation and c ≃ N_* / (N_* – 1) for LSDV estimation. Suppose there are only two observations per cluster, due to only two individuals per household or two time periods in a panel setting, so N_g = N_* = 2. Then c ≃ 2 / (2 – 1) = 2 for LSDV estimation, leading to CRVE that is twice that from within estimation. Within estimation leads to the correct finite-sample correction.

In Stata, the within command xtreg y x, fe vce(robust) gives the desired CRVE. The alternative LSDV commands regress y x i.id_clu, vce(cluster id_clu) and, equivalently, regress y x, absorb(id_clu) vce(cluster id_clu) use the wrong degrees-of-freedom correction. If a CRVE is needed, then use xtreg. If there is reason to instead use regress i.id, then the cluster-robust standard errors should be multiplied by the square root of [N − (K − 1)] / [N − G − (K − 1)], especially if N_g is small.

The inclusion of cluster-specific dummy variables increases the dimension of the CRVE but does not lead to a corresponding increase in its rank. To see this, stack the dummy variable dh_ig for cluster g into the N_g × 1 vector dh_g. Then , by the OLS normal equations, leading to the rank of falling by one for each cluster-specific effect. If there are k regressors varying within cluster and G – 1 dummies then, even though there are K + G – 1 parameters β, the rank of is only the minimum of K and G – 1. And a test that α₁, …, α_G are jointly statistically significant is a test of G – 1 restrictions (since the intercept or one of the fixed effects needs to be dropped). So even if the cluster-specific fixed effects are consistently estimated (that is, if N_g → ∞), it is not possible to perform this test if K < G – 1, which is often the case.

If cluster-specific effects are present, then the pairs cluster bootstrap must be adapted to account for the following complication. Suppose Cluster 3 appears twice in a bootstrap resample. Then if clusters in the bootstrap resample are identified from the original cluster-identifier, the two occurrences of Cluster 3 will be incorrectly treated as one large cluster rather than two distinct clusters.

In Stata, the bootstrap option idcluster ensures that distinct identifiers are used in each bootstrap resample. Examples are regress y x i.id_clu, vce(boot, cluster(id_clu) idcluster(newid) reps(400) seed(10101)) and, more simply, xtreg y x, fe vce(boot, reps(400) seed(10101)), as in this latter case Stata automatically accounts for this complication.

C. Feasible GLS with Fixed Effects

When cluster-specific fixed effects are present, more efficient FGLS estimation can become more complicated. In particular, if asymptotic theory relies on G → ∞ with N_g fixed, the α_g cannot be consistently estimated. The within estimator of β is nonetheless consistent, as α_g disappears in the mean-differenced model. But the resulting residuals are contaminated because they depend on both and , and these residuals will be used to form a FGLS estimator. This leads to bias in the FGLS estimator so one needs to use bias-corrected FGLS unless N_g → ∞. The correction method varies with the model for Ω_g = V[u_g], and currently there are no Stata user-written commands to implement these methods.

For panel data, a commonly used model specifies an AR(p) model for the errors u_ig in Equation 17. If fixed effects are present, then there is a bias (of order N_g^-1) in estimation of the AR(p) coefficients. Hansen (2007b) obtains bias-corrected estimates of the AR(p) coefficients and uses these in FGLS estimation. Hansen (2007b) in simulations shows considerable efficiency gains in bias-corrected FGLS compared to OLS.

Brewer, Crossley, and Joyce (2013) considers a DiD model with individual-level U.S. panel data with N = 750,127, T = 30, and a placebo state-level law, so clustering is on state with G = 50. It finds that bias-corrected FGLS for AR(2) errors, using the Hansen (2007b) correction, leads to higher power than FE estimation. In its example, ignoring the bias correction does not change results much, perhaps because T = 30 is reasonably large.

For balanced clusters with Ω_g the same for all g, say Ω_g = Ω_*, and for N_g small, then the FGLS estimator in Equation 13 can be used without need to specify a model for Ω_*. Instead we can let have ij^th entry , where are the residuals from initial OLS estimation. These assumptions may be reasonable for a balanced panel. Two complications can arise. First, even without fixed effects there may be many off-diagonal elements to estimate and this number can be large relative to the number of observations. Second, the fixed effects lead to bias in estimating the off-diagonal covariances. Hausman and Kuersteiner (2008) presents fixes for both complications.

D. Testing the Need for Fixed Effects

FE estimation is often accompanied by a loss of precision in estimation, as only within-cluster variation is used (recall we regress on ). Furthermore, the coefficient of a cluster-invariant regressor is not identified because then . Thus, it is standard to test whether it is sufficient to estimate by OLS or FGLS, without cluster-specific fixed effects.

The usual test is a Hausman test based on the difference between the FE estimator, , and the RE estimator, . Let β₁ denote a subcomponent of β, possibly just the coefficient of a single regressor of key interest; at most, β₁ contains the coefficients of all regressors that are not invariant within cluster or, in the case of panel data, that are not aggregate time effects that take the same value for each individual. The chi-squared distributed test statistic is

where is a consistent estimate of .

Many studies use the standard form of the Hausman test. This obtains under the strong assumption that is fully efficient under the null hypothesis. This requires the unreasonably strong assumptions that α_g and ε_ig in Equation 16 are iid, requiring that neither α_g and ε_ig is heteroskedastic and that ε_ig has no within-cluster correlation. As already noted, these assumptions are likely to fail and one should not use default standard errors. Instead, a CRVE should be used. For similar reasons, the standard form of the Hausman test should not be used.

Wooldridge (2010, p. 332) instead proposes implementing a cluster-robust version of the Hausman test by the following OLS regression

where w_g denotes the subcomponent of x_ig that varies within cluster and . If H₀ : = γ = 0 is rejected using a Wald test based on a cluster-robust estimate of the variance matrix, then the fixed effects model is necessary. The Stata user-written command xtoverid, due to Schaffer and Stillman (2010), implements this test.

An alternative is to use the pairs cluster bootstrap to obtain , in each resample obtaining and , leading to B resample estimates of . We are unaware of studies comparing these two cluster-robust versions of the Hausman test.

A. Factors Determining What to Cluster Over

There are two guiding principles that determine what to cluster over.

First, given defined in Equations 7 and 9, whenever there is reason to believe that both the regressors and the errors might be correlated within cluster, we should think about clustering defined in a broad enough way to account for that clustering. Going the other way, if we think that either the regressors or the errors are likely to be uncorrelated within a potential group, then there is no need to cluster within that group.

Second, is an average of G terms that gets closer to only as G gets large. If we define very large clusters, so that there are very few clusters to average over in Equation 11, then the resulting can be a very poor estimate of . This complication, and discussion of how few is “few,” is the subject of Section VI.

These two principles mirror the bias-variance tradeoff that is common in many estimation problems—larger and fewer clusters have less bias but more variability. There is no general solution to this tradeoff, and there is no formal test of the level at which to cluster. The consensus is to be conservative and avoid bias and to use bigger and more aggregate clusters when possible, up to and including the point at which there is concern about having too few clusters.

For example, suppose your data set included individuals within counties within states, and you were considering whether to cluster at the county level or the state level. We have been inclined to recommend clustering at the state level. If there was within-state cross-county correlation of the regressors and errors, then ignoring this correlation (for example, by clustering at the county level) would lead to incorrect inference. In practice, researchers often cluster at progressively higher (that is, broader) levels and stop clustering when there is relatively little change in the standard errors. This seems to be a reasonable approach.

There are settings where one may not need to use cluster-robust standard errors. We outline several though note that in all these cases it is always possible to still obtain cluster-robust standard errors and contrast them to default standard errors. If there is an appreciable difference, then use cluster-robust standard errors.

If a key regressor is randomly assigned within clusters, or is as good as randomly assigned, then the within-cluster correlation of the regressor is likely to be zero. Thus, there is no need to cluster standard errors even if the model’s errors are clustered. In this setting, if there are additionally control variables of interest, and if these are not randomly assigned within cluster, then we may wish to cluster our standard errors for the sake of correct inference on the control variables.

If the model includes cluster-specific fixed effects, and we believe that within-cluster correlation of errors is solely driven by a common shock process, then we may not be worried about clustering. The fixed effects will absorb away the common shock, and the remaining errors will have zero within-cluster correlation. More generally, control variables may absorb systematic within-cluster correlation. For example, in a state-year panel setting, control variables may capture the state-specific business cycle.

However, as already noted in Section IIIA, the within-cluster correlation is usually not fully eliminated. And even if it is eliminated, the errors may still be heteroskedastic. Stock and Watson (2008) shows that, surprisingly, applying the usual White (1980) heteroskedastic-consistent variance matrix estimate to the FE estimator leads to inconsistent estimation of if N_g is small (though it is correct if N_g = 2). It derives a bias-corrected formula for heteroskedastic-robust standard errors. Alternatively, and more simply, the CRVE is consistent for even if the errors are only heteroskedastic, though this estimator of is more variable.

Finally, as already noted in Section IID, we can always build a parametric model of the correlation structure of the errors and estimate by FGLS. If we believe that this parametric model captures the salient features of the error correlations, then default FGLS standard errors can be used.

B. Clustering Due to Survey Design

Clustering routinely arises due to the sampling methods used in complex surveys. Rather than randomly draw individuals from the entire population, costs are reduced by sampling only a randomly selected subset of primary sampling units (such as a geographic area), followed by random selection, or stratified selection, of people within the chosen primary sampling units.

The survey methods literature uses methods to control for clustering that predate the cluster-robust approach of this paper. The loss of estimator precision due to clustered sampling is called the design effect: “The design effect or Deff is the ratio of the actual variance of a sample to the variance of a simple random sample of the same number of elements” (Kish 1965, p. 258). Kish and Frankel (1974) gives the variance inflation formula Equation 6 in the nonregression case of estimation of the mean. Pfeffermann and Nathan (1981) considers the more general regression case. The CRVE is called the linearization formula, and the common use of G – 1 as the degrees of freedom used in hypothesis testing comes from the survey methods literature; see Shah, Holt, and Folsom (1977), which predates the econometrics literature.

Applied economists routinely use data from complex surveys, controlling for clustering by using a cluster-robust variance matrix estimate. At the minimum, one should cluster at the level of the primary sampling unit though often there is reason to cluster at a broader level, such as clustering on state if regressors and errors are correlated within state.

The survey methods literature additionally controls for two other features of survey data—weighting and stratification. These methods are well-established and are incorporated in specialized software as well as in some broad-based packages such as Stata. Bhattacharya (2005) provides a comprehensive treatment in a GMM framework.

If sample weights are provided, then it is common to perform weighted least squares. Then the CRVE for is that given in Equation 15 with replaced by W_g. The need to weight can be ignored if stratification is on only the exogenous regressors and we assume correct specification of the model so that in our sample E[y|X] = X. In that special case, both weighted and unweighted estimators are consistent, and weighted OLS actually may be less efficient if, for example, model errors are iid; see, for example, Solon, Haider, and Wooldridge (2013). Another situation in which to use weighted least squares, unrelated to complex surveys, is when data for the ig^th observation is obtained by in turn averaging N_ig observations and N_ig varies.

Stratification of the sample can enable more precise statistical inference. These gains can be beneficial in the nonregression case, such as estimating the monthly national unemployment rate. The gains can become much smaller once regressors are included because these can partially control for stratification; see, for example, the application in Bhattacharya (2005). Econometrics applications therefore usually do not adjust standard errors for stratification, leading to conservative inference due to some relatively small overestimation of the standard errors.

A. Multiway Cluster-Robust Variance Matrix Estimate

The cluster-robust estimate of defined in Equations 10 and 11 can be generalized to clustering in multiple dimensions. In a change of notation, suppress the subscript for cluster and more simply denote the model for an individual observation as

(19)

Regular one-way clustering is based on the assumption that E[u_iu_j|x_i,x_j] = 0, unless observations i and j are in the same cluster. Then Equation 11 sets [i, j in same cluster], where . In multiway clustering, the key assumption is that E[u_iu_j|x_i,x_j] = 0, unless observations i and j share any cluster dimension. Then the multiway cluster robust estimate of replaces Equation 11 with

(20)

This method relies on asymptotics that are in the number of clusters of the dimension with the fewest number of clusters. This method is thus most appropriate when each dimension has many clusters.

Theory for two-way cluster robust estimates of the variance matrix is presented in Cameron, Gelbach, and Miller (2006, 2011), Miglioretti and Heagerty (2006), and Thompson (2006, 2011). See also empirical panel data applications by Acemoglu and Pischke (2003), which clustered at individual level and at region×time level, and by Petersen (2009), which clustered at firm level and at year level. Cameron, Gelbach, and Miller (2011) presents an extension to multiway clustering. Like one-way cluster-robust, the method can be applied to estimators other than OLS.

For two-way clustering, this robust variance estimator is easy to implement given software that computes the usual one-way cluster-robust estimate. First, obtain three different cluster-robust “variance” matrices for the estimator by one-way clustering in, respectively, the first dimension, the second dimension, and by the intersection of the first and second dimensions. Then add the first two variance matrices and, to account for double-counting, subtract the third. Thus

(21)

where the three component variance estimates are computed using Equations 10 and 11 for the three different ways of clustering.

We spell this out in a step-by-step fashion.

• (1) Identify your two ways of clustering. Make sure you have a variable that identifies each way of clustering. Also create a variable that identifies unique “group 1 by group 2” combinations. For example, suppose you have individual-level data spanning many U.S. states and many years, and you want to cluster on state and on year. You will need a variable for state (for example, California), a variable for year (for example, 1990), and a variable for state-by-year (California and 1990).

• (2) Estimate your model, clustering on “group 1.” For example, regress y on x, clustering on state. Save the variance matrix as .

• (3) Estimate your model, clustering on “group 2.” For example, regress y on x, clustering on year. Save the variance matrix as .

• (4) Estimate your model, clustering on “group 1 by group 2.” For example, regress y on x, clustering on state-by-year. Save the variance matrix as .

• (5) Create a new variance matrix . This is your new two-way cluster robust variance matrix for .

• (6) Standard errors are the square root of the diagonal elements of this matrix.

If you are interested in only one coefficient, you can also just focus on saving the standard error for this coefficient in Steps 2–4 above, and then create .

In taking these steps, you should watch out for some potential pitfalls. With perfectly multicollinear regressors such as inclusion of dummy variables—some of which are redundant—a statistical package may automatically drop one or more variables to ensure a nonsingular set of regressors. If the package happens to drop different sets of variables in Steps 2, 3, and 4, then the resulting s will not be comparable, and adding them together in Step 5 will give a nonsense result. To prevent this issue, manually inspect the estimation results in Steps 2, 3, and 4 to ensure that each step has the same set of regressors, the same number of observations, etc. The only things that should be different are the reported standard errors and the reported number of clusters.

B. Implementation

Unlike the standard one-way cluster case, is not guaranteed to be positive semidefinite, so it is possible that it may compute negative variances. In some applications with fixed effects, may be non positive-definite, but the subcomponent of associated with the regressors of interest may be positive-definite. This may lead to an error message, even though inference is appropriate for the parameters of interest. Our informal observation is that this issue is most likely to arise when clustering is done over the same groups as the fixed effects. Few clusters in one or more dimensions can also lead to being a nonpositive-semidefinite matrix. Cameron, Gelbach, and Miller (2011) presents an eigendecomposition technique used in the time series HAC literature that zeroes out negative eigenvalues in to produce a positive semidefinite variance matrix estimate.

The Stata user-written command cmgreg, available at the authors’ websites, implements multiway clustering for the OLS estimator with, if needed, the negative eigen-value adjustment. The Stata add-on command ivreg2, due to Baum, Schaffer, and Still-man (2007), implements two-way clustering for OLS, IV and linear GMM estimation. Other researchers have also posted code, available from searching the Web.

Cameron, Gelbach, and Miller (2011) applies the two-way method to data from Hersch (1998) that examines the relationship between individual wages and injury risk measured separately at the industry level and the occupation level. The log-wage for 5,960 individuals is regressed on these two injury risk measures, with standard errors obtained by two-way clustering on 211 industries and 387 occupations. In that case, two-way clustering leads to only a modest change in the standard error of the industry job risk coefficient compared to the standard error with one-way clustering on industry. Because industry job risk is perfectly correlated within industry, by Result 6 we need to cluster on industry if there is any within-industry error correlation. By similar logic, the additional need to cluster on occupation depends on the within-occupation correlation of job industry risk, and this correlation need not be high. For the occupation job risk coefficient, the two-way and one-way cluster (on occupation) standard errors are similar. Despite the modest difference in this example, two-way clustering avoids the need to report standard errors for one coefficient clustering in one way and for the second coefficient clustering in the second way.

Cameron, Gelbach, and Miller (2011) also applies the two-way cluster-robust method to data on volume of trade between 98 countries with 3,262 unique country pairs. In that case, two-way clustering on each of the countries in the country pair leads to standard errors that are 36 percent larger than one-way clustering and 230 percent more than heteroskedastic-robust standard errors. Cameron and Miller (2012) studies such dyadic data in further detail. It notes that two-way clustering does not pick up all the potential correlations in the data. Instead, more general cluster-robust methods, including one proposed by Fafchamps and Gubert (2007), should be used.

C. Feasible GLS

Similar to one-way clustering, FGLS is more efficient than OLS, provided an appropriate model for Ω = E[uu′ |X] is specified and is consistently estimated.

The random effects model can be extended to multiway clustering. For individual i in clusters g and h, the two-way random effects model specifies

where the errors α_g, δ_h, and ε_igh are each assumed to be iid distributed with mean 0. For example, Moulton (1986) considered clustering due to grouping of regressors (schooling, age, and weeks worked) in a log earnings regression, and estimated a model with common random shock for each year of schooling, for each year of age, and for each number of weeks worked.

The two-way random effects model can be estimated using standard software for (nested) hierarchical linear models. See, for example, Cameron and Trivedi (2009, Chapter 9.5.7) for Stata command xtmixed (command mixed from Version 13 on). For estimation of a many-way random effects model, see Davis (2002), which modeled film attendance data clustered by film, theater, and time.

The default standard errors after FGLS estimation require that Ω is correctly specified. For two-way and multiway random effects models, FGLS estimation entails transforming the data in such a way that there is no obvious method for computing a variance matrix estimate that is robust to misspecification of Ω. Instead, if there is concern about misspecification of Ω, then one needs to consider FGLS with richer models for Ω and assume that these are correctly specified—see Rabe-Hesketh and Skrondal (2012) for richer hierarchical models in Stata—or do less efficient OLS with two-way cluster-robust standard errors.

D. Spatial Correlation

Cluster-robust variance matrix estimates are closely related to spatial-robust variance matrix estimates.

In general, for Model 19, in Equation 20 has the form

(22)

where w(i, j) are weights. For cluster-robust inference, these weights are either one (cluster in common) or zero (no cluster in common). But the weights can also decay from one to zero, as in the case of the HAC variance matrix estimate for time series where w(i, j) decays to zero as |i – j| increases.

For spatial data, it is assumed that model errors become less correlated as the spatial distance between observations grows. For example, with state-level data the assumption that model errors are uncorrelated across states may be relaxed to allow correlation that decays to zero as the distance between states gets large. Conley (1999) provides conditions under which Equations 10 and 22 provide a robust variance matrix estimate for the OLS estimator, where the weights w(i, j) decay with the spatial distance. The estimate (which Conley also generalizes to GMM models) is often called a spatial-HAC estimate rather than spatial-robust, as proofs use mixing conditions (to ensure decay of dependence) as observations grow apart in distance. These conditions are not applicable to clustering due to common shocks, which leads to the cluster-robust estimator with independence of observations across clusters.

Driscoll and Kraay (1998) considers panel data with T time periods and N individuals, with errors potentially correlated across individuals (and no spatial dampening), though this correlation across individuals disappears for observations that are more than m time periods apart. Let it denote the typical observation. The Driscoll-Kraay spatial correlation consistent (SCC) variance matrix estimate can be shown to use weight w(it, js) = 1 – d(it, js) / (m + 1) in Equation 22, where the summation is now over i, j, s, and t, and d(it, js) = | t – s| if |t – s| ≥ m and d(it, js) = 0 otherwise. This method requires that the number of time periods T → ∞, so is not suitable for short panels while N may be fixed or N → ∞. The Stata add-on command xtscc, due to Hoechle (2007), implements this variance estimator.

An estimator proposed by Thompson (2006) allows for across-cluster (in its example firm) correlation for observations close in time in addition to within-cluster correlation at any time separation. The Thompson estimator can be thought of as using w(it,js) = 1[i, j share a firm, or d(it, js) ≤ m]. Foote (2007) contrasts the two-way cluster-robust and these other variance matrix estimators in the context of a macroeconomics example. Petersen (2009) contrasts various methods for panel data on financial firms where there is concern about both within firm correlation (over time) and across firm correlation due to common shocks.

Barrios et al. (2012) considers state-year panel data on individuals in states over years with state-level treatment and outcome (earnings) that is correlated spatially across states. This spatial correlation can be ignored if the state-level treatment is randomly assigned. But if the treatment is correlated over states (for example, adjacent states may be more likely to have similar minimum wage laws) then one can no longer use standard errors clustered at the state level. Instead, one should additionally allow for spatial correlation of errors across states. The authors additionally contrast traditional model-based inference with randomization inference.

In practice, data can have cluster, spatial, and time series aspects, leading to hybrids of cluster-robust, spatial-HAC, and time-series HAC estimators. Furthermore, it may be possible to parameterize some of the error correlation. For example, for a time series AR(1) error it may be preferable to use based on an AR(1) model rather than . To date, empirical practice has not commonly modeled these combined types of error correlations. This may become more common in the future.

A. The Basic Problems with Few Clusters

There are two main problems with few clusters. First, OLS leads to “overfitting,” with estimated residuals systematically too close to zero compared to the true error terms. This leads to a downward-biased cluster-robust variance matrix estimate. The second problem is that even with bias-correction the use of fitted residuals to form the estimate of b leads to overrejection (and confidence intervals that are too narrow) if the critical values are from the standard normal or even the T(G – 1) distribution.

For the linear regression model with independent homoskedastic normally distributed errors, similar problems are easily controlled. An unbiased variance matrix is obtained by estimating the error variance σ² by rather than . This is the “fix” in the OLS setting for the first problem. The analogue to the second problem is that the N[0,1] distribution is a poor approximation to the true distribution of the Wald statistic. In the iid case, the Wald statistic w can be shown to be exactly T(N – K) distributed. For nonnormal homoskedastic errors, the T(N – K) is still felt to provide a good approximation provided N is not too small. Both of these problems arise in the clustered setting albeit with more complicated manifestations and fixes.

For independent heteroskedastic normally distributed errors there are no exact results. MacKinnon and White (1985) considers several adjustments to the heteroskedastic-consistent variance estimate of White (1980) including one called HC2 that is unbiased in the special case that errors are homoskedastic. Unfortunately, if errors are actually heteroskedastic, as expected, the HC2 estimate is no longer unbiased for —an unbiased estimator depends on the unknown pattern of heteroskedasticity and on the design matrix X. And there is no way to obtain an exact T distribution result for w even if errors are normally distributed. Other proposed solutions for testing and forming confidence intervals include using a T distribution with data-determined degrees of freedom and using a bootstrap with asymptotic refinement.

In the following subsections, we consider extensions of these various adjustments to the clustered case where the problems can become even more pronounced.

Before proceeding, we note that there is no specific point at which we need to worry about few clusters. Instead, “more is better.” Current consensus appears to be that G = 50 is enough for state-year panel data. In particular, Bertrand, Duflo, and Mullainathan (2004, Table 8) finds in its simulations that for a policy dummy variable with high within-cluster correlation, a Wald test based on the standard CRVE with critical value of 1.96 had rejection rates of 0.063, 0.058, 0.080, and 0.115 for number of states (G) equal to, respectively, 50, 20, 10, and 6. The simulations of Cameron, Gelbach, and Miller (2008, Table 3), based on a quite different data generating process but again with standard CRVE and critical value of 1.96, had rejection rates of 0.068, 0.081, 0.118, and 0.208 for G equal to, respectively, 30, 20, 10, and 5. In both cases, the rejection rates would also exceed 0.05 if the critical value was from the T(G – 1) distribution.

The preceding results are for balanced clusters. Cameron, Gelbach, and Miller (2008, Table 4, Column 8) considers unbalanced clusters when G = 10. The rejection rate with unbalanced clusters, half of size N_g = 10 and half of size 50, is 0.183, appreciably worse than rejection rates of 0.126 and 0.115 for balanced clusters of sizes, respectively, 10 and 100. Recent papers by Carter, Schnepel, and Steigerwald (2013) and Imbens and Kolesar (2012) provide theory that also indicates that the effective number of clusters is reduced when N_g varies across clusters; see also the simulations in MacKinnon and Webb (2013). Similar issues may also arise if the clusters are balanced, but estimation is by weighted LS that places different weights on different clusters. Cheng and Hoekstra (2013) documents that weighting can result in overrejection in the panel DiD setting of Bertrand, Duflo, and Mullainathan (2004).

To repeat a key message, there is no clear-cut definition of “few.” Depending on the situation, “few” may range from less than 20 clusters to less than 50 clusters in the balanced case and even more clusters in the unbalanced case.

B. Solution 1: Bias-Corrected Cluster-Robust Variance Matrix

A weakness of the standard CRVE with residual is that it is biased for , since . The bias depends on the form of Ω_g but will usually be downward. Several corrected residuals to use in place of in Equation 11 have been proposed. The simplest, already mentioned, is to use or where c is defined in Equation 12. One should at least use either of these corrections.

Bell and McCaffrey (2002) uses

(24)

where H_gg = X_g(X′X)^-1 X′g. This transformed residual leads to unbiased CRVE in the special case that E[u_gu′_g] = σ²I. This is a cluster generalization of the HC2 variance estimate of MacKinnon and White (1985), so we refer to it as the CR2VE.

Bell and McCaffrey (2002) also uses

(25)

This transformed residual leads to CRVE that can be shown to equal the delete-one-cluster jackknife estimate of the variance of the OLS estimator. This jackknife correction leads to upward-biased CRVE if in fact E[u_gu′_g] = σ²I. This is a cluster generalization of the HC3 variance estimate of MacKinnon and White (1985), so we refer to it as the CR3VE.

Angrist and Pischke (2009, Chapter 8) and Cameron, Gelbach, and Miller (2008) provide a more extensive discussion and cite more of the relevant literature. This literature includes papers that propose corrections for the more general case that E[u_gu′_g] ≠ σ²I but has a known parameterization, such as an RE model, and extension to generalized linear models.

Angrist and Lavy (2002) applies the CR2VE Correction 24 in an application with G = 30 to 40 and finds that the correction increases cluster-robust standard errors by between 10 and 50 percent. Cameron, Gelbach, and Miller (2008, Table 3) finds that the CR3VE Correction 24 has rejection rates of 0.062, 0.070, 0.092, and 0.138 for G equal to, respectively, 30, 20, 10, and 5. These rejection rates are a considerable improvement on 0.068, 0.081, 0.118, and 0.208 for the standard CRVE, but there is still considerable overrejection for very small G.

The literature has gravitated to using the CR2VE adjustment rather than the CR3VE adjustment. This reduces but does not eliminate overrejection when there are few clusters.

C. Solution 2: Cluster Bootstrap with Asymptotic Refinement

In Section IIF, we introduced the bootstrap as it is usually used, to calculate standard errors that rely on regular asymptotic theory. Here, we consider a different use of the bootstrap, one with asymptotic refinement that may lead to improved finite-sample inference.

We consider inference based on G → ∞ (formally, has a limit normal distribution). Then a two-sided Wald test of nominal size 0.05, for example, can be shown to have true size 0.05 + O(G^–1) when the usual asymptotic normal approximation is used. For G → ∞ this equals the desired 0.05 but for finite G this differs from 0.05. If an appropriate bootstrap with asymptotic refinement is instead used, the true size is 0.05 + O(G^–3/2). This is closer to the desired 0.05 for large G, as G^–3/2 < G^–1. Hopefully, this is also the case for small G, something that is established using appropriate Monte Carlo experiments. For a one-sided test or a nonsymmetric two-sided test the rates are instead, respectively, 0.05 + O(G^–1/2) and 0.05 + O(G^–1).

Asymptotic refinement can be achieved by bootstrapping a statistic that is asymptotically pivotal, meaning the asymptotic distribution does not depend on any unknown parameters. The estimator is not asymptotically pivotal as its distribution depends on , which in turn depends on unknown variance parameters in V[u|X]. The Wald t-statistic defined in Equation 23 is asymptotically pivotal as its asymptotic distribution is N[0,1], which does not depend on unknown parameters.

D. Solution 3: Improved Critical Values Using a T-Distribution

The simplest small-sample correction for the Wald statistic is to base inference on a T distribution rather than the standard normal, with degrees of freedom at most the number of clusters G. Recent research has proposed methods that lead to using degrees of freedom much less than G, especially if clusters are unbalanced.

E. Special Cases

With few clusters, additional results can be obtained if there are many observations in each group. In DiD studies, the few-clusters problem arises if few groups are treated, even if G is large. And the few-clusters problem is more likely to arise if there is multiway clustering.

A. Cluster-Robust F-tests

Consider Wald joint tests of several restrictions on the regression parameters. Except in the special case of linear restrictions and OLS with iid normal errors, asymptotic theory yields only a chi-squared distributed statistic, such as W, that is χ²(h) distributed where h is the number of (linearly independent) restrictions.

Alternatively, we can use the related F statistic, F = W / h. This yields the same p-value as the chi-squared test if we treat F as being F(h, ∞) distributed. In the cluster case, a finite-sample adjustment instead treats F as being F(h, G – 1) distributed. This is analogous to using the T(G – 1) distribution rather than N[0,1] for a test on a single coefficient.

In Stata, the finite-sample adjustment of using the T(G – 1) for a t-test on a single coefficient, and using the F(h,G – 1) for an F-test, is only done after OLS regression with command regress. Otherwise, Stata reports critical values and p-values based on the N[0,1] and χ²(h) distributions.

Thus, Stata does no finite-cluster correction for tests and confidence intervals following instrumental variables estimation commands, nonlinear model estimation commands, or even after command regress in the case of tests and confidence intervals using commands testnl and nlcom. The discussion in Section VI was limited to inference after OLS regression, but it seems reasonable to believe that for other estimators one should also base inference on the T(G – 1) and F(h,G – 1) distributions, and even then tests may overreject when there are few clusters.

Some of the few-cluster methods of Section VI can be extended to tests of more than one restriction following OLS regression. The Wald test can be based on the bias-adjusted variance matrices CR2VE or CR3VE, rather than CRVE. For a bootstrap with asymptotic refinement of a Wald test of H₀ : Rβ = r, in the b^th resample we compute . Extension of the data-determined degrees of freedom method of Section VID to tests of more than one restriction requires, at a minimum, extension of Theorem 4 of Bell and McCaffrey (2002) from the case that covers β, where β is a single component of β, to Rβ. An alternative ad hoc approach would be to use the distribution where is an average (possibly weighted by estimator precision) of v* defined in Equation 26 computed separately for each exclusion restriction.

For the estimators discussed in the remainder of Section VII, the rank of is again the minimum of G – 1 and the number of parameters (K). This means that at most G – 1 restrictions can be tested using a Wald test, in addition to the usual requirement that h ≤ K.

B. Instrumental Variables Estimators

The cluster-robust variance matrix estimate for the OLS estimator extends naturally to the IV estimator, the two-stage least squares (2SLS) estimator, and the linear GMM estimator.

The following additional adjustments must be made when errors are clustered. First, a modified version of the Hausman test of endogeneity needs to be used. Second, the usual inference methods when instruments are weak need to be adjusted. Third, tests of overidentifying restrictions after GMM need to be based on an optimal weighting matrix that controls for cluster correlation of the errors.

C. Nonlinear Models

For nonlinear models, there are several ways to handle clustering. We provide a brief summary; see Cameron and Miller (2011) for further details.

For concreteness, we focus on logit regression. Recall that in the cross-section case y_i takes value 0 or 1 and the logit model specifies that , where Λ(Z) = e^z / (1 + e^z).

D. Cluster-Randomized Experiments

Increasingly, researchers are gathering their own data, often in the form of field or laboratory experiments. When analyzing data from these experiments, they will want to account for the clustered nature of the data. And so when designing these experiments, they should also account for clustering. Fitzsimons et al. (2012) uses a wild cluster bootstrap in an experiment with 12 treated and 12 control clusters.

Traditional guidance for computing power analyses and minimum detectable effects (see Duflo, Glennerster, and Kremer 2007, pp. 3918–22; Hemming and Marsh 2013) are based on assumptions of either independent errors or, in a clustered setting, a random effects common-shock model. Ideally, one would account for more general forms of clustering in these calculations (the types of clustering that motivate cluster-robust variance estimation), but this can be difficult to do ex ante. If you have a data set that is similar to the one you will be analyzing later, then you can assign a placebo treatment and compute the ratio of cluster-robust standard errors to default standard errors. This can provide a sense of how to adjust the traditional measures used in design of experiments.

A. Individual-Level Cross-Section Data: One Sample

In our first application, we use data on 65,685 individuals from the year 2012. The model is

(33)

where y_ig is log-wage, d_g is a randomly generated dummy “policy” variable equal to one for one-half of the states and zero for the other half, and z_ig is a set of individual-level controls (age, age squared, and education in years). Estimation is by OLS and by FGLS controlling for state-level random effects.

The policy variable d_g is often referred to as a “placebo” treatment and should be statistically significant in only 5 percent of tests performed at significance level 0.05.

Table 1 reports the estimated coefficient of the policy variable, along with its standard error computed in several different ways, when there are 51 clusters (states).

OLS results given in the first column of Table 1 are obtained using Stata command regress. The default standard error is misleadingly small (se = 0.0042), leading to the dummy variable being very highly statistically significant (t = –0.0226 / 0.0042 = –5.42) even though it was randomly generated independently of log-wage. The White heteroskedastic-robust standard error, from regress option vce(robust), is similar in magnitude. From Section IVA, White standard errors should not be used if N_g is small but here N_g is large. The cluster-robust standard error (se = 0.0229) using option vce(cluster state) is 5.5 times larger and leads to the more sensible result that the regressor is statistically insignificant (t = –0.99). In results not presented in Table 1, the cluster-robust standard errors of the other regressors—age, age squared, and education—were, respectively, 1.2, 1.2, and 2.3 times the default. So ignoring clustering again understates the standard errors. As expected, a pairs cluster bootstrap (without asymptotic refinement) using option vce(boot, cluster(state)), yields very similar cluster-robust standard error.

Note that Formula 6 suggests that the cluster-robust standard errors are 4.9 times the default , close to the observed multiple of 5.5. Formula 6 may work especially well in this example as taking the natural logarithm of wage leads to model error that is close to homoskedastic and equicorrelation is a good error model for individuals clustered in regions.

FGLS estimates for a random effects model with error process defined in Equation 16 are given in the second column of Table 1. These were obtained using command xtreg, re after xtset state. The cluster-robust standard error defined in Equation 15 and computed using option vce(robust) is 0.0214/0.0199 = 1.08 times larger than the default. The pairs cluster bootstrap, implemented using option vce(boot), yields a similar cluster-robust standard error.

In principle, FGLS can be more efficient than OLS. In this example, there is a modest gain in efficiency with the cluster-robust standard error equal to 0.0214 for FGLS compared to 0.0229 for OLS.

Finally, to illustrate the potential pitfalls of pairs cluster bootstrapping for standard errors when there are few clusters, discussed in Section VIC, we examine a modification with only six states broken into treated (AZ, LA, MD) and control (DE, PA, UT). For these six states, we estimate a model similar to that in Table 1. Then with default se = 0.0128. We then perform a pairs cluster bootstrap with 999 replications. The bootstrap se = 0.0622 is similar to the cluster-robust se = 0.0577. However, several problems arise. First, 28 replications cannot be estimated, presumably due to no variation in treatment in the bootstrap samples. Second, a kernel density estimate of the bootstrapped s reveals that their distribution is very multimodal and has limited density near the middle of the distribution. Considering these results, we would not feel comfortable using the pairs cluster bootstrap in this data set with these few clusters. Better to base inference on a wild cluster bootstrap.

This example highlights the need to use cluster-robust standard errors even when model errors are only weakly correlated within cluster (the intraclass correlation of the residuals in this application is 0.018), if the regressor is substantially correlated within cluster (here perfectly correlated within cluster), and/or cluster sizes are large (ranging here from 519 to 5,866).

B. Individual-Level Cross-Section Data: Monte Carlo

We next perform a Monte Carlo exercise to investigate the performance of various cluster-robust methods as the number of clusters becomes small. The analysis is based on the same cross-section regression as in the previous subsection.

In each replication, we generate a data set by sampling (with replacement) states and all their associated observations. For quicker computation of the Monte Carlo simulation, within each state we use only a 20 percent subsample of individuals, so there are on average approximately 260 observations per cluster.

We explore the effect of the number of clusters G by performing varying simulations with G equal to 6, 10, 20, 30, or 50. Given a sample of states, we assign a dummy “policy” variable to one-half of the states. We run OLS regressions of log-wage on the policy variable and the same controls as used for the Table 1 regressions.

In these simulations, we perform tests of the null hypothesis that the slope coefficient of the policy variable is zero. Table 2 presents rejection rates that with millions of replications should equal 0.05, since we are testing a true hypothesis at a nominal 5 percent level. For G = 6 and 10, we perform 4,000 simulations, so we expect that 95 percent of these simulations will yield estimated test size in the range (0.043, 0.057) if the true test size is 0.05. For larger G, there are 1,000 simulations and the 95 percent simulation interval is instead (0.036, 0.064).

We begin with lengthy discussion of the many clusters case. These results are given in the final column (G = 50) of Table 2. Rows 1–9 report sizes for Wald tests based on where se is computed in various ways while Rows 10–15 report sizes for tests using various bootstraps with an asymptotic refinement. Basic Stata commands yield the standard errors in Rows 1–4 and 9, while the remaining rows require additional coding.

Row 1 presents the size of tests using heteroskedastic-robust standard errors, obtained using Stata command regress, vce(robust). Ignoring clustering leads to great overrejection due to considerable underestimation of the standard error. Using Formula 6 for this 20 percent subsample yields a standard error inflation factor of . So t = 1.96 using the heteroskedastic-robust standard error is really t = 1.96 / 2.38 = 0.82. And using standard normal critical values, an apparent p = 0.05 is really p = 0.41 since Pr[|z| > 0.82] = 0.41. This crude approximation is fairly close to p = 0.498 obtained in this simulation.

Results using cluster-robust standard errors, presented in Rows 2–4 and obtained from regress, vce(cluster state), show that even with 50 clusters the choice of distribution to use in obtaining critical value makes a difference. The rejection rate is closer to 0.05 when T(G – 1) critical values are used than when N[0,1] critical values are used. Using T(G – 2) in Row 4, suggested by the study of Donald and Lang (2007), leads to slight further improvement, but there is still overrejection.

Results using the bias adjustments CR2 and CR3 discussed in Section VIB, along with T(G – 1) critical values, are presented in Rows 5–6. Bias adjustment leads to further decrease in the rejection rates, toward the desired 0.05.

Rows 7 and 8 use critical values from the T distribution with the data-determined degrees-of-freedom of Section VID, equal to 17 on average when G = 50 (see Rows 14 and 17). This leads to further improvement in the Monte Carlo rejection rate.

Bootstrap standard errors obtained from a standard pairs cluster bootstrap, implemented using command regress, vce(boot, cluster(state)) are used in Row 9. For G = 50, the rejection rate is essentially the same as that in Row 3, as expected since this bootstrap has no asymptotic refinement.

Rows 10–15 implement the various percentile-t bootstraps with asymptotic refinement presented in Section VIC. Only 399 bootstraps are used here as any consequent bootstrap simulation error averages out over the many Monte Carlo replications. But if these bootstraps were used just once, as in an empirical study, a percentile-t bootstrap should use at least 999 replications. Row 10 can be computed using the bootstrap: command, see our posted code, while Rows 11–15 require additional coding. For G = 50, the various bootstraps give similar results with rejection rates of a bit more than 0.06.

Rows 16–19 give the effective number of clusters. The Imbens and Kolesar (2013) measure v* in Equation 26, denoted IK, and the Carter, Schnepel, and Steigerwald (2013) measure G* in Equation 27, denoted CSS, are both equal to 17 on average when G = 50. For the IK measure, across the 1,000 simulations the 5^th percentile is 9.6 and the 95^th percentile is 29.5.

We next examine settings with fewer clusters than G = 50. Then most methods lead to rejection rates even further away from the nominal test size of 0.05.

Consider the case G = 6. Rows 2–4 and 8 compute the same Wald test statistic but use different degrees of freedom in computing p-values. This makes an enormous difference when G is small, as the critical value for a Wald test at level 0.05 rises from 2.571 to 2.776 and 3.182 for, respectively, the T(5), T(4), and T(3) distributions, and from Row 16 the IK degrees of freedom averages 3.3 across the simulations. The CSS degrees of freedom is larger than the IK as, from Section VID, it involves an approximation that only disappears as G becomes large.

Using a bias-corrected CRVE also makes a big difference. It appears from Rows 6 and 7 that it is best to use the CR3 bias-correction with T(G – 1) critical values, and the CR2 bias-correction with T(v*) critical values where v* is the Imbens and Kolesar (2013) calculated degrees of freedom.

A downside to using cluster-robust standard errors is that they provide an estimate of the standard deviation of that is more variable than the default or heteroskedastic-robust standard errors. This introduces a potential bias—variance tradeoff. To see whether this increased variability is an issue, we performed 1,000 Monte Carlo replications using the full cross-section microdata set, resampling the 50 states with replacement. The standard deviation of the cluster-robust standard error across the 1,000 replications equaled 12.3 percent of the mean cluster-robust standard error, while the standard deviation of the heteroskedastic-robust standard error equaled 4.5 percent of its mean. So while the CRVE is less biased than heteroskedastic-robust (or default), it is also more variable. But the increased variability is relatively small, especially compared to the very large bias that can arise if clustering is not controlled for.

Rows 10–15 present various bootstraps with asymptotic refinement. From Row 10, the pairs cluster bootstrap performs extremely poorly for G ≤ 10.

Results for the wild cluster bootstrap using a Rademacher 2 point distribution are presented in Rows 11–13. From Section VIB, this bootstrap yields only 2⁶ = 64 possible data sets when G = 6, and hence at most 64 unique values for w*. This leads to indeterminacy for the test p-value. Suppose the p-value is in the range [a,b]. Then H₀ is rejected in Row 11 if a < 0.05, in Row 12 if (a+b)/2 < 0.05, and in Row 13 if b < 0.05. The indeterminacy leads to substantially different results for G as low as six, though not for G ≥ 10.

The wild cluster bootstrap using the Webb 6 point distribution, see Row 14, does not have this complication when G = 6. And it yields essentially the same results as those using the Rademacher 2 point distribution when G ≥ 10.

Comparing Row 15 to Row 12, we see that imposing the null hypothesis in performing the wild bootstrap does not change the rejection rate very much in this set of simulations when G ≥ 10, although it appears to matter when G = 6. By comparison, we have found a more substantial difference when simulating from the d.g.p. of Cameron et al. (2008).

In summary, the various wild cluster bootstraps lead to test size that is closer to 0.05 than using a standard Wald test with cluster-robust standard errors and T(G – 1) critical values. But the test size still exceeds 0.05 and the bias adjustments in Rows 6 and 7 appear to do better.

This example illustrates that at a minimum one should use a cluster-robust Wald test with T(G – 1) critical values. Especially when there are few clusters, it is better still to use bias-adjusted cluster-robust standard errors or to use a wild cluster bootstrap. In this Monte Carlo experiment, with few clusters the test size was closest to the nominal size using a Wald test with cluster-robust standard errors computed using the CR2 correction and T(v*) critical values, or using the larger CR2 correction with T(G – 1) critical values.

C. State-Year Panel Data: One Sample

We next turn to a panel difference-in-difference application motivated by Bertrand et al. (2004). The underlying data are again individual-level data from the CPS but now obtained for each of the years 1977 to 2012.

In applications where the policy regressor of interest is only observed at the state-year level, it is common to first aggregate the individual-level data to the state-year level before OLS regression. Several methods are used. We use the following method:

The model estimated for 51 states from 1977–2012 is

(34)

where is the average log-wage in year t and state s (after partialling out individual level covariates), α_s and δ_t are state and year dummies, and d_ts is a random “policy” variable that turns on and stays on for the last 18 years for one-half of the states. Here G = 51, T = 36, and N = 1836.

The individual level covariates (age, age squared, and years of education) are partialled out using a two-step estimation procedure presented in Hansen (2007b). Define D_ts to be state-by-year dummies. First, we OLS regress log wage (y_its) on state-by-year dummies D_ts and on the individual level covariates. And second, in Equation 34 equals the estimated coefficients of the D_ts dummies.

To speed up bootstraps, and to facilitate computation of the CR2 residual adjustment, we additionally partial out the state fixed effects and year fixed effects in Equation 34 by the standard Frisch-Waugh method. We separately regress and d_ts on the state dummies and the year dummies. Then β is estimated by regressing the residuals of on the residuals of d_ts with no constant. As noted below, regression using residuals leads to slightly different standard errors due to different degrees of freedom used in calculating the CRVE.

Table 3 presents results for the policy dummy regressor, which should have coefficient zero since it is randomly assigned.

We begin with Model 1, OLS controlling for state and year fixed effects. Using default or White-robust standard errors (Rows 1–2) leads to a standard error of 0.0037 that is much smaller than the cluster-robust standard error of 0.0119 (Row 3), where clustering is on state. Similar standard errors are obtained using the CR2 correction (Rows 4 and 5) and bootstrap without asymptotic refinement (Row 6). Note that from Rows 10 and 11 the IK and CSS degrees of freedom are calculated to be, respectively, G – 1 and G, an artifact of having balanced cluster sizes and a single regressor that is symmetric across clusters.

The inclusion of state and year fixed effects complicates computation of the degrees of freedom (df) adjustment used in computing the CRVE. The Row 3, Column 1 results are obtained from regression of residual on residual without intercept, using command regress, noconstant vce(cluster(state)). Then, from Formula 12 df = [G / (G −1)] × [GT / (GT −1)]. If instead we directly regressed log-wage on the state and year fixed effects and the K regressors, again using regress, vce(cluster(state)), then df = [G / (G −1)] × [GT / (GT − T – G −1)] and the cluster-robust standard error equals 0.0122 rather than 0.0119. And if instead we directly estimate the log-wage equation using xtreg, vce(robust) after xtset state, then df = G / (G −1) and the cluster-robust standard error equals 0.0120. In this example with large T and G, these adjustments make little difference, but for small G or T one should use df = G / (G −1) as explained in Section IIIB.

The corresponding p-values for tests of the null hypothesis that β = 0, following OLS regression, are given in Column 4 of Table 3. Default and heteroskedastic-robust standard errors lead to erroneously large t-statistics (of 0.0156 / 0.0037 = 4.22), so p = 0.000 and the null hypothesis is incorrectly rejected. Using various standard errors that control for clustering (Rows 3–6) leads to p ≃ 0.20 so that the null is not rejected. Rows 7–9 report p-values from several percentile-t bootstraps that again lead to rejection of H₀ : β = 0.

For illustrative purposes, we also compute standard errors allowing for two-way clustering, see Section V, with clustering on both state and year. These are computed using the user-generated Stata add-on program cgmreg.ado. Clustering on year is necessary if both the regressor and the model errors are correlated across states in a given year. For this application, the result (s.e. = 0.01167) is very similar to that from clustering on state alone (s.e. = 0.01185). In some other panel applications, the two-way cluster robust standard errors can be substantially larger than those from clustering on state alone.

The main lesson from the Model 1 OLS results is that even after inclusion of state fixed effects one needs to use cluster-robust standard errors that cluster on state. The inclusion of state fixed effects did not eliminate the within-state correlation of the error. In this example, the correct cluster-robust standard errors are 3.6 times larger than the default.

Model 2 again uses OLS estimation but drops the state fixed effects from Equation 34. Dropping these fixed effects leads to much less precise estimation as the cluster-robust standard error (Row 3) increases from 0.0119 to 0.0226 and this cluster-robust standard error is now 0.0226 / 0.0062 = 3.6 times the default, compared to a ratio of 3.6 when state fixed effects were included. Note that inclusion of state fixed effects (Model 1) did soak up some of the within-state error correlation as expected, but there still remained substantial within-cluster correlation of the error so that cluster-robust standard errors need to be used.

For Model 2, the comparisons of the various standard errors and p-values are qualitatively similar to those for Model 1 so are not discussed further.

Model 3 estimates the same model as Model 1, except that the state and year fixed effects are directly estimated, and estimation is now by FGLS allowing for an AR(1) process for the errors. Because there are 36 years of data, the bias correction of Hansen (2007b)—see Section IIIC—will make little difference and is not used here. Estimation uses Stata command xtreg, pa corr(ar 1) after xtset state.

Comparing Rows 1 and 3, again, even with inclusion of state fixed effects, one should obtain standard errors that cluster on state, using xtreg, pa option vce(robust)). The difference is not as pronounced as for OLS, with FGLS cluster-robust standard error that is 0.0084/0.0062 = 1.4 times the default.

FGLS estimation has led to substantial gain in efficiency, with cluster-robust standard error (Row 3) for FGLS of 0.0084 compared to 0.0119 for OLS.

This example illustrates that, even with state fixed effects included in a state year panel, inference should be based on cluster-robust standard errors. Furthermore, there can be substantial efficiency gains to estimating by FGLS rather than OLS.

D. State-Year Panel Data: Monte Carlo

We next perform a Monte Carlo exercise to investigate the performance of various cluster-robust methods as the number of clusters becomes small. The analysis is based on the same state-year panel regression as in the previous subsection, with each state-year observation based on log-wages after partialling out the individual level covariates.

In each simulation, we draw a random set of G states (with replacement), where G takes values 6, 10, 20, 30, and 50. When a state is drawn, we take all years of data for that state. We then assign our DiD “policy” variable to half the states with the policy turned on midway through the time period. In these simulations, we perform tests of the null hypothesis that the slope coefficient of the policy variable is zero. As in Table 2, for G = 6 and 10 we perform 4,000 simulations, so we expect that 95 percent of these simulations will yield estimated test size in the range (0.043, 0.057). For larger G, there are 1,000 simulations and the 95 percent simulation interval is instead (0.036, 0.064).

We begin with the last column of Table 4, with G = 50 states. All tests aside from that based on default standard errors (Row 1) have rejection rates that are not appreciably different from 0.05 once we allow for simulation error.

As the number of clusters decreases, it becomes clear that one should use the T(G – 1) or T(G – 2) distribution for critical values and even this leads to mild overrejection with low G. The pairs cluster percentile-t bootstrap fails with few clusters, with rejection rate of only 0.005 when G = 6. For low G, the wild cluster percentile-t bootstrap has similar results with either two-point or six-point weights with very slight overrejection.