Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
173,592 characters · 31 sections · 142 citation commands
Cluster-Robust Inference: A Guide to Empirical Practice
\thispagestyle{empty}
\setcounter{page}{1} \onehalfspacing
Ideally, the observations in a sample would be independent of each other and would each contribute roughly the same amount of information about the parameter(s) of interest. From the earliest days of econometrics, it has been recognized that this ideal situation often does not apply to time\kern 0.04167em-series data, because there may be serial correlation. But it has taken much longer for econometricians to realize that it generally does not apply to cross\kern 0.04167em-section data either. The first important step, following White_1980, was to allow for heteroskedasticity of unknown form, and for quite some time this was the default in empirical work that used cross\kern 0.04167em-section data. More recently, however, it has become common for investigators to drop the assumption of independence as well as the assumption of homoskedasticity.
There are many ways in which cross\kern 0.04167em-section observations might be dependent, and sometimes it is possible to model this dependence explicitly. For example, there is a large literature on spatial econometrics and statistics, in which each observation is associated with a point in space, and the correlation between any two observations is assumed to depend (usually in a rather simple parametric way) on the distance between them. See, among many others, Anselin_1988, \citet*{Gelfand_2010}, and Corrado_2012. However, there are a great many cases in which either the “distance” between any pair of observations cannot be measured, or the correlation between them is not related to distance in any way that can readily be modeled.
A more widely applicable approach, on which we focus in this paper, is to employ cluster-robust inference in the context of least-squares estimation. This approach has become increasingly popular over the past quarter century and is now used routinely in a great deal of empirical microeconomic work. The idea is to divide the sample into $G$ disjoint clusters. Depending on the nature of the data, the clusters might correspond to classrooms, schools, families, villages, hospitals, firms, industries, years, cities, counties, states, or countries. This list is by no means exhaustive. Any pattern of heteroskedasticity and/or dependence is allowed within each cluster, but it is assumed that there is independence across clusters and that the assignment of observations to clusters is known.
Under these assumptions, it is easy to compute cluster-robust standard errors that can be used to produce asymptotically valid inferences; see (ref). However, these inferences may not be at all reliable in finite samples. Hypothesis tests may reject far more often than they should. Less commonly, they may reject far less often. In consequence, the actual coverage of confidence intervals may differ greatly from their nominal coverage. Therefore, in practice, using cluster-robust inference often requires a good deal of care.
There are several recent survey papers on cluster-robust inference, including CM_2015, JGM-CJE, Esarey_2019, and MW-survey. \citet*{CGH_2018} surveys a broader class of methods for various types of dependent data. Although there will inevitably be some overlap with these papers, our aim is to provide a guide to empirical practice rather than a survey of the extant literature. We therefore apologize for any missing references and refer the reader to the survey papers just mentioned for more complete bibliographies. Our guide is closely based on the econometric theory and simulation evidence that is currently available. When the theory is clear and the evidence is strong, we make definitive recommendations for empirical practice. However, when the theory is less clear or the evidence is weak, our recommendations are more guarded.
This guide does not discuss models with clustered data estimated by instrumental variables (IV). For such models, neither the current state of econometric theory nor the available simulation evidence allows us to make recommendations with any confidence. The number of over-identifying restrictions and the strength of the instruments can greatly affect the reliability of finite\kern 0.04167em-sample IV inference, and dealing with these issues may often be even more important than dealing with the issues associated with clustering. There is an enormous literature on the topic of weak instruments; see \citet*{Andrews_2019} for a recent survey. That paper suggests that, when the disturbances of a regression model are independent and homoskedastic, it is generally possible to obtain reliable (although perhaps imprecise) inferences even when the instruments are quite weak. However, it also states that this is not the case, in general, when there is heteroskedasticity and/or clustering.
In (ref), we obtain the (true) variance matrix for the coefficient estimators in a linear regression model with clustered data. The form of this matrix depends on critical assumptions about the score vectors for each cluster. In practice, inference must be based on a cluster-robust variance estimator, or CRVE, which estimates the unknown variance matrix. We discuss the three CRVEs that are commonly encountered.
(ref) deals with the important and sometimes controversial issue of when to use cluster-robust inference. It also illustrates how complicated patterns of intra-cluster correlation can arise in the context of a simple factor model, introduces the concept of leverage at the cluster level, discusses the role of cluster fixed effects, and describes several procedures for deciding the level at which to cluster.
(ref) concerns the key issue of asymptotic inference. It explains how to obtain asymptotically valid inferences and discusses what determines how reliable, or unreliable, these inferences are likely to be in practice. In many cases, bootstrap inference tends to be more reliable than asymptotic inference. (ref) describes two methods for bootstrap inference, namely, the pairs cluster bootstrap and the restricted version of the wild cluster bootstrap, which is called the WCR bootstrap. The former has the advantage of being applicable to a wide variety of econometric models, while the latter is only applicable to regression models with clustered data, for which it typically performs better. For clustered linear regression models, both of these methods can be remarkably inexpensive to implement, even for very large samples. We recommend that, in most cases, the WCR bootstrap be among the methods employed for inference.
(ref) goes on to discuss some related inferential procedures. The first of these uses an alternative critical value estimated from the data, and the second is randomization inference, which can work well in certain cases where even the WCR bootstrap fails. (ref) discusses what an empirical investigator should report in order to convince the reader that results are reliable. (ref) presents an empirical example that uses individual data to study the effects of the minimum wage on the labor supply of teenagers. (ref) provides a summary of the main points of the paper. This is presented in the form of a short checklist or guide for empirical researchers on what to do in practice, with references to relevant sections.
Throughout the paper, we deal with the linear regression model $y_i = {\bm{x}}_i^\top\kern -0.08333em{\bm\beta} + u_i$, which, if the data have been divided into $G$ disjoint clusters, can be rewritten as
Here ${\bm{X}}_g$ is an $N_g\times k$ matrix of exogenous regressors, ${\bm\beta}$ is a $k$-vector of coefficients, ${\bm{y}}_g$ is an $N_g$-vector of observations on the regressand, and ${\bm{u}}_g$ is an $N_g$-vector of disturbances (or error terms). Thus ${\bm{X}}_g$, ${\bm{y}}_g$, and ${\bm{u}}_g$ stack the ${\bm{x}}_i^\top$, $y_i$, and $u_i$, respectively. In many cases, the regressors will include cluster fixed effects; see (ref). Since the \th{g} cluster has $N_g$ observations, the sample size is $N = \sum_{g=1}^G N_g$. The ${\bm{X}}_g$ may of course be stacked into an $N\times k$ matrix ${\bm{X}}$\kern -0.08333em, and likewise the ${\bm{y}}_g$ and ${\bm{u}}_g$ may be stacked into $N$-vectors ${\bm{y}}$ and ${\bm{u}}$, so that \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}} can be rewritten in the usual way as ${\bm{y}} = {\bm{X}}\!{\bm\beta} + {\bm{u}}$.
It is assumed that the data are generated by \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}} with ${\bm\beta}={\bm\beta}_0$. Under this assumption, the OLS estimator of ${\bm\beta}$ is
It follows that
where ${\bm{s}}_g = {\bm{X}}_g^\top {\bm{u}}_g$ denotes the $k\times1$ score vector corresponding to the \th{g} cluster. For a correctly specified model, ${\rm E}({\bm{s}}_g)={\bm{0}}$ for all $g$. From the rightmost expression in \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}}, the distribution of the OLS estimator $\hat{\bm\beta}$ depends on ${\bm{u}}$ only through the distribution of the score vectors ${\bm{s}}_g$. Ideally, the sum of the ${\bm{s}}_g$, suitably normalized, would be well approximated by a multivariate normal distribution with mean zero.
Because we can always divide the sample into $G$ clusters in any way we like, \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}} is true for any distribution of the disturbance vector ${\bm{u}}$. Dividing the sample into clusters only becomes meaningful if we further assume that
where the variance matrix of the scores for the \th{g} cluster, ${\bm{\Sigma}}_g$, is a $k\times k$ symmetric, positive semidefinite matrix. The second assumption in \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}} is the key one. It states that the scores for every cluster are uncorrelated with the scores for every other cluster. In contrast, the first assumption imposes no real limitations, so that the ${\bm{\Sigma}}_g$ matrices may display any patterns of heteroskedasticity and/or within-cluster dependence. Indeed, one motivation for using cluster-robust inference is that it is robust against both heteroskedasticity and intra-cluster dependence without imposing any restrictions on the (unknown) form of either of them.
For now, we will simply assume that \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}} holds for some specified division of the observations into clusters. Although the choice of clustering structure is often controversial, or at least somewhat debatable, the structure is almost always assumed known in both theoretical and applied work. The important issue of how to choose the clustering structure will be discussed below in (ref).
It follows immediately from \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}} that an estimator of the variance of $\hat{\bm\beta}$ should be based on the usual sandwich formula,
Of course, this matrix cannot be computed, because we need to estimate the ${\bm{\Sigma}}_g$. This can be done in several ways, as we discuss in (ref).
As \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}} makes clear, it is the properties of the score vectors that matter for inference. Of course, those properties are inherited from the properties of the disturbances and the regressors. If ${\bm{\Omega}}_g = {\rm E} ({\bm{u}}_g{\bm{u}}_g^\top |{\bm{X}})$ denotes the conditional variance matrix of ${\bm{u}}_g$, then
Thus, instead of making assumptions directly about the ${\bm{\Sigma}}_g$, as we did in \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}}, it may be more illuminating to make assumptions about the ${\bm{\Omega}}_g$ and the ${\bm{X}}_g$. If ${\rm E}({\bm{u}}_g{\bm{u}}_{g'}^\top |{\bm{X}}) = {\bm{0}}$ for all $g' \neq g$, then the second assumption in \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}} will hold. It will also hold if the regressors are exogenous and uncorrelated across clusters even when the disturbances are not.
Since the score vector ${\bm{s}}_g$ can be written as $\sum_{i=1}^{N_g} {\bm{s}}_{gi} = \sum_{i=1}^{N_g} {\bm{X}}^\top_{gi} u_{gi}$, where ${\bm{X}}_{gi}$ is the \th{i} row of ${\bm{X}}_g$ and $u_{gi}$ is the \th{i} element of ${\bm{u}}_g$, the outer product of the score vector with itself is seen to be
When ${\rm E}(u_{gi}^2 | {\bm{X}}) = \sigma^2$ and ${\rm E}(u_{gi} u_{gj} | {\bm{X}}) = 0$ for $i \neq j$, then ${\rm E} ( {\bm{s}}_g{\bm{s}}_g^\top | {\bm{X}} ) = \sigma^2 {\bm{X}}_g^\top\kern -0.08333em {\bm{X}}_g$. In that case, we would replace ${\bm{\Sigma}}_g$ with $\sigma^2 ({\bm{X}}_g^\top{\bm{X}}_g )$ in \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}} and obtain the classic result that $\operatorname{Var} (\hat{\bm\beta} | {\bm{X}} ) = \sigma^2({\bm{X}}^\top {\bm{X}})^{-1}$.
Taking expectations in \hyperref[{eq:opscores}]{\tagform@{\ref*{eq:opscores}}} and defining the covariance matrix ${\bm{\Sigma}}_{g,ij} = {\rm E} ( {\bm{s}}_{gi}{\bm{s}}^\top_{gj} )$, we find that, in general, ${\bm{\Sigma}}_g = \sum_{i=1}^{N_g}\sum_{j=1}^{N_g}{\rm E} ({\bm{s}}_{gi} {\bm{s}}_{gj}^\top) =\sum_{i=1}^{N_g}\sum_{j=1}^{N_g}{\bm{\Sigma}}_{g,ij}$. In the special case where the score vectors ${\bm{s}}_{gi}$ are uncorrelated within each cluster, i.e.\ where ${\bm{\Sigma}}_{g,ij}={\bm{0}}$ for $i \neq j$, we find that ${\bm{\Sigma}}_g = \sum_{i=1}^{N_g}{\rm E} ({\bm{s}}_{gi}{\bm{s}}_{gi}^\top) = \sum_{i=1}^{N_g} {\bm{\Sigma}}_{g,ii}$. The difference between these two expressions for ${\bm{\Sigma}}_g$ is
The rightmost expression in \hyperref[{eq:diffop}]{\tagform@{\ref*{eq:diffop}}} is just the summation of the $N_g^2 - N_g$ matrices that correspond to the off-diagonal elements of ${\bm{\Sigma}}_g$. It equals zero whenever there is no intra-cluster correlation, but in general it is $O(N_g^2)$. Therefore, incorrectly assuming that the scores are not correlated within clusters potentially leads to much larger errors of inference when clusters are large than when they are small. For sufficiently large values of $N_g$, these errors may be large even when all of the ${\bm{\Sigma}}_{g,ij}$ for $i\ne j$ are very small JGM_2016.
The famous “Moulton factor” Moulton_1986 gives the ratio of the true variance of an OLS coefficient, from \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}}, to the variance based on the classic formula $\sigma^2({\bm{X}}^\top{\bm{X}})^{-1}$ under the assumption that both the disturbances and the regressor of interest (after other regressors have been partialed out) are equi-correlated within clusters; see (ref). If the scores were scalars with intra-cluster correlation $\rho_s$, and the cluster sizes were constant, say $N_g = M$\kern -0.08333em, then the Moulton factor would be $1+(M-1)\rho_s$. The second term is proportional to the number of observations per cluster, so the mistakes made by not clustering can be enormous when clusters are large.
Since the disturbances in \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}} are neither independent nor homoskedastic, it seems relevant to consider GLS estimation, even though OLS estimation is almost always used in practice. If we were willing to specify a simple parametric form for the ${\bm{\Omega}}_g$ matrices, then we could use feasible GLS. For example, if we assumed that the disturbances were equi-correlated within each cluster, that would be equivalent to specifying a random-effects model; see (ref). In practice, however, the regressors in \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}} very often include cluster fixed effects ((ref)), and the latter remove whatever intra-cluster correlation a random-effects specification induces. So we would need to specify a more complicated model if we wanted to use feasible GLS. In any case, specifying a parametric form for the intra-cluster correlations would imply making assumptions much stronger than those in \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}}, and this would violate the principal objective of cluster-robust inference, namely, to be robust to arbitrary and unknown dependence and heteroskedasticity within clusters.
The natural way to estimate \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}} is to replace the ${\bm{\Sigma}}_g$ matrices by their empirical counterparts, which are the outer products of the empirical score vectors $\hat{\bm{s}}_g = {\bm{X}}_g^\top \hat{\bm{u}}_g$ with themselves. If, in addition, we multiply by a correction for degrees of freedom, we obtain
At present, this is by far the most widely used CRVE in practice. Observe that, when $G=N$\kern -0.08333em, CV$_{\kern -0.08333em1}$ reduces to the familiar HC$_1$ estimator MW_1985 that is robust only to heteroskedasticity of unknown form.
The empirical score vectors $\hat{\bm{s}}_g$ are not always good estimators of the ${\bm{s}}_g$. CV$_{\kern -0.08333em1}$ attempts to compensate for this by including a degrees-of-freedom factor. Two alternative CRVEs, proposed in BM_2002, instead replace the empirical score vectors $\hat{\bm{s}}_g$ by modified score vectors that use transformed residuals. The first of these is
where $\grave{\bm{s}}_g = {\bm{X}}_g^\top{\bm{M}}_{gg}^{-1/2}\kern 0.08333em\hat{\bm{u}}_g$, with ${\bm{M}}_{gg} = {\bf I}_{N_g} - {\bm{X}}_g({\bm{X}}^\top{\bm{X}})^{-1}\kern -0.08333em{\bm{X}}_g^\top$. Thus ${\bm{M}}_{gg}$ is the \th{g} diagonal block of the projection matrix ${\bm{M}}_{\bm{X}}$, which satisfies $\hat{\bm{u}} = {\bm{M}}_{\bm{X}}{\bm{u}}$, and ${\bm{M}}_{gg}^{-1/2}$ is its inverse symmetric square root. The CV$_2$ estimator reduces to the familiar HC$_2$ estimator when $G=N$. If the variance matrix of every ${\bm{u}}_g$ were proportional to an identity matrix, then CV$_{\kern -0.08333em2}$ would actually be unbiased PT_2018.
The second alternative CRVE is
where $\acute{\bm{s}}_g = {\bm{X}}_g^\top{\bm{M}}_{gg}^{-1}\kern 0.08333em\hat{\bm{u}}_g$. As we discuss in (ref), CV$_{\kern -0.08333em3}$ is actually a jackknife estimator which generalizes the familiar HC$_3$ estimator of MW_1985.
As written in \hyperref[{eq:CV2}]{\tagform@{\ref*{eq:CV2}}} and \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}}, both CV$_{\kern -0.08333em2}$ and CV$_{\kern -0.08333em3}$ are computationally infeasible for large samples, because they involve the $N_g\times N_g$ matrices ${\bm{M}}_{gg}$. However, \citet*{NAAMW_2020} proposes a more efficient algorithm for both of them, and \citet*{MNW-influence} provides an even more efficient one for CV$_{\kern -0.08333em3}$ by exploiting the fact that it is a jackknife estimator; see (ref). Moreover, when one of the regressors is a fixed-effect dummy for cluster $g$, the ${\bm{M}}_{gg}$ matrices are singular. This problem can be avoided, and some computer time saved, by partialing out the fixed-effect dummies as discussed in (ref) and PT_2018.
It seems plausible that both CV$_{\kern -0.08333em2}$ and CV$_{\kern -0.08333em3}$ should perform better in finite samples than CV$_{\kern -0.08333em1}$, because the modified score vectors $\grave{\bm{s}}_g$ and $\acute{\bm{s}}_g$ ought to provide better approximations to the ${\bm{s}}_g$ than do the $\hat{\bm{s}}_g$. We would expect tests based on CV$_{\kern -0.08333em3}$ to be more conservative than ones based on CV$_{\kern -0.08333em2}$, just as ones based on HC$_3$ are more conservative than ones based on HC$_2$, because the $\acute{\bm{s}}_g$ are “shrunk” more than the $\grave{\bm{s}}_g$. Simulation evidence dating back to BM_2002 suggests that CV$_{\kern -0.08333em3}$ typically yields the most reliable tests, but that they can sometimes under-reject; see also \citet*{MNW-bootknife}.
We cannot hope to obtain reliable inferences when using clustered data unless we know the actual clustering structure, at least to a good approximation. Thus, before specifying any clustering structure, we need to think about how intra-cluster correlations may arise and why independence across clusters may, or may not, be plausible for that structure.
\citet*{AAIW_2017} distinguishes between two alternative approaches to inference, referred to as “model-based” and “design-based.” The model-based approach is the traditional one, according to which every sample is treated as a random outcome, or realization, from some data-generating process (DGP), which in our context is the model \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}, and the objective is to draw inferences about parameters in the DGP, which are interpreted as features of the population. In particular, the DGP is the source of randomness and hence an important determinant of the clustering structure.
The “design-based” approach to inference, which is analyzed in detail in AAIW_2017, is conceptually different from the model-based framework. It involves thinking about the population, estimand, sampling scheme, parameters of interest, and even the notion of statistical uncertainty in a different manner. For example, according to the design-based approach, statistical uncertainty does not derive from a DGP, but rather is induced solely by the sampling uncertainty coming from sampling from a fixed population AAIW_2020. Therefore, under this approach, the sampling process is an important determinant of the clustering structure.
In some cases, the two approaches make use of similar inferential procedures, but in others they employ quite different ones. Both approaches may be informative about the choice of whether to cluster and at which level, although their motivations may differ. We follow most of the existing literature and focus exclusively on the model-based approach in the remainder of this paper.
Intra-cluster correlations of the disturbances and regressors, and hence of the scores, can arise for many reasons. By making the assumptions in \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}} and using cluster-robust inference, we avoid the need to model these correlations. Nevertheless, it can be illuminating to consider such models in order to learn about intra-cluster dependence and their consequences. The simplest and most popular model is the random-effects, or error-components, model
where $u_{gi}$ is the disturbance for observation $i$ within cluster $g$, $\varepsilon_{gi} \sim {\rm iid}(0, \omega^2)$ is an idiosyncratic shock for observation $i$, $\varepsilon_g \sim {\rm iid}(0, 1)$ is a cluster-wide shock for cluster $g$, and the two shocks are independent. The model \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}} implies that the variance matrix ${\bm{\Omega}}_g$ of the $u_{gi}$ for cluster $g$ has a very simple form with diagonal elements equal to $\lambda^2 + \omega^2$ and off-diagonal elements equal to $\lambda^2$. Thus the disturbances within every cluster are equi-correlated, with correlation coefficient $\lambda^2/(\lambda^2 + \omega^2)$.
Although the random-effects model \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}} has a long and distinguished history, it is almost certainly too simple. As we discuss in (ref), it is very common in modern empirical practice to include a set of cluster fixed effects among the regressors in \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}. In the case of \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}}, these fixed effects are simply estimates of the $\lambda\kern 0.04167em \varepsilon_g$, and by including them we therefore remove all of the intra-cluster correlation. Thus, if the random-effects model \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}} were correct, there would be no need to worry about cluster-robust inference whenever the regressors included cluster fixed effects. Note also that inclusion of cluster fixed effects usually comes at the price of larger standard errors on the coefficients of interest.
In practice, however, it usually seems to be the case that we need both cluster fixed effects and a CRVE; see (ref). This implies that whatever process is generating the intra-cluster correlations must be more complicated than \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}}. A simple example is the (very standard) factor model
which differs from \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}} in one important respect. The effect of the cluster-wide shock $\varepsilon_g$ on $u_{gi}$ is given not by a coefficient $\lambda$ but by a weight, or factor loading, $\lambda_{gi}$. These factor loadings could be either fixed parameters or random variables. They determine the extent to which observation $i$ within cluster $g$ is affected by the cluster-wide shock $\varepsilon_g$.
As an example, if the observations were for individual students, the clusters denoted classrooms, and the outcome were student achievement, then $\varepsilon_{gi}$ would measure unobserved student-specific characteristics, $\varepsilon_g$ would measure unobserved teacher quality (and perhaps other features of the class), and $\lambda_{gi}$ would measure the extent to which the disturbance term for student $i$ is affected by teacher quality. Clearly, the $\lambda_{gi}$ do not need to be the same for all $i$. Similar motivating examples based on \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}} can easily be given in many fields, including labor economics, health economics, development economics, and financial economics.
To verify that the factor model in \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}} generates dependence within clusters, it suffices to derive the second-order moments of the $u_{gi}$. We find that ${\rm E} (u_{gi}) = 0$ and $\operatorname{Var} (u_{gi}) = \lambda_{gi}^2 + \omega^2$. The cluster dependence is characterized by $\operatorname{Cov} (u_{gi},u_{gj}) = \lambda_{gi} \lambda_{gj}$, which differs across $(i,j)$ pairs and is zero only when the factor loadings are zero. In the context of the classroom example, the intra-cluster covariances would be zero only if the teacher had no effect on student achievement. Moreover, the correlations would be fully captured by classroom fixed effects if and only if $\lambda_{gi}$ were the same for all $i$.
The factor model in \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}} is discussed in terms of the disturbances $u_{gi}$ rather than the scores. There are at least two simple cases in which the same model structure, and in particular the same within-cluster correlation structure, applies to the scores. The first is when a regressor is generated by a model similar to \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}}, but possibly with different parameters. The second is when a regressor only varies at the cluster level, as is often the case for dummy variables, especially treatment dummies.
The model \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}} has only one clustering dimension, but the idea does not apply exclusively to cross\kern 0.04167em-section data. For example, if the observations also had a time dimension, we could replace each of the $\varepsilon_g$ by a time\kern 0.04167em-series process at the cluster level. This would yield a pattern where, within a cluster, observations that were closer together in time would be more correlated than observations that were further apart. In (ref), we generate placebo regressors in this way. For panel data, it is possible that, in addition to correlation within cross\kern 0.04167em-sectional units across time periods, there may be correlation within time periods across cross\kern 0.04167em-sectional units. This leads to two\kern 0.04167em-way clustering, which is discussed in (ref).
In principle, we might be able to estimate a factor model like \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}} and use it to obtain feasible GLS estimates, as mentioned at the end of (ref). However, this would be relatively complicated and rather arbitrary, since there are many plausible ways in which the details of \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}} could be specified. Moreover, any factor model would necessarily impose far stronger restrictions than the weak assumptions given in \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}}. Thus estimating any sort of factor model would inevitably require far more effort, and surely result in much more fragile inferences, than simply employing OLS estimation together with a CRVE.
Investigators very often include fixed effects at the cluster level among the regressors. There are generally good reasons for doing so. The cluster fixed effects implicitly model a large number of possibly omitted explanatory variables without assuming, implausibly, that the omitted variables are uncorrelated with the included ones.
It is sometimes believed that fixed effects remove any within-cluster dependence and hence eliminate the need to use a CRVE. However, as was pointed out by Arellano_1987, that is in fact only true under very special circumstances. Including cluster fixed effects in any regression model forces the intra-cluster sample average to be zero for each cluster. In particular, including cluster fixed effects transforms the factor model \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}} into
where the averages are taken across observations within each cluster, so that, for example, $\bar{u}_g = N_g^{-1}\sum_{i=1}^{N_g}u_{gi}$. The intra-cluster covariance for \hyperref[{eq:FEscore}]{\tagform@{\ref*{eq:FEscore}}} is
which is zero if and only if $\lambda_{gi}$ is the same for all $i$. In other words, the random-effects model \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}} is the only model within the class of factor models \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}} for which including cluster fixed effects can remove all intra-cluster dependence. Some dependence necessarily remains whenever there is any variation in factor loadings across observations within clusters.
Furthermore, \hyperref[{eq:FEcov}]{\tagform@{\ref*{eq:FEcov}}} strongly suggests that, whether or not a regression model includes cluster fixed effects, the scores will tend to be clustered whenever within-cluster dependence can be approximated by a factor model like \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}}. Including fixed effects will almost always reduce the intra-cluster correlations, but rarely will it entirely eliminate them. Because even very small intra-cluster correlations can have a large effect on standard errors when the clusters are large (see \hyperref[{eq:diffop}]{\tagform@{\ref*{eq:diffop}}} and the discussion that follows) it generally seems unwise to assume that cluster fixed effects make it unnecessary to use a CRVE.
In view of these arguments, it has become quite standard in modern empirical practice both to include cluster fixed effects (and perhaps other fixed effects as well) and also to employ cluster-robust inference. The empirical example in (ref) is typical in these respects. Of course, cluster fixed effects cannot be included when the regressor of interest is a treatment dummy and treatment is at the cluster level, since the treatment dummy and the fixed effects would be perfectly collinear. This problem does not arise for difference\kern 0.04167em-in-differences (DiD) regressions, because only some observations in the treated clusters are treated. In recent empirical work with non-staggered adoption of treatment, the regressions almost always include at least two sets of fixed effects, one for time periods and one for cross-sectional units, with clustering typically by the latter; see (ref).
In many cases, there is more than one level at which we could cluster. For example, with data on educational outcomes, we may be able to cluster by classroom, by school, or perhaps by school district. With data that are coded geographically, we may be able to cluster by county, by state, or even by region. Choosing the right level at which to cluster is not always easy, and choosing the wrong level can have serious consequences.
Suppose, for concreteness, that there are two possible levels of clustering, coarse and fine, with one or more fine clusters nested within each of the coarse clusters. When there are $G$ coarse clusters, the middle matrix in \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}} is $\sum_{g=1}^G {\bm{\Sigma}}_g$. If each coarse cluster contains $M_g$ fine clusters indexed by $h$, then ${\bm{\Sigma}}_g$ can be written as
where ${\bm{\Sigma}}_{g,h_1h_2}$ denotes the covariance matrix of the scores for fine clusters $h_1$ and $h_2$ within coarse cluster $g$. Under the assumption of fine clustering, ${\bm{\Sigma}}_{g,h_1h_2} = {\bm{\Sigma}}_{gh}$ when $h_1=h_2=h$ and ${\bm{\Sigma}}_{g,h_1h_2} = {\bm{0}}$ when $h_1\ne h_2$, so that the middle matrix in \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}} reduces to $\sum_{g=1}^G \sum_{h=1}^{M_g} {\bm{\Sigma}}_{gh}$.
From \hyperref[{eq:midcoarse}]{\tagform@{\ref*{eq:midcoarse}}}, the difference between the middle matrices for coarse and fine clustering is
The finest possible level of clustering is no clustering at all. In that case, the right-hand side of \hyperref[{eq:middiff}]{\tagform@{\ref*{eq:middiff}}} reduces to the right-hand side of \hyperref[{eq:diffop}]{\tagform@{\ref*{eq:diffop}}}, because the fine clusters within each coarse cluster are just the individual observations.
Under the assumption of fine clustering, the terms on the right-hand side of \hyperref[{eq:middiff}]{\tagform@{\ref*{eq:middiff}}} are all equal to zero. Under the assumption of coarse clustering, however, at least some of them are non-zero, and \hyperref[{eq:middiff}]{\tagform@{\ref*{eq:middiff}}} must therefore be estimated. If we cluster at the fine level when coarse clustering is appropriate, the CRVE is inconsistent. On the other hand, if we cluster at the coarse level when fine clustering is appropriate, the CRVE has to estimate \hyperref[{eq:middiff}]{\tagform@{\ref*{eq:middiff}}} even though it is actually zero. This makes the CRVE less efficient than it should be, leading to loss of power, or, equivalently, to confidence intervals that are unnecessarily long, especially when the number of coarse clusters is small.
Using simulation methods, MW-survey investigates the consequences on hypothesis tests of clustering at an incorrect level. Clustering at too fine a level generally leads to serious over-rejection, which becomes worse as the sample size increases with the numbers of clusters at all levels held constant. This is exactly what we would expect; see the discussion following \hyperref[{eq:diffop}]{\tagform@{\ref*{eq:diffop}}}. Clustering at too coarse a level also leads to both some over-rejection and some loss of power, especially when the number of clusters is small.
Two rules of thumb are commonly suggested for choosing the right level of clustering. The simplest is just to cluster at the coarsest feasible level CM_2015. This may be attractive when the number of coarse clusters $G$ is reasonably large, but it can be dangerous when $G$ is small, or when the clusters are heterogeneous in size or other features; see (ref).
A more conservative rule of thumb is to cluster at whatever level yields the largest standard error(s) for the coefficient(s) of interest MHE_2008. This rule will often lead to the same outcome as the first one, but not always. When $G$ is small, cluster-robust standard errors tend to be too small, sometimes much too small ((ref)). Hence, the second rule of thumb is considerably less likely to lead to severe over-rejection than the first one. However, because it is conservative, it can lead to loss of power (or, equivalently, confidence intervals that are unnecessarily long).
When the regressor of interest is a treatment dummy, and the level at which treatment is assigned is known, then it generally makes sense to cluster at that level \citep*{BDM_2004}. If treatment is assigned by cluster, whether for all observations in each cluster or just for some of them, as in the case of DiD models, then the scores will be correlated within the treated clusters whenever there is any intra-cluster correlation of the disturbances. Thus it never makes sense to cluster at a level finer than the one at which treatment is assigned. If we are certain that clusters are treated at random, then it also does not make sense to cluster at a coarser level. However, when there are two or more possible levels of clustering, it may not be realistic to assume that treatments are independent across finer-level clusters within the coarser-level ones. Unless we are certain that this is actually the case, it may be safer to cluster at a coarser level than the one at which treatment was supposedly assigned. For example, it may make sense to cluster by school instead of by classroom even when treatment was supposed to be assigned by classroom.
Instead of using a rule of thumb, we can test for the correct level of clustering. The best-known such test is an ingenious but indirect one proposed in IM_2016. It requires the model to be estimated separately for every coarse cluster, something that is not possible when the regressor of interest is invariant within some clusters, as is typically the case for treatment models and DiD models. It is also invalid if the parameter of interest has different meanings for different clusters; see (ref). When it is valid, the test statistic compares the observed variation of the estimates across clusters with an estimate of what that variation would be if clustering were actually at a finer level.
\citet*{MNW-testing} proposes direct tests called score\kern 0.04167em-variance tests, which compare the variance of the scores for two nested levels of clustering. For example, when there is just one coefficient, the empirical analog of \hyperref[{eq:middiff}]{\tagform@{\ref*{eq:middiff}}} is a scalar that, when divided by the square root of an estimate of its variance, is asymptotically distributed as ${\rm N}(0,1)$. The null hypothesis is that the (true) standard errors are the same for fine and coarse clustering, and the (one\kern 0.04167em-sided) alternative is that they are larger for the latter than for the former. MNW-testing also proposes wild (cluster) bootstrap implementations of these tests to improve their finite\kern 0.04167em-sample properties.
Cai_2021 takes a different approach, proposing a test for the level of clustering that is based on randomization inference ((ref)). This test is designed for settings with a small number of coarse clusters and a small number of fine clusters within each of them.
It seems natural to cluster at the coarse level when a test rejects the null hypothesis, and to cluster at the fine level when it does not. However, choosing the level of clustering in this way is a form of pre\kern 0.04167em-testing, which can lead to estimators with distributions that are poorly approximated by asymptotic theory, even in large samples LP_2005. Using a pre\kern 0.04167em-test in this way will inevitably lead to over-rejection when there is actually coarse clustering, the standard errors for coarse clustering are larger than the ones for fine clustering, and the test incorrectly fails to reject the null hypothesis. Thus using such a test is less conservative than relying on the second rule of thumb discussed above. On the other hand, we may feel more comfortable with the second rule of thumb when it agrees with the outcomes of one or more tests for the level of clustering.
As will be explained in (ref), asymptotic inference depends on being able to apply laws of large numbers and central limit theorems to functions of the (empirical) score vectors. How well those theorems work depends on how homogeneous the score vectors are across clusters. When they are quite heterogeneous, asymptotic inference may be problematic; see (ref). It is therefore desirable to measure the extent of cluster-level heterogeneity.
Classic measures of observation-level heterogeneity are leverage and influence \citep*{BKW_1980,CH_1986}. These are generalized to cluster-level measures in MNW-influence. One possible consequence of heterogeneity is that the estimates may change a lot when certain clusters are deleted. When this is the case, a cluster is said to be influential. In order to identify individually influential clusters, first construct the matrices ${\bm{X}}_g^\top{\bm{X}}_g$ and the vectors ${\bm{X}}_g^\top{\bm{y}}_g$, for $g=1,\ldots,G$. Then
is the vector of least squares estimators when cluster $g$ is deleted. It should not be expensive to compute $\hat{\bm\beta}^{(g)}$ for every cluster using \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}}. Note, however, that we cannot partial out regressors other than cluster fixed effects (see below) prior to computing the $\hat{\bm\beta}^{(g)}$, because the latter would then depend indirectly on the observations for the \th{g} cluster.
When there is a parameter of particular interest, say $\beta_j$, then it will often be a good idea to report the $\hat\beta_j^{(g)}$ for $g=1,\ldots,G$ in either a histogram or a table. If $\hat\beta_j^{(h)}$ differs a lot from $\hat\beta_j$ for some cluster $h$\kern -0.08333em, then cluster $h$ is evidently influential. In a few extreme cases, there may be a cluster $h$ for which it is impossible to compute $\hat\beta_j^{(h)}$\kern -0.08333em. If so, then the original estimates should probably not be believed. This will happen, for example, when cluster $h$ is the only treated one, and we will see in (ref) that inference is extremely unreliable in that case.
The $\hat{\bm\beta}^{(g)}$ are of interest even when there is no reason to expect any clusters to be influential. As MNW-bootknife shows, an alternative way to write CV$_{\kern -0.08333em3}$ is
This is the matrix version of the classic jackknife variance estimator given in Efron_81 and others. Unless all clusters are very small, \hyperref[{eq:jackvar}]{\tagform@{\ref*{eq:jackvar}}} is enormously faster to compute than \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}}. The Stata command summclust, which is described in detail in MNW-influence, calculates CV$_{\kern -0.08333em3}$ standard errors based on \hyperref[{eq:jackvar}]{\tagform@{\ref*{eq:jackvar}}}.
As pointed out in BKW_1980 and CH_1986, it is often valuable to identify high-leverage observations as well as influential ones. It is perhaps even more valuable to identify high-leverage clusters MNW-influence. Loosely speaking, a high-leverage cluster is one whose regressors contain a lot of information. At the observation level, high-leverage observations are associated with a high value of $h_i$, the \th{i} diagonal element of ${\bm{H}} = {\bm{P}}_{\bm{X}} = {\bm{X}} ({\bm{X}}^\top{\bm{X}})^{-1}{\bm{X}}^\top$. The analog of $h_i$ in the cluster case is the $N_g\times N_g$ matrix ${\bm{H}}_g = {\bm{X}}_g({\bm{X}}^\top{\bm{X}})^{-1} {\bm{X}}_g^\top$\kern -0.08333em. Since it is not feasible to report the ${\bm{H}}_g$, we suggest that investigators instead report their traces, which are
These are easy to compute because we have already calculated $({\bm{X}}^\top{\bm{X}})^{-1}$ and the ${\bm{X}}_g^\top{\bm{X}}_g$. For any cluster that contains just one observation, $L_g$ reduces to the usual measure of leverage at the observation level. High-leverage clusters can be identified by comparing the $L_g$ to their own average, which is $k/G$. If, for some $h$, $L_h$ is substantially larger than $k/G$, then cluster $h$ has high leverage. This can happen either because $N_h$ is much larger than $G/N$ or because the matrix ${\bm{X}}_h$ is somehow extreme relative to the other ${\bm{X}}_g$ matrices, or both. For example, $L_h$ is likely to be much larger than $k/G$ if cluster $h$ is one of just a few treated clusters.
Regression models often include cluster fixed effects. It is computationally attractive to partial them out before estimation begins, using for example the areg procedure in Stata. When one of the regressors is a fixed-effect dummy for cluster $g$, the matrices ${\bm{X}}^\top{\bm{X}} - {\bm{X}}_g^\top{\bm{X}}_g$ are singular. However, the problem solves itself if we partial out the fixed-effect dummies and replace ${\bm{X}}$ by $\tilde{\bm{X}}$ and ${\bm{y}}$ by $\tilde{\bm{y}}$, the matrix and vector of deviations from cluster means. For example, the \th{gj} element of $\tilde{\bm{y}}$ is $y_{gj} - N_g^{-1}\sum_{i=1}^{N_g} y_{gi}$. Since this depends only on observations for cluster $g$, the jackknife CV$_{\kern -0.08333em3}$ estimator \hyperref[{eq:jackvar}]{\tagform@{\ref*{eq:jackvar}}} remains valid.
In (ref), we discuss what quantities investigators should report in any empirical analysis that involves cluster-robust inference. In addition to measures of influence and leverage, these may include measures of partial leverage (the analog of leverage for a single coefficient) and summary statistics based on either leverage, partial leverage, or the effective number of clusters. An example is provided in (ref).
An interesting way to assess the validity of alternative standard errors is to run “placebo regressions.” The idea, first suggested in BDM_2004, is to start with a model and dataset, then generate a completely artificial regressor at random, add it to the model, and perform a $t$-test of significance. This is repeated a large number of times, and the rejection frequency is observed. The artificial regressor is often a dummy variable that is referred to as a “placebo law” or “placebo treatment.” Using such a dummy variable is natural because, for any level of intra-cluster correlation of the disturbances, the intra-cluster correlation of the scores is greatest for regressors that do not vary within clusters. However, any artificial regressor that is not completely uncorrelated within clusters can potentially be used.
Because a placebo regressor is artificial, we would expect valid significance tests at level $\alpha$ to reject the null close to $\alpha\%$ of the time when the experiment is repeated many times. Following the lead of BDM_2004 by using models for log-earnings based on age, education, and other personal characteristics, together with data taken from the Current Population Survey, several papers \citep*{JGM_2016,MW-JAE,Brewer_2018} find that not clustering, or clustering at below the state level, leads to rejection rates far greater than $\alpha$. In (ref), we find similar results for the datasets used in our empirical example. Our findings, and those of the papers cited above, all suggest that using a state\kern 0.04167em-level CRVE is important for survey data that samples individuals from multiple states. If we fail to do so, we will find, with probability much higher than $\alpha$, that nonsense regressors apparently belong in the model.
Since the empirical score vectors are $\hat{\bm{s}}_g = {\bm{X}}_g^\top\hat{\bm{u}}_g$, a placebo\kern 0.04167em-regressor experiment should lead to over-rejection whenever both the regressor and the residuals display intra-cluster correlation at a coarser level than the one at which the standard errors are clustered. As in (ref), suppose there are two potential levels of clustering, fine and coarse, with the fine clusters nested within the coarse clusters. If the placebo regressor is clustered at the coarse level, we would expect significance tests based on heteroskedasticity-robust standard errors to over-reject whenever the residuals are clustered at either level. Similarly, we would expect significance tests based on finely-clustered standard errors to over-reject whenever the residuals are clustered at the coarse level. (ref) in (ref) displays both of these phenomena.
Placebo regressions can provide useful guidance as to the correct level of clustering. However, using the rejection rates for placebo regressions with different levels of clustering as informal tests is really a form of pre\kern 0.04167em-testing. Thus, like using the formal tests discussed in (ref), doing this seems very likely to yield less conservative inferences than simply relying on the second rule of thumb.
Up to this point, we have assumed that there is clustering in only one dimension. However, there could well be clustering in two or more dimensions. With data that have both a spatial and a temporal dimension, there may be clustering by jurisdiction and also by time period. In finance, there is often clustering by firm and by year. Thus, instead of \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}, we might have
where the vectors ${\bm{y}}_{gh}$ and ${\bm{u}}_{gh}$ and the matrix ${\bm{X}}_{gh}$ contain, respectively, the rows of ${\bm{y}}$, ${\bm{u}}$, and ${\bm{X}}$ that correspond to both the \th{g} cluster in the first clustering dimension and the \th{h} cluster in the second one. The $GH$ clusters into which the data are divided in \hyperref[{eq:modelgh}]{\tagform@{\ref*{eq:modelgh}}} represent the intersection of the two clustering dimensions.
If there are $N_g$ observations in the \th{g} cluster for the first dimension, $N_h$ observations in the \th{h} cluster for the second dimension, and $N_{gh}$ observations in the \th{gh} cluster for the intersection, the number of observations in the entire sample is $N=\sum_{g=1}^G N_g=\sum_{h=1}^H N_h=\sum_{g=1}^G \sum_{h=1}^H N_{gh}$, where $N_{gh}$ might equal 0 for some values of $g$ and $h$. The scores for the clusters in the first dimension are ${\bm{s}}_g={\bm{X}}_g^\top{\bm{u}}_g$, for the clusters in the second dimension ${\bm{s}}_h={\bm{X}}_h^\top{\bm{u}}_h$, and for the intersections ${\bm{s}}_{gh}={\bm{X}}_{gh}^\top{\bm{u}}_{gh}$. If, by analogy with \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}}, we assume that
then the variance matrix of the scores is seen to be
The last condition in \hyperref[{eq:varmats}]{\tagform@{\ref*{eq:varmats}}} means that the scores are assumed to be independent whenever they do not share a cluster along either dimension. The third term in \hyperref[{eq:scorevar}]{\tagform@{\ref*{eq:scorevar}}} must be subtracted in order to avoid double counting. It is important to distinguish between two\kern 0.04167em-way clustering and clustering by the intersection of the two dimensions. If we assumed the latter instead of the former, then all three terms on the right-hand side of \hyperref[{eq:scorevar}]{\tagform@{\ref*{eq:scorevar}}} would be equal, and consequently ${\bm{\Sigma}} = \sum_{g=1}^G \sum_{h=1}^H {\bm{\Sigma}}_{gh}$. Thus these assumptions are radically different.
An estimator of the variance matrix of $\hat{\bm\beta}$ is
Here $\hat{\bm{\Sigma}}$ is an estimator of \hyperref[{eq:scorevar}]{\tagform@{\ref*{eq:scorevar}}}, with the empirical scores defined in the usual way; for example, $\hat{\bm{s}}_g = {\bm{X}}_g^\top\hat{\bm{u}}_g$. In practice, each of the matrices on the right-hand side of the second equation in \hyperref[{eq:betavar}]{\tagform@{\ref*{eq:betavar}}} is usually multiplied by a scalar factor, like the one in \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}}, designed to correct for degrees of freedom. Because the third term is subtracted, the matrix $\hat{\bm{\Sigma}}$ may not always be positive definite. This problem can be avoided by omitting the third term, which is asymptotically valid under some assumptions \citep*{Davezies_2021,MNW_2021}. Another possibility is to use an eigenvalue decomposition \citep*{CGM_2011}, although this merely forces the variance matrix to be positive semidefinite.
The idea of two\kern 0.04167em-way clustering can, of course, be generalized to three\kern 0.04167em-way clustering, four-way clustering, and so on. However, the algebra rapidly becomes daunting. If there were three clustering dimensions, for example, the analog of \hyperref[{eq:scorevar}]{\tagform@{\ref*{eq:scorevar}}} would have seven terms.
Two\kern 0.04167em-way clustering seems to have been suggested first in MH_2006 and rediscovered independently by \citet*{CGM_2011} and Thompson_2011. Although two\kern 0.04167em-way clustering has been widely used in empirical work, the asymptotic theory to justify it is much more challenging than the theory for the one\kern 0.04167em-way case, and this theory is still under active development \citep*{Chiang_2020,Chiang_2021JBES,Davezies_2021,MNW_2021,Menzel_2021}. In view of this, and because of the technical difficulties involved, we will focus mainly on one\kern 0.04167em-way clustering in the remainder of the paper.
For the regression model \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}, inference is commonly based on the $t$-statistic,
where the hypothesis to be tested is ${\bm{a}}^\top{\bm\beta}={\bm{a}}^\top{\bm\beta}_0$, with ${\bm{a}}$ a known $k$-vector. Here $\hat{\bm{V}}$ denotes one of CV$_{\kern -0.08333em1}$, CV$_{\kern -0.08333em2}$, or CV$_{\kern -0.08333em3}$, given in \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}}, \hyperref[{eq:CV2}]{\tagform@{\ref*{eq:CV2}}}, and either \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}} or \hyperref[{eq:jackvar}]{\tagform@{\ref*{eq:jackvar}}}, respectively. In many cases, just one element of ${\bm{a}}$, say the \th{j}, equals 1, and the remaining elements equal 0, so that \hyperref[{eq:tstat}]{\tagform@{\ref*{eq:tstat}}} is simply $\hat\beta_j - \beta_{j0}$ divided by its standard error. When there are $r>1$ linear restrictions, which can be written as ${\bm{R}}{\bm\beta} = {\bm{r}}$ with ${\bm{R}}$ an $r\times k$ matrix, inference can be based on the Wald statistic,
Of course, when $r=1$, the $t$-statistic \hyperref[{eq:tstat}]{\tagform@{\ref*{eq:tstat}}} is just the signed square root of a particular Wald statistic with ${\bm{R}}={\bm{a}}^\top$ and ${\bm{r}} = {\bm{a}}^\top{\bm\beta}_0$.
By letting the sample size become arbitrarily large, one can frequently obtain a tractable asymptotic distribution for any test statistic of interest, including \hyperref[{eq:tstat}]{\tagform@{\ref*{eq:tstat}}} and \hyperref[{Waldstat}]{\tagform@{\ref*{Waldstat}}}. Ideally, this would provide a good approximation to the actual distribution. With clustered data, there is more than one natural way to let the sample size become large, because we can make various assumptions about what happens to $G$ and the $N_g$ as we let $N$ tend to infinity. Which assumptions it is appropriate to use, and how well the resulting approximations work, will depend on the characteristics of the sample and the (unknown) DGP.
In order for inferences based on the statistics \hyperref[{eq:tstat}]{\tagform@{\ref*{eq:tstat}}} and \hyperref[{Waldstat}]{\tagform@{\ref*{Waldstat}}} to be asymptotically valid, two key asymptotic results must hold. First, a central limit theorem (CLT) must apply to the sum of the score vectors ${\bm{s}}_g$ in \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}}. In the limit, after appropriate normalization, the vector $\sum_{g=1}^G {\bm{s}}_g$ needs to follow a multivariate normal distribution with variance matrix $\sum_{g=1}^G {\bm{\Sigma}}_g$. Second, again after appropriate normalization, a law of large numbers (LLN) must apply to the matrices $\sum_{g=1}^G \hat{\bm{s}}_g\hat{\bm{s}}_g^\top$, $\sum_{g=1}^G \grave{\bm{s}}_g\grave{\bm{s}}_g^\top$, or $\sum_{g=1}^G \acute{\bm{s}}_g\acute{\bm{s}}_g^\top$ in the middle of the variance matrix estimators \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}}, \hyperref[{eq:CV2}]{\tagform@{\ref*{eq:CV2}}}, or \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}}, so that they converge to $\sum_{g=1}^G {\bm{\Sigma}}_g$. We refer to “appropriate normalization” here rather than specifying the normalization factors explicitly because, with clustered data, the issue of normalization is a very tricky one; see (ref). For asymptotic inference to be reliable, we need both the CLT and the LLN to provide good approximations.
There are currently two quite different types of assumptions on which the asymptotic theory of cluster-robust inference can be based. The most common approach, and we believe usually the most appropriate one, is to let the number of clusters tend to infinity. We refer to this as the “large number of clusters” approach and discuss it in (ref). An alternative approach is to hold the number of clusters fixed and let the number of observations within each cluster tend to infinity. We refer to this as the fixed-$G$ or “small number of large clusters” approach and discuss it in (ref). Some of the material in (ref) is quite technical, but it helps to explain when and why asymptotic inference can fail.
Inference based on asymptotic theory often performs well, but it can perform poorly in some commonly-encountered situations that are discussed in (ref). We therefore do not recommend relying only on CV$_{\kern -0.08333em1}$ and asymptotic theory. Because the bootstrap methods to be discussed in (ref) can work much better than asymptotic methods when the latter do not work well, we recommend that they be used almost all the time, at least to verify that both approaches yield similar results. In particular, we recommend using one or more variants of the wild cluster restricted, or WCR, bootstrap ((ref)) as a matter of routine.
The simplest assumption about how the sample size goes to infinity is that every cluster has a fixed number of observations, say $M$\kern -0.08333em. Then $N=MG$, and both $N$ and $G$ go to infinity at the same rate. Thus the appropriate normalizing factor for the parameter estimator is either $\sqrt{G}$ or $\sqrt{N}$. In this case, it is not difficult to show that $\sqrt{G}(\hat{\bm\beta} - {\bm\beta}_0)$ is asymptotically multivariate normal with variance matrix equal to the probability limit of $G$ times the right-hand side of \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}}. Moreover, the latter can be estimated consistently by $G$ times the CV$_{\kern -0.08333em1}$, CV$_{\kern -0.08333em2}$, or CV$_{\kern -0.08333em3}$ matrices. The first proof for this case of which we are aware is in White_1984; see also Hansen_2007.
In actual samples, clusters often vary greatly in size, so it is usually untenable to assume that every cluster has the same number of observations. The assumption that $G$ is proportional to $N$ may be relaxed by allowing $G$ to be only approximately proportional to $N$\kern -0.08333em, so that $G/N$ is roughly constant as $N\to\infty$. This implies that all the clusters must be small. In this case, the quality of the asymptotic approximations is not likely to be harmed much by moderate variation in cluster sizes. If a sample has, say, 500 clusters that vary in size from 10 to 50 observations, we would expect asymptotic inference to perform well unless there is some other reason (unrelated to cluster sizes) for it to fail.
\citet*{DMN_2019} and HansenLee_2019 take a more flexible approach, with primitive conditions that restrict the variation in the $N_g$ relative to the sample size. These conditions allow some clusters to be “small” and others to be “large” in the sense that some but not all $N_g \to \infty$ as $N \to \infty$. Although a key assumption is that $G\to\infty$ (i.e., is “large”), the appropriate normalization factor for $\hat{\bm\beta} - {\bm\beta}_0$ is usually not $\sqrt{G}$. Instead, this factor depends in a complicated way on the regressors, the relative cluster sizes, the intra-cluster correlation structure, and interactions among these; some examples of different normalizing factors are given in the papers cited above. For this reason, the key result that the $t$-statistic defined in \hyperref[{eq:tstat}]{\tagform@{\ref*{eq:tstat}}} is asymptotically distributed as standard normal is derived assuming that the rate at which $\hat{\bm\beta} - {\bm\beta}_0$ tends to zero is unknown. Of course, this result also justifies using the $t(G-1)$ distribution, which is more conservative, is widely used, and is derived from the theory discussed in (ref).
The application of a CLT to $\sum_{g=1}^G {\bm{s}}_g$, appropriately normalized, requires a restriction on the amount of heterogeneity that is allowed. Otherwise, just a few clusters might dominate the entire sample in the limit, thus violating the Lindeberg or Lyapunov conditions for the CLT. The necessary restrictions on the heterogeneity of clusters may be expressed in terms of two key parameters. The first of these parameters is the number of moments that is assumed to exist for the distributions of ${\bm{s}}_{gi}$ (uniformly in $g$ and $i$). We denote this parameter by $\gamma > 2$. When more moments exist, the distributions of ${\bm{s}}_{gi}$ are closer to the normal distribution, and hence the sample will feature fewer outliers or other highly leveraged observations or clusters.
In the clustered regression model, the variance of the scores, which is often referred to as the Fisher information matrix, is given by $\mathcal{J}_N = \sum_{g=1}^G \operatorname{Var} ({\bm{s}}_g )$. When appropriately normalized, $\mathcal{J}_N$ converges to a nonzero and finite matrix $\mathcal{J}$\kern -0.08333em. The rate of convergence $\eta_N$ is defined implicitly by $\eta_N^{-1} \mathcal{J}_N \to \mathcal{J}$\kern -0.08333em. This rate is the second key parameter. One interpretation of $\eta_N$ can be found in \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}}, from which the stochastic order of magnitude of $\hat{\bm\beta} - {\bm\beta}_0$ is seen to be $O_P(\eta_N^{1/2}\!/N)$. In general, $\eta_N \geq N$, with the equality holding whenever there is no intra-cluster correlation. The larger the value of $\eta_N$, the more slowly does $\hat{\bm\beta}$ converge to ${\bm\beta}_0$.
The conditions required on the heterogeneity of clusters to apply a CLT can be stated in terms of the parameters $\gamma$ and $\eta_N$. Specifically, when expressed in our notation, Assumption 3 of DMN_2019 states the following condition:
Because $\eta_N = o ( N^2)$ for consistency of $\hat{\bm\beta}$, the condition in \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}} makes it clear that we cannot allow a single cluster to dominate the sample, in the sense that its size is proportional to $N$\kern -0.08333em. More generally, \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}} shows that there is a tradeoff between information accumulation and variation in cluster sizes, as measured by the largest cluster size. To interpret this tradeoff, we will consider three different implications of \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}}.
First, when $\gamma$ increases, so that more moments exist, condition \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}} becomes less strong. In particular, when the scores are nearly normally distributed, in the sense that all their moments exist, then $\gamma=\infty$, which implies that $-2\gamma/(2\gamma-2) = -1$. In this case, \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}} reduces to $\eta_N^{-1/2} \sup_g N_g \to 0$, so that the size of the largest cluster must increase more slowly than the square root of the rate at which the Fisher information matrix converges. When $\gamma<\infty$, so that there are fewer moments, then the rate at which $\sup_g N_g$ is allowed to increase becomes smaller.
Second, suppose that the scores are uncorrelated, or more generally that a CLT applies to $N_g^{-1/2}{\bm{s}}_g$, as assumed in \citet*{BCH_2011} and IM_2010,IM_2016; see (ref). In this case, $\operatorname{Var} ({\bm{s}}_g)= O(N_g)$, so that $\eta_N = N$. Condition \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}} is then $N^{-(\gamma-2)/(2\gamma-2)} \sup_g N_g \to 0$. If the scores are nearly normal, then it reduces further to $N^{-1/2} \sup_g N_g \to 0$, so that the size of the largest cluster must increase no faster than the square root of the sample size. Once again, the fewer moments there are, the more slowly $\sup_g N_g$ is allowed to increase.
Third, suppose that the scores are generated by the factor model in \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}}, or by the simpler random-effects model \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}}. Then $\operatorname{Var} ({\bm{s}}_g)=O(N_g^2)$. If, in addition, $\inf_g N_g$ and $\sup_g N_g$ are of the same order of magnitude, then $\eta_N = N \sup_g N_g$ and the condition \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}} collapses to $N^{-1} \sup_g N_g \to~0$, regardless of the number of finite moments.
One possibly surprising implication of the above discussion is that, when there is more intra-cluster correlation, so that $\eta_N$ is relatively large, then greater heterogeneity of cluster sizes is allowed. That is, a higher degree of intra-cluster correlation implies a faster rate of convergence, $\eta_N$, of the Fisher information matrix, which in turn allows a larger $\sup_g N_g$ in \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}}. The intuition is that greater intra-cluster correlation reduces the effective cluster size, as measured by the amount of independent information a cluster contains. In the extreme case in which all observations in the \th{g} cluster are perfectly correlated, the size of the cluster is effectively 1 and not $N_g$. Note, however, that large clusters are implicitly weighted more heavily than small clusters even in this extreme case.
Although it is impossible to verify condition \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}} in any finite sample, investigators can always observe the $N_g$. The discussion of \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}} above suggests that asymptotic inference tends to be unreliable when the $N_g$ are highly variable, especially when a very few clusters are unusually large. This is exacerbated when the distribution of the data is heavy-tailed (has fewer moments), but is mitigated when the clusters have an approximate factor structure. Circumstances in which asymptotic inference can be unreliable are discussed in more detail in (ref).
A few authors have assumed that $G$ remains fixed (i.e., is “small”) as $N\to\infty$, while the cluster sizes diverge (i.e., are “large”). Notably, BCH_2011 proves that, for CV$_{\kern -0.08333em1}$, the $t$-statistic \hyperref[{eq:tstat}]{\tagform@{\ref*{eq:tstat}}} follows the $t(G-1)$ distribution asymptotically. An analogous result for the Wald statistic \hyperref[{Waldstat}]{\tagform@{\ref*{Waldstat}}} is discussed in (ref). These results provide useful approximations, but they are proven under some very strong assumptions. In particular, all the clusters are assumed to be the same size $M$\kern -0.08333em. In addition, the pattern of dependence within each cluster is assumed to be such that a CLT applies to the normalized score vectors $M^{-1/2}{\bm{s}}_g$ for all $g=1,\ldots,G$, as $M\to\infty$.
This second assumption is crucial, as it limits the amount of dependence within each cluster and requires it to diminish quite rapidly as $M\to\infty$. Although BCH_2011 discusses a particular model for which this requirement holds, it rules out simple DGPs such as the factor model \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}}, with or without cluster fixed effects. It even rules out the random-effects model \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}}, which is the most common model of intra-cluster correlation. For all these models, no CLT can possibly apply to the vector $M^{-1/2} {\bm{s}}_g$. Consider, for example, the factor model \hyperref[{eq:factor}]{\tagform@{\ref*{eq:factor}}}, and suppose that ${\bm{x}}_g=1$. In this case,
Because of the double summation, the second term on the right-hand side of \hyperref[{clt:variance}]{\tagform@{\ref*{clt:variance}}} clearly does not converge as $M \to \infty$ unless additional, and very strong, assumptions are made.
Another, quite different, approach to inference when $G$ is fixed is developed in IM_2010. The parameter of interest is a scalar, say $\beta$, which can be thought of as one element of ${\bm\beta}$. The key idea is to estimate $\beta$ separately for each of the $G$ clusters. This yields estimators $\hat\beta_g$ for $g=1,\ldots,G$. Inference is then based on the average, say $\bar\beta$, and standard error, say $s_{\kern -0.08333em\hat\beta}$, of the $\hat\beta_g$. IM_2010 shows that the test statistic $\sqrt{G}(\bar\beta-\beta_0) /s_{\kern -0.08333em\hat\beta}$ is approximately distributed as $t(G-1)$ when all clusters are large and a CLT applies to $N_g^{-1/2} {\bm{s}}_g$ for each $g$. However, as we saw above, this assumption cannot hold even for the simple random-effects or factor models.
A practical problem with this procedure is that $\beta$ may not be estimable for at least some clusters. For models of treatment effects at the cluster level, this will actually be the case for every cluster. For difference\kern 0.04167em-in-differences (DiD) models with clustering at the jurisdiction level, it will be the case for every jurisdiction that is never treated. IM_2016 suggests a way to surmount this problem by combining clusters into larger ones that allow $\beta$ to be estimated for each of them. Even when $\beta$ itself can be estimated for each cluster, however, the full model may not be estimable. This can happen, for example, when there are fixed effects for categorical variables, and not all categories occur in each cluster. In such cases, the interpretation of $\beta$ may differ across clusters.
Although the estimators and test statistics proposed in IM_2010,IM_2016 differ from the more conventional ones studied in BCH_2011, both approaches lead to $t$-statistics that follow the $t(G-1)$ distribution asymptotically. This distribution has in fact been used in Stata as the default for CV$_{\kern -0.08333em1}$-based inference for many years. For small values of $G$, using it can lead to noticeably more accurate, and more conservative, inferences than using the $t(N-k)$ or normal distributions. However, as we discuss in the next subsection and in (ref), inferences based on the $t(G-1)$ distribution are often not nearly conservative enough, especially when $G$ is small.
Whenever we rely on asymptotic theory, we need to be careful. What is true for infinitely large samples may or may not provide a good approximation for any actual sample. Unless the very strong assumptions discussed in (ref) are satisfied, we cannot expect to obtain reliable inferences when $G$ is small.
Unfortunately, there is no magic number for $G$ above which asymptotic inference can be relied upon. It is sometimes claimed that asymptotic inference based on CV$_{\kern -0.08333em1}$ is reliable when $G \geq 50$, or even (partly in jest) when $G \geq 42$ MHE_2008, but this is not true. In very favorable cases, inference based on CV$_{\kern -0.08333em1}$ and the $t(G-1)$ distribution can be fairly reliable when $G=20$, but in unfavorable ones it can be unreliable even when $G=200$ or more; see (ref). Moreover, there is evidence BM_2002,MNW-bootknife that inference based on CV$_{\kern -0.08333em3}$ tends to be more reliable, sometimes much more reliable, than inference based on other CRVEs.
What determines whether a case is favorable or unfavorable for a given $G$ is mostly the heterogeneity of the cluster score vectors. Unfortunately, the latter cannot be observed directly. We can observe the empirical score vectors, but they can sometimes differ greatly from the true ones ((ref)). We can also observe cluster sizes, which the discussion in (ref) focused on. These are often particularly important, but any form of heterogeneity can have serious consequences. This includes both heteroskedasticity of the disturbances at the cluster level and systematic variation across clusters in the distribution of the regressors. In general, the number of clusters $G$ and the extent to which the distribution of the scores varies across clusters will determine the quality of the asymptotic approximation.
In principle, a poor asymptotic approximation could lead $t$-tests based on the $t(G-1)$ distribution either to under-reject or over-reject. We have never observed $t$-tests based on CV$_{\kern -0.08333em1}$ or CV$_{\kern -0.08333em2}$ to under-reject in any simulation experiments, but we have observed ones based on CV$_{\kern -0.08333em3}$ to do so. Unless $G$ is fairly large, it is not difficult to find cases in which $t$-tests based on any of these CRVEs reject much more than their nominal level. Wald tests of several restrictions typically perform even worse; see (ref).
As we discussed in (ref), the condition \hyperref[{RateCondition}]{\tagform@{\ref*{RateCondition}}} imposes a restriction on the size of the largest cluster relative to the sample size. Thus, the quality of the asymptotic approximation will surely diminish as the size of the largest cluster increases relative to the average cluster size, and over-rejection will consequently increase. This conjecture is supported by simulation evidence in MW-JAE and DMN_2019, as well as by analytic results based on Edgeworth expansions in the latter paper.
There are at least two situations in which cluster-robust $t$-tests and Wald tests are at risk of over-rejecting to an extreme extent, namely, when one or a few clusters are unusually large, or when only a few clusters are treated. In both of these cases, one cluster, or just a few of them, have high leverage, in the sense that omitting one of these clusters has the potential to change the OLS estimates substantially; see (ref). Since both of these situations can occur even when $G$ is not small, all users of cluster-robust inference need to be on guard for them.
The first case in which conventional inference fails is when one or a very few clusters are much larger than the others. This implies that the distributions of the score vectors for those clusters are much more spread out than the ones for the rest of the clusters. An extreme example is studied in DMN_2019. When half the sample is in one large cluster, rejection rates for $t$-tests based on CV$_{\kern -0.08333em1}$ actually increase as $G$ increases, approaching 50% at the 5% level for $G=201$. Unfortunately, this extreme case is empirically relevant. Because roughly half of all incorporations in the United States are in Delaware, empirical studies of state laws and corporate governance encounter precisely this situation whenever they cluster at the state level HS_2020. MNW-bootknife studies a case with cluster sizes proportional to incorporations, where almost 53% of the observations are in the largest cluster. Although $t$-tests based on all three CRVEs over-reject, those based on CV$_{\kern -0.08333em3}$ do so much less severely than ones based on CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em2}$.
Not all forms of heterogeneity are harmful. In particular, having some extremely small clusters in a sample generally does not cause any problems, so long as there is not too much heterogeneity in the remainder of the sample. For example, suppose that a sample consists of, say, 25 large clusters, each with roughly 200 observations, and 15 tiny clusters, each with just one or a handful of observations. Except in very unusual cases, the coefficient estimates and their $t$-statistics would hardly change if we were to drop the tiny clusters, so this sample is better thought of as having 25 equal-sized clusters. The asymptotic approximations would perform just about the same whether or not the tiny clusters were included. Of course, if we changed this example so that there were 5 large clusters and 15 tiny ones, then asymptotic inference would surely be very problematic, because there would effectively be just 5 clusters.
The second case in which conventional inference fails is when the regressor of interest is a treatment dummy, and treatment occurs only for observations in a small number of clusters. In such cases, the empirical score vectors for the treated clusters, even when they have been modified, can provide very poor estimates of the actual score vectors.
Suppose that $d_{gi}$ is the value of the treatment dummy for observation $i$ in cluster $g$, and let $s^d_g$ denote the element of ${\bm{s}}_g$ corresponding to the dummy. Consider first the extreme case in which only some or all of the observations in the first cluster are treated. Then $s^d_g = \sum_{i=1}^{N_g} d_{gi} u_{gi}$ is equal to $\sum_{i=1}^{N_1} d_{1i} u_{1i}$ for $g=1$ and to 0 for all $g \neq 1$. Thus the scores corresponding to the treatment dummy equal zero for the control clusters. Moreover, because the treatment regressor must be orthogonal to the residuals, the empirical score $\hat s^d_1 = 0$. Since the actual score $s^d_1 \neq 0$, this implies that \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}} provides a dreadful estimate of \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}}, at least for the elements corresponding to the coefficient on the treatment dummy. In consequence, the CV$_{\kern -0.08333em1}$ standard error of this coefficient can easily be too small by a factor of five or more. When more than one cluster is treated, the problem is not as severe, because the $\hat s^d_g$ now sum to zero over the observations in all the treated clusters. This causes them to be too small, but not to the same extent as when just one cluster is treated; see MW-JAE,MW-EJ.
How well the empirical scores mimic the actual scores depends on the sizes of the treated and control clusters, the values of other regressors, and the number of treated observations within the treated clusters. Thus all these things affect the accuracy of cluster-robust standard errors and the extent to which $t$-statistics based on them over-reject. As the number of treated clusters, say $G_1$, increases, the problem often goes away fairly rapidly. But increasing $G$ when $G_1$ is small and fixed does not help and may cause over-rejection to increase. For models where all observations in each cluster are either treated or not, having very few control clusters is just as bad as having very few treated clusters. The situation is more complicated for DiD models, however; see MW-TPM.
It seems plausible that the modified empirical score vectors used by CV$_{\kern -0.08333em2}$ and CV$_{\kern -0.08333em3}$ may mimic the actual score vectors more accurately than the unmodified ones used by CV$_{\kern -0.08333em1}$. In fact, there is some evidence MNW-bootknife that $t$-tests based on CV$_{\kern -0.08333em2}$ over-reject less severely than ones based on CV$_{\kern -0.08333em1}$ when there are few treated clusters, and that $t$-tests based on CV$_{\kern -0.08333em3}$ over-reject considerably less severely. However, even though CV$_{\kern -0.08333em3}$ can perform much better than CV$_{\kern -0.08333em1}$, it still tends to over-reject when $G_1$ is very small.
Most of the literature on cluster-robust inference has focused on $t$-tests, but the cluster-robust Wald tests defined in \hyperref[{Waldstat}]{\tagform@{\ref*{Waldstat}}} also generally over-reject in finite samples. In fact, they tend to do so more severely as $r$, the number of restrictions, increases; see PT_2018 and JGM-fast. As is well known, this phenomenon occurs for Wald tests of all kinds. The problem may well be unusually severe in this case, however, because all CRVEs have rank at most $G$ and, in many cases, only $G-1$. On the other hand, although the true variance matrix $\operatorname{Var}(\hat{\bm\beta})$ is unknown in finite samples, it will normally have full rank $k$. As $r$ increases, the inverse of ${\bm{R}}^\top\hat{\bm{V}}{\bm{R}}$ thus seems likely to provide an increasingly poor approximation to the inverse of ${\bm{R}}^\top\!\operatorname{Var}(\hat{\bm\beta}){\bm{R}}$\kern 0.04167em.
BCH_2011 proves a result that helps to mitigate this problem. Under the very strong assumptions discussed in (ref), where $G$ is fixed, every cluster is the same size $M$, and the amount of within-cluster dependence is limited, the paper shows that $r(G-1)/(G-r)$ times the appropriate quantile of the $F(r,G-r)$ distribution provides an asymptotic critical value for the Wald statistic \hyperref[{Waldstat}]{\tagform@{\ref*{Waldstat}}} based on CV$_{\kern -0.08333em1}$. Equivalently, $(G-r)/(r(G-1))$ times $W$ is shown to follow the $F(r,G-r)$ distribution asymptotically in this special case. However, the WCR bootstrap ((ref)) can provide a much better approximation.
Although cluster-robust inference is most commonly used with the linear regression model \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}, it can actually be employed for a wide variety of models estimated by maximum likelihood or the generalized method of moments (GMM); see HansenLee_2019.
Consider a model characterized by the log-likelihood function
where ${\bm{\theta}}$ is the $k\times 1$ parameter vector to be estimated, and $\ell_{gi}({\bm{\theta}})$ denotes the contribution to the log-likelihood made by the \th{i} observation within the \th{g} cluster. Let $\hat{\bm{\theta}}$ denote the vector that maximizes \hyperref[{logl}]{\tagform@{\ref*{logl}}}, ${\bm{s}}_{gi}({\bm{\theta}})$ the $k\times 1$ vector of the first derivatives of $\ell_{gi}({\bm{\theta}})$ (that is, the score vector), and ${\bm{H}}_{gi}({\bm{\theta}})$ the $k\times k$ Hessian matrix of the second derivatives. Further, let $\hat{\bm{s}}_g = \sum_{i=1}^{N_g} {\bm{s}}_{gi}(\hat{\bm{\theta}})$ and $\hat{\bm{H}} = \sum_{g=1}^G \sum_{i=1}^{N_g} {\bm{H}}_{gi}(\hat{\bm{\theta}})$. Then HansenLee_2019 shows (using somewhat different notation) that the cluster-robust variance estimator for the maximum likelihood estimator $\hat{\bm{\theta}}$ is
The resemblance between \hyperref[{mlvar}]{\tagform@{\ref*{mlvar}}} and the CV$_{\kern -0.08333em1}$ variance matrix in \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}} is striking. Indeed, since the Hessian is proportional to ${\bm{X}}^\top{\bm{X}}$ for the linear regression model, CV$_{\kern -0.08333em1}$ without the leading scalar factor is really just a special case of \hyperref[{mlvar}]{\tagform@{\ref*{mlvar}}}.
The variance matrix estimator \hyperref[{mlvar}]{\tagform@{\ref*{mlvar}}} can be used for a wide variety of models estimated by maximum likelihood. In fact, Stata has been using it for various models, including probit and logit, for some years. HansenLee_2019 provides a similar result for GMM estimation, which is also very widely applicable. More recently, the fixed-$G$ approach discussed in (ref) has been applied to GMM estimation by Hwang_2021. It leads to a novel inferential procedure that involves modifying the usual asymptotic $t$ and $F$ statistics, but it requires that cluster sizes be approximately equal.
Unfortunately, very little currently seems to be known about the finite\kern 0.04167em-sample properties of tests based on \hyperref[{mlvar}]{\tagform@{\ref*{mlvar}}} or its GMM analog. They are probably worse than those of tests based on \hyperref[{Waldstat}]{\tagform@{\ref*{Waldstat}}}. It seems quite plausible that bootstrapping, to which we turn in the next section, would help, and Hwang_2021 discusses one bootstrap method for GMM estimation. However, bootstrapping nonlinear models tends to be computationally expensive, and the properties of applicable bootstrap procedures are largely unknown at the present time.
Instead of basing inference on an asymptotic approximation to the distribution of a statistic of interest, it is often more reliable to base it on a bootstrap approximation. In (ref), we briefly review some key concepts of bootstrap testing and bootstrap confidence intervals. Then, in (ref), we discuss bootstrap methods for regression models with clustered data. These methods, in particular the wild cluster restricted (WCR) bootstrap to be discussed in (ref), can be surprisingly inexpensive to compute and often lead to much more reliable inferences than asymptotic procedures. We therefore recommend that at least one variant of the WCR bootstrap be used almost all the time.
Suppose we are interested in a test statistic $\tau$, which might be a $t$-statistic or a Wald statistic. Instead of using $P$ values or critical values taken from an asymptotic distribution, we can use ones from the empirical distribution function (EDF) of a large number of bootstrap test statistics. This EDF often provides a good approximation to the unknown distribution of $\tau$. In order to obtain the EDF, we need to generate $B$ bootstrap samples and use each of them to compute a bootstrap test statistic, say $\tau^*_b$, for $b=1,\ldots,B$.
Precisely how the bootstrap samples are generated is critical, and we will discuss some methods for doing so in the next two subsections. The choice of $B$ also matters. Ideally, it should be reasonably large DM_2000 and satisfy the condition that $\alpha(B+1)$ is an integer for any $\alpha$ (the level of the test) that may be of interest RM_2007. In general, the computational cost of generating the bootstrap test statistics is proportional to $B$ times $N$\kern -0.08333em, so that bootstrapping can be expensive. However, as we discuss in (ref), surprisingly inexpensive methods are available for linear regression models with clustered disturbances. Unless computational cost is an issue, $B=9,\kern -0.08333em999$ and even $B=99,\kern -0.08333em999$ are generally good choices.
The EDF of the $\tau^*_b$ often provides a better approximation to $F(\tau)$, the distribution of $\tau$, than does its asymptotic distribution. This can sometimes be shown formally, but generally only under strong assumptions and at the cost of a great deal of algebra DMN_2019. For the model \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}, however, the intuition is quite simple. In many cases, the poor finite\kern 0.04167em-sample properties of test statistics based on CV$_{\kern -0.08333em1}$ arise because $\sum_{g=1}^G \hat{\bm{s}}_g\hat{\bm{s}}_g^\top$, provides a poor approximation to $\sum_{g=1}^G{\bm{\Sigma}}_g$. Often the bootstrap analog of the former provides a similarly poor approximation to the bootstrap analog of the latter. If so, then it is plausible that the empirical distribution of the $\tau^*_b$ will differ from the asymptotic distribution of the $\tau^*_b$ in roughly the same way as the distribution of $\tau$ differs from its asymptotic distribution. In that case, the EDF of the bootstrap test statistics should provide a reasonably good approximation to $F(\tau)$.
The EDF of the $\tau^*_b$ may be obtained by sorting the $\tau^*_b$ from smallest to largest. Number $(1-\alpha)(B+1)$ then provides an estimate of the $1-\alpha$ quantile, which may be used as the critical value for an upper-tail test at level $\alpha$. Identical inferences will be obtained by calculating the upper-tail bootstrap $P$ value,
and rejecting the null hypothesis whenever $\hat{P}^*(\tau) < \alpha$. Here $\tau$ could be either the Wald statistic \hyperref[{Waldstat}]{\tagform@{\ref*{Waldstat}}} or the absolute value of the $t$-statistic \hyperref[{eq:tstat}]{\tagform@{\ref*{eq:tstat}}}.
Setting $\tau = |t_a|$ in \hyperref[{bootp}]{\tagform@{\ref*{bootp}}} imposes symmetry on the bootstrap distribution of $t_a$. In many cases, it makes sense to do this, because cluster-robust $t$-statistics for linear regression models with exogenous regressors are often symmetrically distributed around the origin, at least to a good approximation. When $\tau = |t_a|$, the quantity $\hat{P}^* (\tau)$ defined in \hyperref[{bootp}]{\tagform@{\ref*{bootp}}} is a symmetric bootstrap $P$ value for a two\kern 0.04167em-sided test of ${\bm{a}}^\top({\bm\beta}-{\bm\beta}_0)=0$; see \hyperref[{eq:tstat}]{\tagform@{\ref*{eq:tstat}}}.
In dynamic models, nonlinear models, and models estimated by instrumental variables, however, it is common for the coefficients of interest to be biased. This causes the associated $t$-statistics to have non-zero means in finite samples. In such cases, it makes sense to use the equal-tail bootstrap $P$ value,
Here we compute upper-tail and lower-tail $P$ values, take the minimum of them, and then multiply by 2 to ensure that the nominal level of the test is correct.
There are many ways to construct a bootstrap confidence interval for a regression coefficient $\beta$ of which we have an estimate $\hat\beta$. A method that is conceptually (but not always computationally) simple is to invert a bootstrap test. This means finding two values of the coefficient, say $\beta_l$ and $\beta_u$, with $\beta_u > \beta_l$ and normally on opposite sides of $\hat\beta$, such that
Here $t(\beta=\beta_c)$, for $c=l$ and $c=u$, is a cluster-robust $t$-statistic for the hypothesis that $\beta=\beta_c$. The desired $1-\alpha$ confidence interval is then $[\beta_l,\, \beta_u]$. When the bootstrap DGP imposes the null hypothesis, the distribution of the bootstrap samples depends on the value of $\beta_c$. Solving the two equations in \hyperref[{eq:confint}]{\tagform@{\ref*{eq:confint}}} therefore requires iteration; see Hansen_1999 and JGM_2015,JGM-fast. In general, these methods tend to be expensive, but this is not the case for the WCR bootstrap to be discussed in (ref).
For bootstrap confidence intervals, it is common to use bootstrap DGPs that do not impose the null hypothesis, because no iteration is then required. The simplest method is just to calculate the standard deviation of the $\hat\beta^*_b$ and use this number, say $s^*(\hat\beta)$, as an estimator of the standard error of $\hat\beta$. The confidence interval is then
where $c_{1-\alpha/2}$ is the $1 - \alpha/2$ quantile of (in this case) the $t(G-1)$ distribution. A better approach, at least in theory, is to use the studentized bootstrap, or percentile\kern 0.04167em-$t$, confidence interval advocated in Hall_1992, which is
where $s_\beta$ is the standard error of $\hat\beta$ from the CRVE, and $c^*_z$ denotes the $z$ quantile of the bootstrap $t$-statistics $\tau_b^*$. Although the higher-order theory in DMN_2019 does not explicitly deal with confidence intervals, it strongly suggests that the intervals \hyperref[{eq:sstarint}]{\tagform@{\ref*{eq:sstarint}}} and \hyperref[{eq:studboot}]{\tagform@{\ref*{eq:studboot}}} should not perform as well as inverting a bootstrap test based on a bootstrap DGP that imposes the null hypothesis. Simulation results in JGM_2015 are consistent with these predictions. However, the intervals \hyperref[{eq:sstarint}]{\tagform@{\ref*{eq:sstarint}}} and \hyperref[{eq:studboot}]{\tagform@{\ref*{eq:studboot}}} have the advantage that they are easy to compute. No iteration is required, and a single set of bootstrap samples can be used to compute confidence intervals for all the parameters of interest.
The most important aspect of any bootstrap procedure is how the bootstrap samples are generated. The only procedure applicable to every model that uses clustered data is the pairs cluster bootstrap, which is also sometimes referred to as the cluster bootstrap, the block bootstrap, or resampling by cluster. The pairs cluster bootstrap works by grouping the data for every cluster into a $[{\bm{y}}_g,{\bm{X}}_g]$ pair and then resampling from the $G$ pairs. Every bootstrap sample is constructed by choosing $G$ pairs at random with equal probability $1/G$.
Although this procedure ensures that every bootstrap sample contains $G$ clusters, the number of observations inevitably varies across the bootstrap samples, unless all cluster sizes are the same. The size of the bootstrap samples can vary greatly, because the largest clusters may be over-represented in some bootstrap samples and under-represented in others. This limits the ability of the bootstrap samples to mimic the actual sample. So does the fact that the ${\bm{X}}^\top{\bm{X}}$ matrix is different for every bootstrap sample.
Because the pairs cluster bootstrap does not impose the null hypothesis, care must be taken when calculating the bootstrap test statistics. If the null hypothesis is that $\beta=\beta_0$, the actual $t$-statistic will have numerator $\hat\beta-\beta_0$, but the bootstrap $t$-statistic must have numerator $\hat\beta^*_b-\hat\beta$, where $\hat\beta^*_b$ is the estimator of $\beta$ for the \th{b} bootstrap sample. In this case, $\hat\beta$ is the parameter value associated with the bootstrap DGP. Because the bootstrap DGP does not impose the null hypothesis, the pairs cluster bootstrap cannot be used to construct the confidence interval \hyperref[{eq:confint}]{\tagform@{\ref*{eq:confint}}}, but it can be used to construct the intervals \hyperref[{eq:sstarint}]{\tagform@{\ref*{eq:sstarint}}} and \hyperref[{eq:studboot}]{\tagform@{\ref*{eq:studboot}}}.
CM_2015 discusses several problems that can arise with the pairs cluster bootstrap and sensibly suggests that investigators should examine the empirical distributions of the bootstrap coefficient estimates and test statistics. For example, if the bootstrap distribution has more than one mode, then it probably does not provide a good approximation to the actual distribution. This can happen when one or two clusters are very different from all the others; see (ref).
FP_2019 proposes several bootstrap procedures that can be thought of as variants of the pairs (not pairs cluster) bootstrap. The first step is to run regressions at either the individual level or the group $\times$ time\kern 0.04167em-period level, then aggregate the residuals so that there is just one residual per cluster, and finally run bootstrap regressions on the resampled residuals. These procedures include parametric methods to correct for the heteroskedasticity generated by variation in the number of observations per group. Remarkably, they can work well even with just one treated cluster. However, this is possible only because, unlike methods based on a CRVE, they do not allow for unrestricted heteroskedasticity.
In general, the pairs cluster bootstrap is expensive to compute. However, a computational shortcut for linear regression models can make it feasible even when $N$ and $B$ are both large; see JGM-fast. Nevertheless, we do not recommend this method for the linear regression model \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}, because, as we discuss in the next subsection, a much better method is available. With nonlinear models such as the probit model, the pairs cluster bootstrap may be attractive even though it can be expensive. However, we cannot recommend it without reservation, because it appears that very little is known about its finite\kern 0.04167em-sample properties. Simulation evidence suggests that, at least in some cases, it can either over-reject or under-reject severely MW-TPM,JGM-fast.
The restricted wild cluster, or WCR, bootstrap was first suggested in CGM_2008. Until recently, the only variant of it ever used in practice is based on CV$_{\kern -0.08333em1}$, and, at time of writing, this classic variant is the only one for which software is readily available. Other variants will be discussed briefly at the end of this subsection.
Suppose that $\tilde{\bm\beta}$ denotes the OLS estimator of ${\bm\beta}$ subject to the restriction ${\bm{a}}^\top{\bm\beta}={\bm{a}}^\top{\bm\beta}_0$, which is to be tested, and $\tilde{\bm{u}}_g = {\bm{y}}_g - {\bm{X}}_g\tilde{\bm\beta}$ denotes the vector of restricted residuals for the \th{g} cluster. Then the traditional way to write the WCR bootstrap DGP is
where the $v_g^{\ast b}$ are independent realizations of an auxiliary random variable $v^*$ with zero mean and unit variance. In practice, the best choice for $v^*$ is usually the Rademacher distribution, in which case $v^*$ equals 1 or $-1$ with equal probabilities DF_2008,DMN_2019. This imposes symmetry on the bootstrap disturbances.
Instead of using \hyperref[{eq:wcr}]{\tagform@{\ref*{eq:wcr}}}, we can generate the bootstrap scores directly as
where $\tilde{\bm{s}}_g = {\bm{X}}_g^\top \tilde{\bm{u}}_g$ is the score vector for the \th{g} cluster evaluated at the restricted estimators. Plugging the ${\bm{s}}_g^{*b}$ into \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}} then yields the bootstrap estimators $\hat{\bm\beta}^*_b$. This method is computationally inexpensive JGM-fast. It also provides an intuitive justification for the WCR bootstrap because, as we stressed in (ref), the finite\kern 0.04167em-sample properties of cluster-robust tests depend mainly on the properties of the scores.
DMN_2019 establishes the asymptotic validity of the classic WCR bootstrap with CV$_{\kern -0.08333em1}$ standard errors and also studies the unrestricted wild cluster, or WCU, bootstrap. The difference between the WCR and WCU bootstraps is that, for the latter, the unrestricted scores $\hat{\bm{s}}_g$ are used instead of the restricted ones $\tilde{\bm{s}}_g$ in the bootstrap DGP \hyperref[{eq:scoreboot}]{\tagform@{\ref*{eq:scoreboot}}}. In general, the WCU bootstrap does not perform as well in finite samples as the WCR one. The paper contains both theoretical results (based on Edgeworth expansions) and simulation results to support this assertion. However, the WCU bootstrap has the advantage that the bootstrap DGP does not depend on the restrictions to be tested. The same set of bootstrap samples can therefore be used to perform tests on any restriction or set of restrictions and/or to construct confidence intervals based on either \hyperref[{eq:sstarint}]{\tagform@{\ref*{eq:sstarint}}} or \hyperref[{eq:studboot}]{\tagform@{\ref*{eq:studboot}}} for any coefficient of interest.
DMN_2019 also proves the asymptotic validity of two variants of the ordinary wild bootstrap, restricted (WR) and unrestricted (WU), for the model \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}. The ordinary wild bootstrap uses $N$ realizations of $v^*$\kern -0.08333em, one for each observation, instead of just $G$. This means that the disturbances for the bootstrap samples are uncorrelated within clusters. Although this implies that the distribution of the $\hat{\bm\beta}^*_b$ cannot possibly match that of $\hat{\bm\beta}$, it often does not prevent the distribution of the $\tau^*_b$ from providing a good approximation to the distribution of $\tau$. With one important exception (see below), the WR bootstrap seems to work less well than the WCR bootstrap, as asymptotic theory predicts. It can also be much more expensive to compute when $N$ is large.
The classic versions of the WCR and WCU bootstraps that use CV$_{\kern -0.08333em1}$ are surprisingly inexpensive to compute. The computations require only the matrices ${\bm{X}}_g^\top{\bm{X}}_g$ and the vectors ${\bm{X}}_g^\top{\bm{y}}_g$; see JGM-fast. The calculations can be made even faster by rewriting the bootstrap test statistic so that it depends on all the sample data in the same way for every bootstrap sample; see \citet*{RMNW}. The only thing that varies across the bootstrap samples is the $G$-vector ${\bm{v}}^{*b}$ of realizations of the auxiliary random variable. There are some initial computations that may be expensive when $N$ is large, but they only have to be done once. After that, the ${\bm{v}}^{*b}$ and the results of the initial computations are used to compute all the bootstrap test statistics.
This fast procedure is implemented in the package boottest for both Stata and Julia; for details, see RMNW. In (ref), there is an illustration of how fast it can be; see the notes to (ref). Importantly, boottest not only computes WCR bootstrap $P$ values for both $t$-tests and Wald tests; it also computes WCR bootstrap confidence intervals based on \hyperref[{eq:confint}]{\tagform@{\ref*{eq:confint}}}. The package has many other capabilities as well.
In many cases (exceptions will be discussed below), the classic CV$_{\kern -0.08333em1}$ variant of the WCR bootstrap yields very accurate inferences. These are generally more accurate than those for the pairs cluster or WCU bootstraps. In addition to DMN_2019, see CGM_2008, MW-JAE, and JGM-fast. Because the ${\bm{X}}_g$ are always the same, the wild cluster bootstrap is able to replicate what is often the main source of heterogeneity, namely, variation in cluster sizes, in every bootstrap sample. This is not the case for the pairs cluster bootstrap, and it surely contributes to the superior accuracy of inferences based on the WCR bootstrap.
Of course, no bootstrap method can work perfectly. Not surprisingly, the performance of the WCR bootstrap using CV$_{\kern -0.08333em1}$ tends to deteriorate as $G$ becomes smaller, as the number of regressors increases, and as the clusters become more heterogeneous. In particular, it can sometimes perform very badly when the number of treated clusters $G_1$ is very small MW-JAE, MW-TPM, MW-EJ. This is true both for pure treatment models, where all the observations in each cluster are either treated or not, and for DiD models, where only some observations in the treated clusters are treated.
Unlike other methods, which generally over-reject severely when there are few treated clusters ((ref)), the WCR bootstrap usually under-rejects in this case. This happens because the distribution of the bootstrap statistics $\tau^*_b$ depends on the value of the actual test statistic $\tau$. The larger is $\tau$, the more dispersed are the $\tau^*_b$. This makes $\hat P^*(\tau)$ in \hyperref[{bootp}]{\tagform@{\ref*{bootp}}} larger than it should be. In the most extreme case of just one treated cluster, rejection frequencies may be essentially zero. In this case, the bootstrap distribution is often bimodal MW-JAE, so that plotting it can provide a useful diagnostic. When there are few treated clusters and the WCR bootstrap $P$ value seems suspiciously large, it may be worth trying the ordinary wild restricted (WR) bootstrap, which can sometimes work surprisingly well in this context MW-EJ. This is a case where randomization inference can be attractive; see (ref).
We recommend using at least one variant of the WCR bootstrap (preferably with at least $B=9,\kern -0.08333em999$) almost all the time. All variants are often remarkably inexpensive to compute, and they often seem to work well. When $G$ is not too small and the clusters are not too heterogeneous, WCR bootstrap $P$ values and confidence intervals may be quite similar to ones based on CV$_{\kern -0.08333em1}$ or CV$_{\kern -0.08333em3}$ with $t(G-1)$ critical values. In that case, it is likely that finite\kern 0.04167em-sample issues are not severe, and there is probably no need to do anything else. When there is a large discrepancy between the results of different methods, however, it would make sense to try several bootstrap methods and maybe some of the methods discussed in (ref).
The classic version of the WCR bootstrap can sometimes work remarkably well even when $G$ is very small. In fact, \citet*{CSS_2021} shows that it can yield exact inferences in certain cases where $N$ is large and $G$ is small. These results are obtained by exploiting the relationship between the WCR bootstrap with Rademacher auxiliary random variables and randomization inference ((ref)). However, they require rather strong homogeneity conditions on the distribution of the covariates across clusters, as well as limits on the amount of dependence within each cluster similar to those in BCH_2011.
At this point, a word of warning is in order. Almost all the simulation results that we have referred to are based on models with one or just a few regressors, and these regressors are typically generated in a fairly simple way. Moreover, almost all existing simulations focus on $t$-statistics rather than Wald statistics. There is evidence that rejection frequencies for all methods increase as the number of regressors increases, and that Wald tests are less reliable than $t$-tests; see DMN_2019 and MNW-bootknife on the former point, PT_2018 on the latter, and JGM-fast on both of them. Thus the classic WCR bootstrap may perform less well in empirical applications with large numbers of regressors than it has typically done in simulations, especially when there is more than one restriction.
When $G$ is small, the WCR bootstrap encounters an important practical problem. For the Rademacher distribution, or any other two\kern 0.04167em-point distribution, the number of possible bootstrap samples is just $2^G$\kern -0.08333em. Webb_2014 proposes a six-point distribution which largely solves this problem, because $6^G$ is much larger than $2^G$\kern -0.08333em. This distribution seems to work almost as well as Rademacher for most values of $G$, and sometimes much better when $G$ is very small. Whenever either $2^G$ (for Rademacher) or $6^G$ (for six-point) is smaller than the chosen value of $B$, we can enumerate all possible bootstrap samples instead of drawing them at random. For example, when $G=16$, there are just $65,\kern -0.08333em536$ Rademacher bootstrap samples to enumerate (one of which is identical to the actual sample). This eliminates simulation randomness from the bootstrap procedure. In fact, boottest uses enumeration by default whenever $B$ is greater than the number of possible bootstrap samples.
Most of the discussion above pertains to the classic variant of the WCR bootstrap, which uses \hyperref[{eq:scoreboot}]{\tagform@{\ref*{eq:scoreboot}}} to generate bootstrap samples and CV$_{\kern -0.08333em1}$ to calculate standard errors. MNW-bootknife introduces three new variants. The simplest of these still uses \hyperref[{eq:scoreboot}]{\tagform@{\ref*{eq:scoreboot}}}, but the actual and bootstrap test statistics employ CV$_{\kern -0.08333em3}$ standard errors, computed using \hyperref[{eq:jackvar}]{\tagform@{\ref*{eq:jackvar}}}. Two somewhat more complicated variants generate the bootstrap scores using an alternative to \hyperref[{eq:scoreboot}]{\tagform@{\ref*{eq:scoreboot}}} that implicitly involves a transformation of the restricted scores similar to the one that yields the $\acute{\bm{s}}_g$ in \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}}. One of these new variants uses CV$_{\kern -0.08333em1}$ standard errors, and the other uses CV$_{\kern -0.08333em3}$ standard errors. Simulation results in MNW-bootknife suggest that all three new variants of the WCR bootstrap outperform the classic variant in a number of circumstances. In particular, the tendency to over-reject seems to increase much more slowly for the new variants than for the classic one as the number of regressors increases and as the variation in cluster sizes increases. However, the new variants still tend to under-reject severely when the number of treated clusters is small. We hope that Stata and/or \texttt{R} packages able to compute these new variants will become available in the near future.
Bootstrap methods are not the only way to obtain inferences more accurate than those given by cluster-robust $t$-tests and confidence intervals using the $t(G-1)$ distribution. Numerous other procedures have been proposed, which broadly fall into two categories that we discuss in the following two subsections. Except in the case of treatment models with very few clusters or very few treated clusters ((ref)), we do not recommend any of these procedures in preference to CV$_{\kern -0.08333em3}$ combined with the $t(G-1)$ distribution or the WCR bootstrap, especially the newest variants of the latter. However, when the recommended methods yield conflicting results, it is surely a good idea to try other methods as well.
Critical values for cluster-robust test statistics can be based on various approximations. The first paper to take this approach seems to be BM_2002. It suggests methods for calculating approximations to CV$_{\kern -0.08333em2}$ or CV$_{\kern -0.08333em3}$ $t$-statistics based on the Student's $t$ distribution with an estimated degrees-of-freedom (d-o\kern 0.04167em-f) parameter. These employ what is called a “Satter\-thwaite approximation” to calculate the d-o\kern 0.04167em-f. This is done under the assumption that the variance matrix of ${\bm{u}}$ is proportional to an identity matrix. The d-o\kern 0.04167em-f parameter is different for every hypothesis to be tested, and it can be much less than $G-1$.
Imbens_2016 proposes a similar procedure for $t$-tests based on CV$_{\kern -0.08333em2}$ under the assumption that the variance matrix of ${\bm{u}}$ corresponds to a random-effects model. As we saw in (ref), such a model implies that the disturbances within each cluster are equi-correlated, and the intra-cluster correlation $\rho$ must be estimated from the residuals. When cluster fixed effects are partialed out, doing so absorbs any random effects, making this approach inapplicable.
AY-exact proposes a related method that uses CV$_{\kern -0.08333em1}$ instead of CV$_{\kern -0.08333em2}$. It involves two steps. In the first step, the CV$_{\kern -0.08333em1}$ standard error for the coefficient of interest is multiplied by a factor greater than one. In the second step, a d-o\kern 0.04167em-f parameter is calculated. The standard error and the d-o\kern 0.04167em-f parameter can then be used to compute either a $t$-statistic and its $P$ value or a confidence interval. A Stata package called edfreg is available.
The three procedures just discussed are described in some detail in MW-EJ. However, the the CV$_2$-based and CV$_3$-based procedures can be very expensive as described there, and NAAMW_2020 provides a better way to compute them. Limited simulation evidence suggests that the performance of Young's method is similar to those of the two methods based on CV$_{\kern -0.08333em2}$. However, this evidence is by no means definitive, because the simulations focus on only a narrow set of treatment models.
Using Hotelling's $T^{\kern 0.08333em2}$ distribution with estimated degrees of freedom, PT_2018 generalizes the CV$_{\kern -0.08333em2}$ procedure of BM_2002 to the case of Wald tests based on \hyperref[{Waldstat}]{\tagform@{\ref*{Waldstat}}}. Simulations suggest that the resulting tests always reject less often than standard Wald tests. They rarely over-reject but often under-reject, and they sometimes do so quite severely. The clubSandwich package for R and the reg$\_$sandwich package for Stata implement the procedures discussed in PT_2018.
Although the procedures discussed in this subsection have some theoretical appeal and seem to work well in many cases, we are not aware of any evidence that they outperform either tests based directly on CV$_{\kern -0.08333em3}$ or the wild cluster bootstrap for a range of models and DGPs. One limitation is that the Student's $t$ distribution is not very flexible. Even when the d-o\kern 0.04167em-f parameter is estimated very well, the best we can hope for is that a test based on this distribution will be accurate for some level of interest. It may well over-reject at some levels and under-reject at others.
Our recommendation is to use the methods discussed in this subsection primarily to confirm (or perhaps cast doubt on) the results of the methods that we recommend when there is concern about the reliability of the latter. Cases of particular concern are ones with few but balanced clusters (say, $G \leq 12$), ones with balanced but few treated (or few control) clusters (say $G_1 \leq 6$ or $G-G_1 \leq 6$), ones with seriously unbalanced cluster sizes (even when $G$ is quite large), ones with treated clusters that are unusually large or small, and ones with any sort of heterogeneity that causes a few clusters to have high leverage; see (ref). It is always reassuring when several methods yield essentially the same inferences.
Randomization inference (RI) was proposed by Fisher_1935 as a distribution-free way to perform hypothesis tests in the context of experiments. RI tests are also called permutation tests. Lehmann_2005 gives a formal introduction, and Imbens_2015 provides a more accessible discussion. The key idea of RI is to compare an outcome $\tau$ that is actually observed with a set of outcomes that might have been observed if treatment had been assigned differently. The outcome could be a sample average, a coefficient estimate, or a test statistic.
Specifically, consider a clustered regression model with treatment at the cluster level, which could be a DiD model where only some observations within the treated clusters receive treatment. Then $\tau$ might be the average treatment effect for some outcome measure. RI procedures may be attractive when it is plausible that treatment was assigned at random and the researcher is interested in the sharp null hypothesis that the treatment has no effect. They can sometimes yield reliable results even when the number of clusters is small and/or the number of treated clusters is small.
Suppose there are $G$ clusters, $G_1$ of which received treatment. The number of ways in which treatment could have been assigned to $G_1$ out of the $G$ clusters is
One of these corresponds to the actual assignment, and the remaining $S = {}_{G}C_{G_1} - 1$ are called re\kern 0.04167em-randomizations. Each re\kern 0.04167em-randomization involves pretending that a particular set of $G_1$ clusters was treated, with the remaining $G-G_1$ serving as controls. The values of the dependent variable do not change across re\kern 0.04167em-randomizations, only the values of the treatment dummy. For every re\kern 0.04167em-randomization, indexed by $j$, we could calculate a test statistic $\tau^*_j$. If the observable characteristics of the clusters were all the same, it would make sense to compare $\tau$ with the empirical distribution of the $\tau^*_j$. To do so, it is natural to calculate the $P$ value for an upper-tail test as either
Here $P_2^*$ implicitly includes the actual assignment to treatment, and $P_1^*$ omits it.
When $\alpha(S+1)$ is an integer for a test at level $\alpha$, both $P$ values in \hyperref[{Pvalues}]{\tagform@{\ref*{Pvalues}}} yield the same result, and the test is exact if the distributions of the $\tau^*_j$ are the same as that of $\tau$. However, $P_1^*$ and $P_2^*$ can differ noticeably when $S$ is small. The latter is more conservative and seems to be more popular. For moderate values of $S$, it is often easy to enumerate all of the possible $\tau^*_j$, but this is infeasible even when $G$ is not particularly large. In such cases, we must choose a number of re\kern 0.04167em-randomizations, say $B=99,\kern -0.08333em999$, at random. In principle, these should be drawn without replacement, but that is not important if $B$ is small relative to $S$.
It seems to be widely believed that tests based on RI are always exact. This is not true. When treatment is not assigned at random, or when the observed characteristics of the treated clusters differ systematically from those of the controls, we cannot expect the distributions of $\tau$ and the $\tau^*_j$ to coincide.
For DiD models, CT_2011 proposes a test that is very similar to an RI test based on the OLS estimate of the coefficient on a treatment dummy. That paper also shows how to obtain a confidence interval by inverting the test. MW-RI studies two procedures for DiD tests, one based on coefficient estimates (called RI-$\beta$) and one based on cluster-robust $t$-statistics (called RI-$t$). Both procedures work very well when the clusters are homogeneous and $G$ is reasonably large, even when $G_1=1$. However, when the treated clusters are systematically larger or smaller than average, neither RI-$\beta$ nor RI-$t$ tests perform well, although the latter typically perform better. In such cases, it appears that $G$ may have to be quite large (much larger than for the WCR bootstrap) before either procedure works really well, even when $G_1$ is not particularly small.
AY_fisher contains an interesting recent application of RI, where the RI-$\beta$ and RI-$t$ procedures are applied to regressions for the results of 53 randomized experiments in published papers. Stata packages called randcmd and randcmdci are available. In many cases, estimates that were originally reported to be significant are not significant according to the RI procedures. Not surprisingly, this is particularly true for regressions with a few high-leverage clusters or observations.
There is evidently a close relationship between RI and the wild cluster bootstrap. Evaluating all possible bootstrap samples by enumeration is quite similar to evaluating all possible re\kern 0.04167em-randomizations. The results of CSS_2021 and MW-RI strongly suggest that homogeneity across clusters is more important for RI than for the WCR bootstrap. For the former, the regressors change as we compute each of the $\tau^*_j$, but the regressand stays the same. For the latter, the regressand changes as we compute each of the $\tau^*_b$, but the regressors stay the same. Thus, the empirical distribution to which $\tau$ is being compared is conditional on the actual regressors (including the clusters actually treated) for the WCR bootstrap, but not for the RI procedures discussed above.
The RI methods discussed so far are not the only ones. \citet*{CRS_2017} proposes “approximate randomization tests” based on the cluster-level estimators of IM_2010; see \citet*{CCKS_2021} for a guide to this approach. In models for treatment effects, every cluster must include both treated and untreated observations. This can be accomplished by merging clusters, but at the cost of making $G$ smaller with a resulting loss in power. Hag_2019a,Hag_2019b develops RI tests that can be used even when $G$ is quite small and there is substantial heterogeneity across clusters. These tests do not require cluster-level estimation, but $G_1$ and $G-G_1$ should both be no less than 4. HS_2020 develops an RI procedure based on RI-$t$ for the case in which one cluster is much larger than any of the others.
An alternative way to deal with the problem of one or very few treated clusters, which can involve an RI-like procedure, is the method of synthetic controls surveyed in Abadie_2021.
Many studies fail to report enough information to convince readers that their empirical results should be believed. Investigators often simply report test statistics or confidence intervals based on CV$_{\kern -0.08333em1}$ together with the $t(G-1)$ distribution. Although results of this type may be reliable, they often will not be. Unless the number of clusters is fairly large and evidence demonstrates a convincing degree of homogeneity across clusters, results from at least one or two other methods should be reported. These might include tests and/or confidence intervals based on CV$_{\kern -0.08333em3}$, $P$ values and/or confidence intervals based on one or more variants of the WCR bootstrap discussed in (ref), or perhaps results from some of the alternative methods discussed in (ref).
It is important to report some key information about the sample as a matter of routine. The fundamental unit of inference when the observations are clustered is not the observation but the cluster; this is evident from \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}} and \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}}. The asymptotic theory discussed in (ref) therefore depends on $G$, not $N$\kern -0.08333em. With the exception of certain very special cases discussed in (ref), asymptotic approximations tend to work poorly when there are few clusters. It is therefore absolutely essential to report the number of clusters, $G$, whenever inference is based on a CRVE. This is even more important than reporting $N$\kern -0.08333em.
Moreover, because the distributional approximations perform best when the scores are homogeneous across clusters, and the most important source of heterogeneity is often variation in cluster sizes, it is extremely important to report measures of this variation. At a minimum, we believe that investigators should always report the median cluster size and the maximum cluster size, in addition to $N$ and $G$. When $G$ is small, or when the distribution of the $N_g$ is unusual, it would be good to report the entire distribution of cluster sizes in the form of either a histogram or a table.
Of course, clusters can be heterogeneous in many ways beyond their sizes. In (ref), we discussed a measure of leverage at the cluster level. Another measure is the partial leverage of each cluster for the coefficient(s) of interest, which generalizes the observation-level notion of partial leverage of CW_1980. If $\acute{\bm{x}}_j$ denotes the vector of residuals from regressing ${\bm{x}}_j$ on all the other regressors, and $\acute{\bm{x}}_{gj}$ is the subvector of $\acute{\bm{x}}_j$ corresponding to the \th{g} cluster, then the partial leverage of cluster $g$ for the \th{j} coefficient is
When a cluster has high partial leverage for the \th{j} coefficient, removing that cluster has the potential to change the estimate of the \th{j} coefficient substantially.
A popular way to quantify the heterogeneity of clusters is to calculate $G_j^*$\kern -0.08333em, the “effective number of clusters” for the \th{j} coefficient, as proposed in \citet*{CSS_2017}. This number is always less than $G$, and it can provide a useful warning when $G_j^*$ is much smaller than $G$. The value of $G_j^*$ depends on an unknown parameter $\rho$, the intra-cluster correlation of the disturbances when they are assumed to be equi-correlated. CSS_2017 suggests setting $\rho=1$, as a sort of worst case. However, since cluster fixed effects absorb all intra-cluster correlation for the random effects model \hyperref[{eq:REmodel}]{\tagform@{\ref*{eq:REmodel}}}, only $\rho=0$ makes sense for models with such fixed effects, and trying to use $\rho\ne0$ can lead to numerical instabilities. MNW-influence discusses these issues and shows how to calculate $G_j^*(\rho)$ efficiently for any value of $\rho$ in models without cluster fixed effects.
Using \hyperref[{partlev}]{\tagform@{\ref*{partlev}}} and the definition of $G_j^*(0)$, it can be shown that the latter is a monotonically decreasing function of the squared coefficient of variation of the partial leverages $L_{gj}$, which is defined as
where $\bar L_j$ is the sample mean of the $L_{gj}$. Thus $G_j^*(0)$ and $V_s(L_{j\bullet})$ provide exactly the same information. When the partial leverage for the \th{j} coefficient is the same for every cluster, $V_s(L_{j\bullet}) = 0$ and $G_j^*(0) = G$. When the $L_{gj}$ differ greatly, $V_s(L_{j\bullet})$ is large, and $G_j^*(0) <\kern -0.08333em< G$.
We suggest that investigators should routinely report the leverages $L_g$ defined in \hyperref[{eq:traceXX}]{\tagform@{\ref*{eq:traceXX}}}, the partial leverages $L_{gj}$ for the coefficient(s) of interest, and the $\hat\beta_j^{(g)}$ for the same coefficient(s), at least when they provide information beyond that in the distribution of the cluster sizes. When $G$ is small, it may be feasible to report all these numbers. Otherwise, it may be feasible to graph them, as in (ref), or to report summary statistics such as $G_j^*(0)$ or $V_s(L_{j\bullet})$. Reporting measures of influence and leverage should help to identify cases in which inference may be unreliable, as well as sometimes turning up interesting, or possibly disturbing, features of the data.
When there is clustering in two or more dimensions ((ref)), we recommend that investigators compute everything just suggested for each of the clustering dimensions and also for their intersection(s). If the results are interesting, they should be reported. For two\kern 0.04167em-way clustering, this might mean reporting, or at least summarizing, up to three sets of leverages, partial leverages, and $\hat\beta_j^{(g)}$\kern -0.08333em.
We illustrate many of our recommendations by revisiting a long-standing empirical question in labor economics, namely, the impact of the minimum wage on young people. In the past two decades, many U.S.\ states have significantly increased their minimum wages. In fact, from 2000 to 2019, every state increased the nominal minimum wage by at least 27%, and six states doubled it. Moreover, recent proposals to increase the national minimum wage to \$15 per hour have reinvigorated the debate about the effects of minimum wages.
Some classic references on the impact of the minimum wage are Mincer_1976 and CK_1994. The latter paper was among the very first and most influential applications of the DiD methodology, which continues to be used in this literature. For example, Seattle_2017 uses a DiD analysis to study the effects of a large increase in the minimum wage in Seattle. Wolfson_2019 and Neumark_2021 survey many recent studies on the impacts of minimum wages. Both conclude that the majority of studies, but not all, find dis-employment effects that are concentrated among teenagers and those with low levels of education. Manning_2021 explores why the evidence on employment effects is mixed.
Instead of using a DiD approach, we exploit state\kern 0.04167em-level differences in the minimum wage and analyze their impacts on labor-market and education outcomes at the individual level. Although we treat the minimum wage as exogenous in our analysis, we hesitate to call our estimates causal. There is reason to believe that state\kern 0.04167em-level minimum wages may be endogenous, because states may be more likely to increase them during good economic times. Moreover, it is possible that the effects of the minimum wage differ depending on the state of the economy. However, we ignore these issues, because our principal objective is to illustrate the importance of clustering for statistical inference.
The model we estimate is
Here $y_{ist}$ is the outcome of interest for person $i$ in state $s$ in year $t$. There are three outcome variables. “Hours” records the usual hours worked per week, which is defined only for employed individuals. “Employed” is a binary variable equal to 1 if person $i$ is employed and to 0 if they are either unemployed or not in the labor force. “Student” is a binary variable equal to 1 if person $i$ is currently enrolled in school and to 0 otherwise. The parameter of interest is $\beta$, which is the coefficient on $\textrm{mw}_{\kern -0.08333em st}$, the minimum wage in state $s$ at time $t$. The row vector ${\bm{Z}}_{ist}$ collects a large set of individual-level controls, including race, gender, age, and education. There are also year and state fixed effects. NW_2007 estimates models similar to \hyperref[{reg:wage}]{\tagform@{\ref*{reg:wage}}} with individual-level data, clustering at the state level.
Data at the individual level from the American Community Survey (ACS) were obtained from IPUMS IPUMS_2020 and cover the years 2005--2019. The minimum wage data were provided by Neumark_2019, and we have collapsed them to state\kern 0.04167em-year averages to match the ACS frequency. Following previous literature, we restrict attention to teenagers aged 16--19. We keep only individuals who are “children” of the respondent to the survey and who have never been married. We drop individuals who had completed one year of college by age 16 and those reporting in excess of 60 hours usually worked per week. We also restrict attention to individuals who identify as either black or white.
We consider six different clustering structures that lead to nine different standard errors for $\hat\beta$. These are no clustering (with HC$_1$ standard errors), one\kern 0.04167em-way clustering at either the state\kern 0.04167em-year, state, or region\footnote{Regions are defined as the U.S.\ Census Divisions, with the following partitioning of states. New England: CT, MA, MN, NH, RI, VT; Middle Atlantic: NJ, NY, PA; South Atlantic: DC, DE, FL, GA, MD, NC, SC, VA, WV; East South Central: AL, KY, MS, TN; East North Central: IL, IN, MI, OH, WI; West North Central: IA, KS, MN, MO, ND, NE, SD; West South Central: AR, LA, OK, TX; Mountain: AZ, CO, ID, MT, NM, NV, UT, WY; Pacific: AK, CA, HI, OR, WA.} level (with both CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em3}$ standard errors), and two\kern 0.04167em-way clustering by state and year or by region and year (with CV$_{\kern -0.08333em1}$ standard errors, the only ones currently available for multi-way clustering). Early empirical research on the impacts of the minimum wage would have used either conventional or HC$_1$ standard errors, but modern studies would almost always cluster at some level.
Since the minimum wage is invariant within state\kern 0.04167em-year clusters, it seems highly likely that the scores are correlated within them, and we therefore consider (at least) state\kern 0.04167em-year clustering. However, after a state has increased its minimum wage, it almost always remains at the new level until it is increased again. This implies that minimum wages must be correlated across years within each state. Unless the disturbances happen to be uncorrelated across years within states, the scores will therefore be correlated, which suggests that state\kern 0.04167em-level clustering may be appropriate. We consider region-level clustering based on the nine census divisions because there may be correlations among nearby states. Largely for completeness, we also consider two\kern 0.04167em-way clustering by either state or region and year.
(ref) presents the estimates of $\beta$ from regression \hyperref[{reg:wage}]{\tagform@{\ref*{reg:wage}}} for all three regressands, along with $t$-statistics and $P$ values for each of the clustering levels considered. The coefficients are negative for the hours and employment regressands and positive for the student one. Under the surely false assumption of no clustering, all of these coefficients are extremely significant, with $P$ values below $10^{-6}$. We do not compute bootstrap $P$ values, because it would be prohibitively costly, and they must all be very close to zero.
When we instead cluster at the state\kern 0.04167em-year level, the $t$-statistics become smaller, especially for the employment model. Nevertheless, clustering at this level still leads us to conclude that all three coefficients are significant, with WCR bootstrap $P$ values ranging from 0.0005 to 0.0145. However, some of our conclusions change radically when we cluster at the state or region levels. For all three outcome variables, the $t$-statistics become smaller, and the $P$ values (especially the CV$_{\kern -0.08333em3}$ and bootstrap ones) become larger. The bootstrap $P$ values often differ substantially from the ones based on the $t(G-1)$ distribution, which is expected given that the clusters are either very unbalanced or small in number ((ref)). Clustering by region always yields larger bootstrap $P$ values than clustering by state.
If instead we adopt the conservative approach of clustering at the level with the largest standard error, we cluster at the state level for the student model, but at the region level for the hours and employment models. Because the number of regions is so small, however, the latter results may not be reliable, even when we use CV$_{\kern -0.08333em3}$ or the WCR bootstrap. In our view, it is probably most reasonable to cluster at the state level in all three cases.
Luckily, the choice between state and region clustering does not matter much. In either case, we find from (ref) that an increase in the minimum wage is associated with a decrease in employment, but the coefficient is not significant at any conventional level using either state or region clustering. It is also associated with an increased probability of being a student, which is significant at the 5% level regardless of clustering level. For hours worked, the coefficient is negative. It is significant at the 5% level for state clustering, but not for region clustering.
The table also reports asymptotic and bootstrap $t$-statistics and $P$ values for two\kern 0.04167em-way clustering, either by state and year or by region and year. The bootstrap method we use combines the one\kern 0.04167em-way WCR bootstrap with the two\kern 0.04167em-way variance matrix \hyperref[{eq:betavar}]{\tagform@{\ref*{eq:betavar}}}. This can be done in two different ways, corresponding to each of the two clustering dimensions. Based on simulation evidence in MNW_2021, we bootstrap by the dimension with the smallest number of clusters. Since two\kern 0.04167em-way clustering does not change most of the results very much, however, there appears to be little reason to use it in this case.
We tentatively conclude that an increase in the minimum wage is associated with a significant decrease in hours worked and a significant increase in the likelihood of being a student. Seattle_2017 find a similar reduction in hours following a minimum wage increase in Seattle, and the effects of the minimum wage on school enrollment are discussed in NW-1995. For employment, we obtain a small negative effect, but it is not close to being statistically significant when we cluster at either the state or region levels. Thus our results are consistent with, and add support for, the conclusions of Manning_2021 about the “elusive” employment effect of the minimum wage.
To confirm the choice of clustering level, we use the score\kern 0.04167em-variance tests of MNW-testing to test for the appropriate level of clustering.\footnote{We do not calculate the tests proposed in IM_2016, because \hyperref[{reg:wage}]{\tagform@{\ref*{reg:wage}}} cannot be estimated on a cluster-by-cluster basis for state\kern 0.04167em-level or state\kern 0.04167em-year clusters without dropping some regressors.} (ref) presents results from six tests for each of the three models. Following a systematic, sequential testing approach, we would test independence vs.\ state\kern 0.04167em-year, then state\kern 0.04167em-year vs.\ state, and finally state vs.\ region. Apart from the possibility of two\kern 0.04167em-way clustering, we would conclude that the appropriate level of clustering is that of the first non-rejected null hypothesis. For the hours model, the score\kern 0.04167em-variance tests marginally favor region clustering over state clustering, but, for the employment and student models, they clearly favor state clustering. Thus, these tests lend support to our decision to cluster at the state level; see the discussion of (ref).
As we discussed in (ref), investigators should be suspicious of any results that are overly dependent on very few clusters. (ref) reports several summary statistics for cluster heterogeneity. Specifically, it reports $G$, the number of clusters and $G_\beta^*(0)$, the effective number of clusters ((ref)), as well as the average, minimum, first and third quartiles, median, and maximum of the $N_g$. These statistics suggest that the state\kern 0.04167em-year clusters are extremely unbalanced. Although there are $G=765$ clusters, the effective number $G_\beta^*(0)$ is smaller than $G$ by a factor of roughly 10 to 12. The maximum cluster size is over six times the median, and the third quartile is nearly twice the median.
From (ref), we also see that the state clusters are extremely unbalanced, based both on their sample sizes and on the values of $G_\beta^*(0)$. The region clusters are fairly balanced in terms of their sample sizes, with the maximum $N_g$ only about four times the minimum. The values of $G_\beta^*(0)$ are also not worrisome. The year clusters (which we only use in two\kern 0.04167em-way clustering) are very balanced in terms of their sample sizes but much less so in terms of $G^*(0)$. We would expect asymptotic tests based on CV$_{\kern -0.08333em1}$ to perform rather poorly in all cases, because there is a lot of cluster heterogeneity for clustering by state\kern 0.04167em-year and state, a small number of clusters for clustering by region, and a small number of effective clusters for clustering by year.
For each of the three levels of one\kern 0.08333em-way clustering, we calculate measures of leverage, partial leverage, and influence using the summclust package in Stata MNW-influence. However, we only report results for state\kern 0.04167em-level clustering. Since state\kern 0.04167em-year clustering is strongly rejected against state\kern 0.04167em-level clustering ((ref)), it does not seem interesting to discuss it. We also do not discuss region-level clustering, because there are just nine regions, and there is not much evidence in favor of clustering at that level.
The empirical distributions of the $N_g$, the $L_g$, and the $L_{g\beta}$ (the partial leverages) are shown in (ref), where the $N_g$ and the $L_g$, which actually sum to $N$ and $k$, respectively, have been rescaled so that they sum to unity. In a perfectly balanced case, where all of the $N_g$ (or $L_g$ or $L_{g\beta}$) are identical, these EDFs would simply be vertical lines at $1/G$. Both datasets are evidently far from being perfectly balanced. The $L_g$ are essentially proportional to the $N_g$ and extremely highly correlated with them; the correlations are 0.9986 for the employment/student dataset and 0.9982 for the hours dataset. This is evident in both panels, where the EDFs for the $N_g$ and the $L_g$ are extremely similar. Thus, for leverage, it appears that there is no heterogeneity in the clusters other than what is implied by the cluster sizes.
In contrast, the partial leverages vary much more across states than do the other two measures. The largest values of the $L_{g\beta}$ are considerably larger than the largest values of the $N_g$ and the $L_g$, especially for the employment/student dataset; note that all of these are for California. Not coincidentally, more of the $L_{g\beta}$ than of the other measures are smaller than $1/G$, the average. The correlations between the $L_{g\beta}$ and the $L_g$ are 0.9150 for the student/employment dataset and 0.9165 for the hours dataset. The correlations between the $L_{g\beta}$ and the $N_g$ are also very similar for the two datasets, at 0.9209 and 0.9163.
The substantial amount of cluster heterogeneity that is evident in both (ref) and (ref) suggests that inference based on CV$_{\kern -0.08333em1}$ and the $t(50)$ distribution may not be reliable. Further evidence on this point will be provided in (ref).
For the student model, we find no evidence of influential clusters. The $\hat \beta^{(g)}$ range from 0.00185 to 0.00246, with $\hat \beta = 0.00221$. On the other hand, for the employment model, there are two noticeable values of $\hat\beta^{(g)}$: Texas has $\hat\beta^{(g)} = -0.00529$ and New York has $\hat\beta^{(g)} = -0.00203$, with $\hat\beta = -0.00367$ and all the remaining $\hat\beta^{(g)}$ in the interval $[-0.00455 ; -0.00302]$. For comparison, the standard error of $\hat\beta$ is $0.00268$.
For the hours model, four values of the $\hat\beta^{(g)}$ stand out. They are Minnesota ($-0.1776$) and Texas ($-0.1770$) at one end and Arizona ($-0.1274$) and New York ($-0.1227$) at the other. For comparison, $\hat\beta = -0.1539$, and all the other $\hat\beta^{(g)}$ lie in the interval $[-0.1713; -0.1406]$. In this case, the standard error of $\hat\beta$ is $0.0623$.
As we discussed in (ref), placebo regressions provide a simple way to check the level at which the residuals are clustered, even when the pattern of intra-cluster correlation is unknown and perhaps very complicated. If a placebo regressor is clustered at, say, the state level, then the empirical scores will also be clustered at that level unless the residuals are clustered only at a finer level. When the residuals display intra-cluster correlation at the state level, we would therefore expect placebo\kern 0.04167em-regression tests with standard errors clustered at that level to reject roughly as often as they should. But we would expect placebo regressions with standard errors clustered at finer levels to over-reject.
We perform two sets of placebo\kern 0.04167em-regression experiments for each of the three equations estimated in (ref). In the first set, the placebo regressor is a DiD-style treatment dummy similar to the ones used in BDM_2004 and the other papers cited in (ref). Treatment is randomly applied to 15, 20, 25, 30, or 35 states. For each state, it begins randomly in any year excluding 2005 (to avoid collinearity with the state fixed effects) and continues through 2019.\footnote{Of course, this sort of DiD model with two\kern 0.04167em-way fixed effects is somewhat obsolete; see CS_2021 and other papers cited therein. But we would still expect to find no effects for placebo treatments beyond those attributable to chance.} Rejection percentages are shown in the top panel of (ref). The first number in each pair is the smallest rejection percentage observed over the five experiments for each equation, and the second number is the largest one. The numbers of treated and control states are deliberately chosen not to be extreme, so as to avoid the issues discussed in (ref).
For the second set of simulation experiments, the placebo regressor is generated by
where the $\epsilon_{ist}$ and the $e_{st}$ are independent standard normals. Thus the $v_{st}$ are 51 separate stationary AR(1) processes, and $x_{ist}$ is a weighted average of $v_{st}$ and $\epsilon_{ist}$. When either $\rho=0$ or $\delta=0$, the $x_{ist}$ are independent. When both $\rho$ and $\delta$ are positive, they are correlated within both state\kern 0.04167em-years and states. They are never correlated across states within regions. The extent to which the $x_{ist}$ are correlated within state\kern 0.04167em-years and states depends on both of the parameters in \hyperref[{xprocess}]{\tagform@{\ref*{xprocess}}}. For simplicity, we report rejection percentages for just two cases. In the first case, $\rho=0.5$ and $\delta=0.9$, so there is a lot of correlation within each state\kern 0.04167em-year. In the second case, $\rho=0.8$ and $\delta=0.5$, so there is less correlation within state\kern 0.04167em-years but more correlation across years within each state.
The first row in (ref) shows that failing to cluster always leads to severe over-rejection, and the next three rows show that clustering at the state\kern 0.04167em-year level also leads to substantial over-rejection. This is true for all methods of inference and for both DGPs, but more so for the DiD-type one. Rows 5 through 7 show that there can be noticeable over-rejection for clustering at the state level when using CV$_{\kern -0.08333em1}$ with $t(G-1)$ critical values, especially for the DiD-type DGP. This is to be expected given the unbalanced cluster sizes at the state level. Both CV$_{\kern -0.08333em3}$ with $t(G-1)$ critical values and the classic WCR bootstrap perform much better. For the AR(1) placebo regressors, both CV$_{\kern -0.08333em3}$ and the WCR bootstrap perform very well indeed when there is clustering at either the state or region level. For both designs, the performance of CV$_{\kern -0.08333em3}$ at the region level is remarkably good. However, the WCR bootstrap tends to outperform it slightly at the state level.
The placebo\kern 0.04167em-regression results in (ref) are largely consistent with those of the score\kern 0.04167em-variance tests. They suggest that the CV$_{\kern -0.08333em3}$ and classic WCR bootstrap results with state\kern 0.04167em-level clustering in (ref) can probably be relied upon, but that the results without clustering or with state\kern 0.04167em-year clustering should not be believed.
We conclude by presenting a brief summary guide. This is essentially a checklist for cluster-robust inference in regression models. The first three items should be dealt with as soon as possible in the process of specifying and estimating a model. The remaining ones should be kept in mind throughout the process of estimation and inference.
{1.5pt} \addcontentsline{toc}{section}{\refname}