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.
140,936 characters · 21 sections · 58 citation commands
Testing for the appropriate level of clustering in linear regression models
\onehalfspacing
Modern empirical econometrics often allows for correlation within clusters of observations, and this can have serious consequences for statistical inference. Theoretical work on cluster-robust inference almost always assumes that the structure of the clusters is known, even though the form of the correlation within clusters is arbitrary. Unless it is obvious that clustering must be at a certain level, however, this can leave empirical researchers in a difficult situation. They must generally rely on rules of thumb, their own intuition, or referees' suggestions to decide how the observations should be clustered. To make this process easier, we propose tests for any given level of clustering (including no clustering as a special case) against an alternative within which it is nested. When two or more levels of clustering are possible, we propose a sequence of such tests.
There has been a great deal of research on cluster-robust inference in the past two decades. \citet*{CM_2015} cover much of the literature up to a few years ago. Esarey_2019 and MW-survey provide more recent surveys. \citet*{CGH_2018} deal with a broader class of methods for various types of dependent data. \citet*{MNW-guide} provide a thorough and detailed guide to empirical practice. Areas that have received particular attention include: asymptotic theory for cluster-robust inference \citep*{DMN_2019,HansenLee_2019}; bootstrap methods with clustered data \citep*{CGM_2008, DMN_2019, RMNW, MNW-bootknife}; and inference with unbalanced clusters \citep*{Imbens_2016, CSS_2017, MW-JAE, DMN_2019, MNW-influence}.
Almost all of this literature assumes that the way in which observations are allocated to clusters is known to the econometrician. This is quite a strong assumption. Imagine that a dataset has many observations taken from individuals in different geographical locations. In order to utilize a cluster-robust variance estimator (CRVE), the researcher needs to specify at what level the clustering occurs. For example, there could possibly be clustering at the zip\kern 0.04167em-code, city, county, state, or country level. Even in this relatively simple setting, there are many possible ways in which a researcher could `cluster' the standard errors.
A few rules of thumb have emerged to cover some common cases. For instance, in the case of nested clusters, such as cities within states, \citet*{CM_2015} advocate clustering at the larger, more aggregate level. In the case of randomized experiments, \citet*{Athey_2017} recommend clustering at the level of randomization. In the case of experiments where treatment is assigned to groups in pairs, with one group treated and one not treated, Chaisemartin_2022 recommend clustering at the pair level rather than the group level. While these rules of thumb can sometimes be very helpful, they may or may not lead to the appropriate clustering level in any particular case.
Getting the level of clustering correct is extremely important. Simulation results in several papers have shown that ignoring clustering in a single dimension can result in rejection frequencies for tests at the 5% level that are actually well over 50% \citep*{BDM_2004, CGM_2008} and confidence intervals that are too narrow by a factor of five or more JGM-CJE. On the other hand, clustering at too coarse a level (say state\kern 0.04167em-level clustering when there is actually city-level clustering) can lead to the problems associated with having few treated clusters, which can be severe MW-JAE,MW-EJ, and can also reduce power MW-survey.
In (ref), we propose two tests for the cluster structure of the error variance matrix in a linear regression model. They test the null hypothesis of a fine level of clustering (or of no clustering at all) against an alternative hypothesis with a coarser level of clustering. The tests are based on the difference between two functions of the scores for the parameter(s) of interest. These functions are essentially the filling in the sandwich for two different cluster-robust variance estimators, one associated with the null level of clustering and one associated with the alternative level. Since the functions estimate the variance of the scores under two different clustering assumptions, we refer to the tests as score\kern 0.04167em-variance, or SV, tests. A procedure for sequential testing, described in (ref), allows for determination of the appropriate level of clustering without inflating the family-wise error rate when there are several possible levels of clustering.
Tests for the appropriate level of clustering have also been proposed by Ibragimov_2016 and recently by Cai_2022. These tests are very different from our tests and very different from each other. We discuss them briefly in (ref).
The model of interest is discussed in (ref). Our score\kern 0.04167em-variance tests are described in (ref), including the bootstrap implementation, the sequential testing procedure, and the use of our tests as pre\kern 0.04167em-tests for inference about regression coefficients. (ref) provides asymptotic theory for the two test statistics, the bootstrap tests, and the sequential testing procedure. In (ref), we consider the common situation in which the regressors that are not of primary interest have been partialed out prior to performing the test. The size and power of the proposed tests are analyzed by Monte Carlo simulations in (ref). An empirical example that deals with clustering by classroom or school using the STAR dataset Finn_1990,Mosteller_1995 is discussed in (ref). Finally, (ref) concludes and offers some guidance for empirical researchers. All mathematical proofs are given in (ref).
We focus on the linear regression model
where ${\bm{y}}$ and ${\bm{u}}$ are, respectively, $N \times 1$ vectors of observations and disturbances (or error terms), and ${\bm{X}}$ is an $N \times k$ matrix of regressors (or covariates). The $k \times 1$ parameter vector ${\bm\beta}$ contains the coefficients on the regressors.
Suppose that the data are divided into $G$ clusters, indexed by $g$, where the \th{g} cluster has $N_g$ observations, so that $N = \sum_{g=1}^G N_g$. Thus, there are $G$ vectors ${\bm{y}}_g$ and ${\bm{u}}_g$ of size $N_g$, along with $G$ matrices ${\bm{X}}_g$, each with $N_g$ rows and $k$ columns. Using this notation, the ordinary least squares (OLS) estimator of ${\bm\beta}$ is
where ${\bm\beta}_0$ denotes the true value of ${\bm\beta}$. Now define the $k \times 1$ score vectors ${\bm{s}}_g = {\bm{X}}_g^\top {\bm{u}}_g$. We assume that these score vectors satisfy ${\rm E}({\bm{s}}_g)={\bm{0}}$ for all $g$ and
where ${\mathbb I} (\cdot)$ denotes the indicator function and ${\bm{\Sigma}}_g$ is a $k \times k$ variance matrix. Although the properties of the ${\bm{\Sigma}}_g$ depend on the properties of the variance matrix of ${\bm{u}}$, we do not explicitly make any assumptions about the latter because our tests are concerned solely with the variances of the score vectors.
It is clear from \hyperref[{betahat}]{\tagform@{\ref*{betahat}}} that the asymptotic distribution of $\hat{\bm\beta}$ depends on the asymptotic distribution of the score vectors. An estimator of the variance matrix of $\hat{\bm\beta}$ is given by the sandwich formula
where $\hat{\bm{\Sigma}}$ is an estimator of the variance matrix of the sum of scores, ${\bm{\Sigma}} = {\rm E} ({\bm{X}}^\top{\bm{u}}{\bm{u}}^\top\!{\bm{X}} )$. The condition \hyperref[{def Sigma}]{\tagform@{\ref*{def Sigma}}} implies that ${\rm E} ({\bm{s}}_g {\bm{s}}_{g'}^\top) = {\bm{0}}$ whenever $g \neq g'$. In this case ${\bm{\Sigma}} =\sum_{g=1}^G {\bm{\Sigma}}_g$, so that the usual estimator for ${\bm{\Sigma}}$ under condition \hyperref[{def Sigma}]{\tagform@{\ref*{def Sigma}}} is
where $\hat{\bm{u}}_g$ contains the residuals for cluster $g$ and $\hat{\bm{s}}_g = {\bm{X}}_g^\top \hat{\bm{u}}_g$ is the $k \times 1$ vector of empirical scores for cluster $g$. The scalar factor $m_{\rm c}$ is a finite\kern 0.04167em-sample correction, the most commonly employed factor being $m_{\rm c}=G/(G-1)\times (N-1)/(N-k)$, which is designed to account for degrees of freedom. Using $\hat{\bm{\Sigma}} = \hat{\bm{\Sigma}}_{\rm c}$ in \hyperref[{covbeta}]{\tagform@{\ref*{covbeta}}} yields CV$_1$, the most widely-used CRVE for $\hat{\bm\beta}$. Asymptotic inference on regression coefficients using CV$_1$ is studied by DMN_2019 and HansenLee_2019.
The fundamental idea of our testing procedure is to compare two estimates of the variance of the coefficient(s) that we want to estimate. We test the null hypothesis that a CRVE based on a “fine” clustering structure is valid against the alternative that the CRVE needs to be based on a “coarser” clustering structure. Since it is only the filling in the sandwich \hyperref[{covbeta}]{\tagform@{\ref*{covbeta}}} that differs across different clustering structures, we are actually comparing two estimates of the variance matrix of the sum of scores. Our procedure is somewhat like the specification test of Hausman_1978. The “fine” CRVE is efficient when there actually is fine clustering, but it is invalid when there is coarse clustering. In contrast, the “coarse” CRVE is inefficient when there actually is fine clustering, but it is valid in both cases.
To make our testing procedure operational, we formulate it in terms of the parameters of the model. To this end, we first define some notation. There are $G$ coarse clusters indexed by $g=1,\ldots,G$. Within coarse cluster $g$, there are $M_g$ fine clusters indexed by $h=1,\ldots,M_g$. In total there are $G_{\rm f} = \sum_{g=1}^G M_g$ fine clusters. Fine cluster $h$ in coarse cluster $g$ contains $N_{gh}$ observations indexed by $i=1,\ldots,N_{gh}$. Coarse cluster $g$ therefore contains $N_g = \sum_{h=1}^{M_g} N_{gh}$ observations, and the entire sample contains $N=\sum_{g=1}^G N_g =\sum_{g=1}^G \sum_{h=1}^{M_g} N_{gh}$ observations. We let ${\bm{X}}_{ghi}$ and $u_{ghi}$ denote the regressors and disturbance for observation $i$ within fine cluster $h$ in coarse cluster $g$. We then define the corresponding score as ${\bm{s}}_{ghi}= {\bm{X}}_{ghi}^\top u_{ghi}$, the score for fine cluster $h$ in coarse cluster $g$ as ${\bm{s}}_{gh}=\sum_{i=1}^{N_{gh}} {\bm{s}}_{ghi}$, and the score for coarse cluster $g$ as ${\bm{s}}_g = \sum_{h=1}^{M_g}{\bm{s}}_{gh}$.
Under the coarse clustering structure, the ${\bm{s}}_g$ satisfy \hyperref[{def Sigma}]{\tagform@{\ref*{def Sigma}}}, so that in particular they are uncorrelated across $g$. Under the fine clustering structure, the ${\bm{s}}_{gh}$ are themselves uncorrelated across $h$, for each $g$. That is, for all $g=1,\ldots ,G$,
where each of the ${\bm{\Sigma}}_{gh}$ is a $k \times k$ matrix. Thus, \hyperref[{def Sigma}]{\tagform@{\ref*{def Sigma}}} and \hyperref[{def Sigma gh}]{\tagform@{\ref*{def Sigma gh}}} embody the assumption that the fine clustering structure is nested within the coarse one. Another possible design would have a one\kern 0.04167em-way clustering structure nested within a two\kern 0.04167em-way one; see (ref).
Now let ${\bm{\Sigma}}_{\rm c}$ and ${\bm{\Sigma}}_{\rm f}$ denote the matrix ${\bm{\Sigma}}$ under the coarse and fine clustering structures, respectively. From \hyperref[{def Sigma}]{\tagform@{\ref*{def Sigma}}} and \hyperref[{def Sigma gh}]{\tagform@{\ref*{def Sigma gh}}}, these matrices are
We consider the null and alternative hypotheses
The hypotheses are expressed in this way, rather than in terms of the difference between the limits of normalized versions of ${\bm{\Sigma}}_{\rm f}$ and ${\bm{\Sigma}}_{\rm c}$, because the appropriate normalizing factors will, in general, be unknown; see DMN_2019.
Our score\kern 0.04167em-variance, or SV, test statistics are based on comparing estimates $\hat{\bm{\Sigma}}_{\rm f}$ and $\hat{\bm{\Sigma}}_{\rm c}$ obtained under fine and coarse clustering, respectively. There are many ways in which one could compare these $k \times k$ matrices. We focus on two quantities of particular interest, which define two test statistics. The first is obtained for $k=1$. This could be after all regressors except one have been partialed out ((ref)), so that interest is focused on a particular coefficient that we are trying to make inferences about. This leads to a test statistic with the form of a $t$-statistic. The second is obtained for $k > 1$, in which case our test statistic is a quadratic form involving all the unique elements of $\hat{\bm{\Sigma}}_{\rm f}$ and $\hat{\bm{\Sigma}}_{\rm c}$, as in \citepos{White_1980} “direct test” for heteroskedasticity. The first test is of course a special case of the second, but we treat it separately because it is particularly simple to compute and may often be of primary interest.
In order to derive the test statistics, we write $\hat{\bm{\Sigma}}_{\rm c}$ and $\hat{\bm{\Sigma}}_{\rm f}$ using common notation. Let $\hat{\bm{s}}_{ghi}$ denote the empirical score for observation $i$ within fine cluster $h$ in coarse cluster $g$, and let $\hat{\bm{s}}_{gh} = \sum_{i=1}^{N_{gh}} \hat{\bm{s}}_{ghi}$ denote the empirical score for fine cluster $h$ in coarse cluster $g$, such that $\hat{\bm{s}}_g = \sum_{h=1}^{M_g} \hat{\bm{s}}_{gh}$. Under coarse clustering, the estimated ${\bm{\Sigma}}_{\rm c}$ matrix in \hyperref[{Sighat}]{\tagform@{\ref*{Sighat}}} is
Similarly, we can write, c.f.\ \hyperref[{def Sigma gh}]{\tagform@{\ref*{def Sigma gh}}} and \hyperref[{Sigmas}]{\tagform@{\ref*{Sigmas}}},
where $m_{\rm f} = G_{\rm f}/(G_{\rm f}-1) \times (N-1)/(N-k)$.
When interest focuses on just one coefficient, so that $k=1$, the matrix ${\bm{X}}$ becomes the vector ${\bm{x}}$, and the empirical scores are scalars. Specifically, $\hat s_{ghi} = x_{ghi} \hat u_{ghi}$ and $\hat s_{gh}=\sum_{i=1}^{N_{gh}} \hat s_{ghi}$ denote the empirical scores for observation $i$ and fine cluster $h$, respectively. Then the matrices \hyperref[{Sigmac}]{\tagform@{\ref*{Sigmac}}} and \hyperref[{Sigmaf}]{\tagform@{\ref*{Sigmaf}}} reduce to the scalars
The quantities given in \hyperref[{Sigmac}]{\tagform@{\ref*{Sigmac}}}, \hyperref[{Sigmaf}]{\tagform@{\ref*{Sigmaf}}}, and \hyperref[{Sigmascalar}]{\tagform@{\ref*{Sigmascalar}}} are all defined in essentially the same way. They simply amount to different choices of empirical scores. If $N_{gh}=1$, then $\hat\sigma^2_{\rm f}$ simplifies to
which is just the sum of the squared empirical scores over all the observations.
Our first test is based on the difference between the two scalars in \hyperref[{Sigmascalar}]{\tagform@{\ref*{Sigmascalar}}}, namely,
Our second test is based on the difference between the $k\times k$ matrices $\hat{\bm{\Sigma}}_{\rm c}$ and $\hat{\bm{\Sigma}}_{\rm f}$. For this test, we consider the vector of contrasts,
where the operator $\operatorname{vech} (\cdot)$ returns a vector, of dimension $k(k+1)/2$ in this case, with all the supra-diagonal elements of the symmetric $k\times k$ matrix argument removed.
In order to obtain test statistics with asymptotic distributions that are free of nuisance parameters, we need to derive the asymptotic means and variances of $\hat\theta$ and $\hat{\bm{\theta}}$, so that we can studentize the statistics in \hyperref[{thetascalar}]{\tagform@{\ref*{thetascalar}}} and \hyperref[{thetaSigma}]{\tagform@{\ref*{thetaSigma}}}. To this end, suppose that we observe the (scalar) scores $s_{gh}$ for fine cluster $h$ in coarse cluster $g$. Then the analog of $\hat\theta$ is the contrast
This is simply the sum of all the cross\kern 0.04167em-products of scores that are in the same coarse cluster but different fine clusters. Under the null hypothesis, $\theta$ clearly has mean zero by \hyperref[{def Sigma gh}]{\tagform@{\ref*{def Sigma gh}}}.
The variance of $\theta$ in \hyperref[{thetaknown}]{\tagform@{\ref*{thetaknown}}} is, under the null hypothesis,
The expectation of any product of scores can only be nonzero, under the null and \hyperref[{def Sigma gh}]{\tagform@{\ref*{def Sigma gh}}}, when their indices are the same in pairs. This implies that either $h_1=\ell_1 \neq h_2=\ell_2$ or $h_1=\ell_2 \neq h_2=\ell_1$. These cases are symmetric, and hence \hyperref[{var1}]{\tagform@{\ref*{var1}}} simplifies to
where $\sigma^2_{gh} = \operatorname{Var} ( s_{gh})$ is used to denote ${\bm{\Sigma}}_{gh}$ in the scalar case.
The sample analog of the right-hand side of \hyperref[{tauvar}]{\tagform@{\ref*{tauvar}}} is $2 \sum_{g=1}^G\sum_{h_1=1}^{M_g} \sum_{h_2\neq h_1}^{M_g}\hat s_{gh_1}^2 \hat s_{gh_2}^2$; see \hyperref[{def Sigma gh}]{\tagform@{\ref*{def Sigma gh}}}. This suggests the variance estimator
This equation avoids the triple summation in \hyperref[{tauvar}]{\tagform@{\ref*{tauvar}}} by squaring the sums of squared empirical scores, which then requires that the second term be subtracted. In deriving \hyperref[{varfast}]{\tagform@{\ref*{varfast}}}, we have ignored the factors $m_{\rm c}$ and $m_{\rm f}$, which are asymptotically irrelevant. If instead we had retained them, there would be no cancellation when subtracting $\hat\sigma^2_{\rm f}$ from $\hat\sigma_{\rm c}$, leading to a much more complicated (and computationally burdensome) expression for $\widehat\operatorname{Var} (\hat\theta)$. Combining \hyperref[{thetascalar}]{\tagform@{\ref*{thetascalar}}} and \hyperref[{varfast}]{\tagform@{\ref*{varfast}}} yields the studentized test statistic
In (ref), we show that $\tau_\sigma$ is asymptotically distributed as ${\rm N}(0,1)$.
When $k>1$, so that $\hat{\bm{\theta}}$ is a vector, the variance estimator analogous to \hyperref[{varfast}]{\tagform@{\ref*{varfast}}} is
Here ${\bm{H}}_k$ is the so\kern 0.04167em-called elimination matrix, which satisfies $\operatorname{vech} (\bm{S})={\bm{H}}_k {\textrm{vec}} (\bm{S})$ for any $k \times k$ symmetric matrix $\bm{S}$ Harville_1997, and $\otimes$ denotes the Kronecker product. A studentized (Wald) statistic is then given by
In (ref), we show that $\tau_\Sigma$ is asymptotically distributed as $\chi^2(k(k+1)/2)$.
In this section, we have proposed two score\kern 0.04167em-variance tests of \hyperref[{hypotheses}]{\tagform@{\ref*{hypotheses}}}. They both involve comparing different variance estimates of the empirical scores, namely, the two scalars in \hyperref[{Sigmascalar}]{\tagform@{\ref*{Sigmascalar}}} for the $\tau_\sigma$ test and the matrices in \hyperref[{Sigmac}]{\tagform@{\ref*{Sigmac}}} and \hyperref[{Sigmaf}]{\tagform@{\ref*{Sigmaf}}} for the $\tau_\Sigma$ test. The former is a special case of the latter, and it can be obtained for the same models as the latter by partialing out all regressors except one, as in (ref). This special case is particularly interesting, because the $\tau_\sigma$ test can be directional, and also because many equations simplify neatly in the scalar case.
As we show in (ref), the finite\kern 0.04167em-sample properties of our asymptotic tests are often good but could sometimes be better, especially when the number of clusters under the alternative is quite small. In such cases, we therefore recommend the use of bootstrap tests based on the statistics \hyperref[{eq:taus}]{\tagform@{\ref*{eq:taus}}} and \hyperref[{eq:tauGs}]{\tagform@{\ref*{eq:tauGs}}}, which often perform much better in finite samples, as we also show in (ref). These bootstrap implementations are described next.
The simplest way to implement a bootstrap test based on any of our test statistics is to compute a bootstrap $P$ value, say $\hat P^*$\kern -.08333em, and reject the null hypothesis when it is less than the level of the test. The bootstrap methods that we propose are based on either the ordinary wild bootstrap Wu_1986, Liu_1988 or the wild cluster bootstrap CGM_2008. These bootstrap methods are normally used to test hypotheses about ${\bm\beta}$, and we are not aware of any previous work in which they have been used to test hypotheses about the variances of parameter estimates. The asymptotic validity of the bootstrap tests that we now describe is established in (ref).
The key idea of the wild bootstrap is to obtain the bootstrap disturbances by multiplying the residuals by realizations of an auxiliary random variable with mean 0 and variance 1. In contrast to many applications of the wild bootstrap, the residuals in this case are unrestricted, meaning that they do not impose a null hypothesis on ${\bm\beta}$. This is because we are not testing any restrictions on ${\bm\beta}$ when testing the level of clustering. In the special case of testing the null of heteroskedasticity, as in (ref), we use the ordinary wild bootstrap. When the null involves clustering, we use the wild cluster bootstrap. Because the test statistics depend only on residuals, the value of ${\bm\beta}$ in the bootstrap DGP does not matter, and so we set it to zero.
The \th{b} wild (cluster) bootstrap sample is thus generated by ${\bm{y}}^{*b} = {\bm{u}}^{*b}$, where the vector of bootstrap disturbances ${\bm{u}}^{*b}$ has typical element given by either $u_{ghi}^{*b} = v_{ghi}^{*b} \hat u_{ghi}$ for the wild bootstrap or $u_{ghi}^{*b} = v_{gh}^{*b} \hat u_{ghi}$ for the wild cluster bootstrap. The auxiliary random variables $v_{ghi}^{*b}$ and $v_{gh}^{*b}$ are assumed to follow the Rademacher distribution, which takes the values $+1$ and $-1$ with equal probabilities. Notice that there is one such random variable per observation for the wild bootstrap and one per cluster for the wild cluster bootstrap. Other distributions can also be used; see DF_2008, DMN_2019, and Webb_2022.
The algorithm for a wild (cluster) bootstrap\kern 0.04167em-based implementation of our tests is as follows. It applies to both $\tau_\sigma$ and $\tau_\Sigma$. For simplicity, the algorithm below simply refers to one test statistic, $\tau$. However, it is easy to perform two or more tests at the same time, using just one set of bootstrap samples for all of them. For example, if there are three possible regressors of interest, we might perform four tests, one with $k=3$ based on $\tau_\Sigma$ and three with $k=1$ based on different versions of $\tau_\sigma$.
As usual, if $\alpha$ is the level of the test, then $B$ should be chosen so that $(1-\alpha)(B+1)$ is an integer RM_2007. Numbers like 999 and 9,999 are commonly used because they satisfy this condition for conventional values of $\alpha$. Power increases in $B$, but it does so very slowly once $B$ exceeds a few hundred DM_2000.
In many applications, there are several possible levels of clustering. In such situations, we suggest a sequential testing procedure. The statistical principle upon which we base our testing procedure is the intersection-union (IU) principle BergerSinclair_1984, whereby a hypothesis is rejected if and only if the hypothesis itself, along with any hypotheses nested within it, are all rejected. The IU principle leads naturally to a bottom-up testing strategy for the level of clustering. BergerSinclair_1984 show that the IU principle does not imply an inflation of the family-wise rejection rate in the context of multiple testing; that is, there is no accumulation of size due to testing multiple hypotheses. We prove a similar result for our sequential procedure below.
Suppose the potential levels of clustering are sequentially nested, and denote their ${\bm{\Sigma}}$ matrices by ${\bm{\Sigma}}_0, {\bm{\Sigma}}_1,\ldots,{\bm{\Sigma}}_p$; see \hyperref[{covbeta}]{\tagform@{\ref*{covbeta}}} and \hyperref[{Sigmas}]{\tagform@{\ref*{Sigmas}}}. Here we assume that ${\bm{\Sigma}}_0$ corresponds to no clustering, c.f.\ (ref), and that, in addition, there are $p$ potential levels of clustering of the data. All these levels of clustering are assumed to be nested from fine to increasingly more coarse clustering.
In this situation, following the IU statistical principle mentioned above, we reject clustering at level $m$ if and only if levels $0,\ldots,m$ are all rejected. That is, the natural testing strategy here is to test clustering at level $m$ against the coarser level $m+1$ sequentially, for $m=0,1,\dots,p-1$, and choose the level of clustering in the first non-rejected test. Algorithmically, we perform the following sequential testing procedure.
We can equivalently state the sequential testing problem in (ref) as a type of estimation problem. Specifically,
Of course, the $\hat{m}$ resulting from (ref) and from \hyperref[{mhat}]{\tagform@{\ref*{mhat}}} will be identical.
Because each individual test will reject a false null hypothesis with probability converging to one, this procedure will (at least asymptotically) never choose a level of clustering that is too fine. In other words, $\hat{m}$ defined in either (ref) or \hyperref[{mhat}]{\tagform@{\ref*{mhat}}} is (nearly) consistent. Precise asymptotic properties of the proposed sequential procedure are established in (ref), and finite\kern 0.04167em-sample performance is investigated by Monte Carlo simulation methods in (ref).
To our knowledge, only two other tests for the appropriate level of clustering have been proposed. Like our $\tau_\sigma$ test, these concern the standard error for a single coefficient. The best-known of them is due to \citet*{Ibragimov_2016}, and we refer to it as the IM test. It is a one\kern 0.04167em-sided test that is derived from the procedure of \citet*{Ibragimov_2010}. The IM test is based on the assumption that $G$ is fixed while the number of observations tends to infinity. In this respect, it differs from our score\kern 0.04167em-variance tests, for which the asymptotic theory in (ref) requires that $G\to\infty$. However, see the discussion in (ref).
There are two versions of the IM test. The first involves estimating the model separately for every coarse cluster. Unfortunately, this is impossible to do for models that involve treatment effects whenever the treatment is invariant within clusters. Even with treatment at the fine\kern 0.04167em-cluster level, it may not be possible to estimate the model for every coarse cluster. This is the case, for example, in the empirical example of (ref). Ibragimov_2016 therefore also propose a two\kern 0.04167em-sample version of their test statistic that can be used for testing the level of clustering for treatment models when entire clusters are treated or not treated. Differences between estimates for treatment and control clusters can be used to estimate the treatment effects and perform a test of fine clustering.
Another test for the appropriate level of clustering for a single coefficient has very recently been proposed by Cai_2022. This test is based on randomization inference. Like the IM test, and unlike our $\tau_\sigma$ test, it is necessarily one\kern 0.04167em-sided and treats $G$ as fixed. It also requires that the fine clusters be large, so that, unlike both the IM test and our tests, it cannot be used to test the null hypothesis of independence at the observation level.
Cai_2022 presents results from a number of simulation experiments for his test, the IM test, and the bootstrap version of our $\tau_\sigma$ test. They suggest that the IM test and our test are much more similar to each other than they are to Cai's test. The IM test always rejects more often than the bootstrap version of our $\tau_\sigma$ test, both when the null hypothesis is true and when it is false. It can over-reject quite severely in some cases, especially when the ratio of fine to coarse clusters is small. In the first version of this paper, we presented some figures comparing rejection frequencies for the IM test and the $\tau_\sigma$ test. In the interests of space, however, we have omitted these results, because they are broadly similar to those from the experiments in Cai_2022.
At this point, what is known about the properties of our tests, the IM test, and Cai's test suggests that none of them is to be preferred in every case. They can all provide useful information about the appropriate level at which to cluster. Two attractive features of our tests, which are not shared by the other two, are that the $\tau_\sigma$ test can be either one\kern 0.04167em-sided or two\kern 0.04167em-sided and that the $\tau_\Sigma$ test is based on more than one coefficient of interest.
The ultimate purpose of using any test for the appropriate level of clustering is to make more reliable inferences about the coefficient(s) of interest, that is, some element(s) of ${\bm\beta}$ in \hyperref[{model}]{\tagform@{\ref*{model}}}. This may or may not involve some sort of formal pre\kern 0.04167em-testing or model averaging procedure.
For simplicity, suppose we are attempting to construct a confidence interval for $\beta_1$, the (scalar) coefficient of interest, when there are just two levels of clustering, fine and coarse. Without a testing procedure, an investigator must choose between fine and coarse clustering on the basis of prior beliefs about which level is appropriate. With a testing procedure like the ones proposed in this paper, an investigator can instead choose the level of clustering based on the outcome of a test. This involves choosing a level $\alpha$ for the test and deciding whether to use a one\kern 0.04167em-sided or a two\kern 0.04167em-sided test. We then form the interval based on coarse clustering when the test rejects, and we form the interval based on fine clustering when it does not reject.
Of course, this procedure can never work as well as the infeasible procedure of simply choosing the correct level of clustering. It inevitably suffers from some of the classic problems associated with pre\kern 0.04167em-testing LP_2005. When there is actually fine clustering, the pre\kern 0.04167em-test will sometimes make a Type I error and reject, leading to an interval that is usually too long. When there is actually coarse clustering, the pre\kern 0.04167em-test will sometimes make a Type II error and fail to reject, leading to an interval that is usually too short. In (ref), we report the results of some simulation experiments that compare confidence intervals based on several alternative procedures.
As Ibragimov_2016 point out, it probably makes sense to report confidence intervals for $\beta_1$ based on all clustering levels that appear plausible. The resulting inferences are then explicitly conditional on the level of clustering. Because our tests, like the other ones discussed in (ref), provide evidence on the plausibility of each level of clustering, they can reduce the number of intervals that need to be reported. For example, if the hypothesis of independence is strongly rejected against one or more clustering structures, then it would not be necessary to report a confidence interval based on a heteroskedasticity-robust standard error. But if the $P$ value for the hypothesis of fine clustering against coarse clustering is neither extremely small nor very large, then it might well seem reasonable to report confidence intervals based on both levels. Whatever intervals an investigator chooses to report, tests for the appropriate clustering level can provide valuable information about which ones are empirically more plausible. These tests may thus be thought of as robustness checks Cai_2022.
In (ref), we derive the asymptotic distributions of the two score\kern 0.04167em-variance test statistics under the null hypothesis and show that they are divergent under the alternative. Then we prove the validity of the bootstrap implementation ((ref)) and prove asymptotic results for the sequential testing procedure ((ref)). We first state and discuss the assumptions needed for our proofs, which may be found in (ref).
(ref) is the assumption of (at most) “fine” clustering, which implies that the null hypothesis in \hyperref[{hypotheses}]{\tagform@{\ref*{hypotheses}}} is satisfied, even without taking the limit. In fact, it is slightly weaker than that, because we do not make the stronger assumption that all observations in any fine cluster are independent of those in a different fine cluster; we only assume that the cluster sums are independent across fine clusters. The moment conditions in (ref) and the multicollinearity condition in (ref) are standard in linear regression models.
Next, the conditions in (ref) rule out degenerate cases. The minimum eigenvalue condition rules out perfect negative correlation between scores within fine clusters. The maximum eigenvalue condition ensures that the variance of a single fine cluster cannot dominate the sum of the variances within a coarse cluster. It is basically satisfied if $M_g > 1$ for all $g$. The latter holds by construction of the test statistics, because any coarse cluster with $M_g=1$ will not contribute to $\hat{\bm{\theta}}$, and hence not to the test statistic.
The conditions in (ref) restrict the amount of heterogeneity of cluster sizes that is allowed under both the null and the alternative. Neither the fine cluster sizes nor the coarse cluster sizes are required to be bounded under these conditions, which allow the cluster sizes to diverge with the sample size. The first condition is used in the proofs to replace residuals with disturbances and for convergence of the variance. The second condition trades off moments and cluster size heterogeneity to rule out the possibility that one cluster dominates the test statistic in the limit in such a way that the central limit theorem does not apply; technically, it is used to verify Lyapunov's condition for the central limit theorem. When $\lambda \to \infty$, the second condition is implied by the first.
The denominators of both terms in (ref) show that these conditions trade off intra-cluster dependence and cluster-size heterogeneity. That is, the greater the amount of intra-cluster correlation, the larger are the denominators in (ref), which allows larger clusters without dominating the limit; a similar tradeoff was found in DMN_2019. Furthermore, more homogeneity in cluster sizes allows for fewer and larger clusters. We illustrate these tradeoffs in the following remarks.
The asymptotic distributions of the test statistics follow immediately from (ref).
We next consider the asymptotic behavior of the test statistics under the alternative. Because (ref) implies that \H{0} is true, we do not make that assumption. Instead, we impose the following conditions:
\goodbreak
\goodbreak
(ref) is the assumption of (at most) coarse clustering. This assumption is very general, and departures from the null could be very small and inconsequential. In order for our tests to be able to detect departures from the null hypothesis, with probability converging to one in the limit, we need to impose sufficient correlation within the coarse clusters. That is, we need ${\bm{\Sigma}}_g = \sum_{h_1=1}^{M_g}\sum_{h_2=1}^{M_g}{\rm E} ({\bm{s}}_{gh_1}{\bm{s}}_{gh_2}^\top)$ to be sufficiently large, in aggregate. This condition is embodied in (ref).
It follows immediately from (ref) that tests based on either of our statistics reject with probability converging to one under the alternative. That is, they are consistent tests.
Of course, power will depend in a complicated way on many aspects of the model and DGP, including the number of large clusters and their sizes, because these will affect the number of correlations that need to be estimated between fine clusters within coarse clusters; see \hyperref[{thetaknown}]{\tagform@{\ref*{thetaknown}}}. Power will also depend on the true values of these correlations. If they are mostly non-zero and non-trivial, then power will be higher with larger coarse clusters.
We now demonstrate the asymptotic validity of the bootstrap implementation of our tests. To this end, let $\tau$ denote either of our statistics, and let the cumulative distribution function of $\tau$ under \H{0} be denoted $P_0 (\tau \leq x)$. The corresponding bootstrap statistic is denoted $\tau^\ast$\kern -.08333em. As usual, let $P^\ast$ denote the bootstrap probability measure, conditional on a given sample, and let ${\rm E}^\ast$ denote the corresponding expectation conditional on a given sample.
First, note that the bootstrap theory requires a slight strengthening of the moment condition since at least four moments are now required. Second, (ref) shows that the bootstrap $P$ values in (ref) are asymptotically valid under (ref) and \H{0}. Third, note that neither the null hypothesis nor (ref) is imposed in (ref). Thus (ref) together show immediately that the bootstrap tests are consistent. We summarize these results in the following corollary.
The next theorem provides theoretical justification for the sequential testing procedure given in (ref).
The results in (ref) show that $\hat m$ defined in (ref) or \hyperref[{mhat}]{\tagform@{\ref*{mhat}}} is asymptotically correct with probability converging to $1-\alpha$ when $m_0 \leq p-1$ and with probability converging to 1 when $m_0 =p$. It is worth emphasizing that the sequential procedure will never “under-estimate” the clustering level, at least asymptotically, because $\hat m < m_0$ with probability converging to 0.
In the literature on cluster-robust inference, a few authors have considered an alternative asymptotic framework, referred to as fixed-$G$ asymptotics, in which the number of clusters is fixed as $N\to\infty$ while cluster sizes diverge; key early papers are Ibragimov_2010 and \citet*{BCH_2011}. However, fixed-$G$ asymptotics are proven under the very strong assumption that a central limit theorem applies to the normalized scores for each cluster. This assumption seriously limits the amount of intra-cluster dependence. For example, it rules out common models such as many types of random-effects and factor models. See MNW-guide for a detailed discussion.
Nonetheless, we now briefly consider an asymptotic framework in which the number of coarse clusters, $G$, is fixed, but there are many fine clusters within each coarse cluster, i.e.\ $M_g \to\infty$ for all $g$. For simplicity, we consider the scalar case with $k=1$. Let $\sigma_g^2 = \operatorname{Var} (\sum_{h=1}^{M_g}s_{gh})=\sum_{h=1}^{M_g}\sigma_{gh}^2$ (under the null) and define the weights $w_g^2 = \lim_{M_g\to\infty} \sigma_g^2 / \operatorname{Var} (\theta)^{1/2}$, where $\operatorname{Var} (\theta )$ is given in \hyperref[{tauvar}]{\tagform@{\ref*{tauvar}}}. Then suppose, for all $g$, that (i) $\sigma_g^{-1}\sum_{h=1}^{M_g}s_{gh}\overset{d} \longrightarrow{\rm N} (0,1)$, (ii) $\sigma_g^{-1}\sum_{h=1}^{M_g}s_{gh}^2\overset{P} \longrightarrow 1$, and (iii) $w_g^2 \in [0,\infty )$. The high-level condition (i) is typical of the fixed-$G$ literature and imposes very strong limitations on the amount of intra-cluster dependence that is allowed. Condition (ii) is a homogeneity assumption, and condition (iii) ensures that one cluster does not dominate the sum in the limit. Under the null hypothesis and these conditions, it can be proven that
where $\chi_{1,g}^2$ for $g=1,\ldots ,G$ denote independent $\chi_1^2$ random variables. Under suitable additional regularity conditions, we conjecture that $\tau_\sigma = \hat\theta / \sqrt{\widehat\operatorname{Var} (\hat\theta )}$ has the same asymptotic distribution as in \hyperref[{chi square}]{\tagform@{\ref*{chi square}}}.
The limiting distribution in \hyperref[{chi square}]{\tagform@{\ref*{chi square}}} is a weighted sum of independent $\chi_1^2$ random variables. Because the weights $w_g^2$ depend on unknown parameters, the distribution is non-pivotal and hence cannot be used for inference. Under the extreme homogeneity condition that the $w_g^2$ are the same for all $g$, the distribution simplifies to $(\chi_G^2-G)/\sqrt{2G}$, which is a centered and normalized $\chi_G^2$. This distribution is pivotal and could be used for inference, although the conditions under which it is derived are extraordinarily strong.
Continuing with this type of fixed-$G$ asymptotic argument, we could instead assume that $N_{gh}\to\infty$ for all $g,h$, while the $M_g$ and $G$ are fixed. That is, the number of observations within each fine cluster diverges, but there are only a fixed number of fine and coarse clusters. This setup is quite similar to the previous one. We conjecture that the asymptotic distribution would again be a weighted sum of $\chi^2_1$ random variables similar to the one in \hyperref[{chi square}]{\tagform@{\ref*{chi square}}}, but the summation would extend over $\sum_{g=1}^G M_g = G_{\rm f}$ elements.
In either case, if the weights $w_g^2$ are not too heterogeneous, the fixed-$G$ limiting distributions could be well approximated by a standard normal distribution, at least when the number of clusters is not very small. In the setup with $M_g \to\infty$, this would be the number of coarse clusters, $G$. In the setup with $N_{gh}\to\infty$, it would be the number of fine clusters, $G_{\rm f}$. Thus, in the end, the normal limit theory obtained under large\kern 0.04167em-$G$ asymptotics in (ref) and (ref) may also provide a good approximation under fixed-$G$ asymptotics. A full analysis of fixed-$G$ asymptotic theory for our model and test statistics would be interesting, but it is beyond the scope of this paper and is consequently left for future work.
As discussed in (ref), it is commonly the case in empirical work that the number of regressors is very large and that most of the regressors are not of primary interest. Comparing large-dimensional CRVE matrices by the methods in (ref) is impractical. Fortunately, it is easy to solve this problem by partialing out the regressors that are not of primary interest prior to performing our tests.
Suppose the full set of regressors is partitioned as ${\bm{X}} =[{\bm{X}}_1,\; {\bm{X}}_2]$, where ${\bm{X}}_1$ denotes the $N\times k_1$ matrix of the regressors of interest and ${\bm{X}}_2$ denotes the $N \times k_2$ matrix of other regressors, with $k = k_1 + k_2$. Similarly, partition ${\bm\beta}^\top =[{\bm\beta}_1^\top\kern -.08333em, \; {\bm\beta}_2^\top ]$, where the coefficients corresponding to the regressors of interest are in the $k_1 \times 1$ parameter vector ${\bm\beta}_1$ and the rest are collected in ${\bm\beta}_2$. If the coefficient vector of interest is actually a linear combination of the elements of ${\bm\beta}_1$ and ${\bm\beta}_2$, we can redefine ${\bm{X}}$ as a nonsingular affine transformation of the original ${\bm{X}}$ matrix, so that ${\bm\beta}_1$ has the desired interpretation.
We regress each column of ${\bm{X}}_1$ on ${\bm{X}}_2$ and define ${\bm{Z}}$ as the matrix of residuals from those $k_1$ regressions. The model \hyperref[{model}]{\tagform@{\ref*{model}}} can then be rewritten as
where ${\bm{M}}_{{\bm{X}}_2}={\bf I}_N -{\bm{X}}_2 ({\bm{X}}_2^\top{\bm{X}}_2)^{-1}{\bm{X}}_2^\top$ is the orthogonal projection matrix that projects off (or partials out) ${\bm{X}}_2$. The regressor matrices ${\bm{Z}}$ and ${\bm{X}}_2$ are orthogonal, and the models \hyperref[{model}]{\tagform@{\ref*{model}}} and \hyperref[{newmodel}]{\tagform@{\ref*{newmodel}}} have exactly the same explanatory power and the same disturbances, ${\bm{u}}$. The coefficient ${\bm\beta}_1$ in \hyperref[{newmodel}]{\tagform@{\ref*{newmodel}}} is identical to the one defined in the previous paragraph, but the coefficient ${\bm{\delta}}$ is different from ${\bm\beta}_2$.
Using the orthogonality between ${\bm{Z}}$ and ${\bm{X}}_2$, the OLS estimate of ${\bm\beta}_1$ is, c.f.\ \hyperref[{betahat}]{\tagform@{\ref*{betahat}}},
where ${\bm\beta}_{1,0}$ is the true value of ${\bm\beta}_1$. The relationship between ${\bm{Z}}$ and ${\bm{X}}$ can be written as
Therefore, the score for ${\bm\beta}_1$ is ${\bm{Z}}_g^\top {\bm{u}}_g = {\bm{Q}}^\top {\bm{X}}_g^\top {\bm{u}}_g = {\bm{Q}}^\top {\bm{s}}_g$. Thus, from \hyperref[{beta1hat}]{\tagform@{\ref*{beta1hat}}} and \hyperref[{Z and Q}]{\tagform@{\ref*{Z and Q}}}, we obtain the following sandwich formula, c.f.\ \hyperref[{covbeta}]{\tagform@{\ref*{covbeta}}},
Under (ref), ${\bm{Q}} \overset{P} \longrightarrow [ {\bf I}_{k_1} , \; -{\bm{\Xi}}_{12}{\bm{\Xi}}_{22}^{-1} ]^\top = {\bm{A}}$, say, so that the middle matrix in \hyperref[{covbeta1}]{\tagform@{\ref*{covbeta1}}} is clearly an estimator of ${\bm{A}}^\top {\bm{\Sigma}} {\bm{A}}$.
The matrix ${\bm{Q}}$ in \hyperref[{Z and Q}]{\tagform@{\ref*{Z and Q}}} and its limit ${\bm{A}}$ can be viewed as mechanisms for dimension reduction. They transform the problem from one involving the $k\times k$ matrix ${\bm{\Sigma}}$ to one involving the $k_1\times k_1$ matrix ${\bm{A}}^\top {\bm{\Sigma}} {\bm{A}}$. The latter is the variance of ${\bm{A}}^\top {\bm{X}}^\top {\bm{u}}$, and it depends on the clustering structure in the same way as ${\bm{\Sigma}}$. For the model \hyperref[{newmodel}]{\tagform@{\ref*{newmodel}}}, we consequently replace the hypotheses in \hyperref[{hypotheses}]{\tagform@{\ref*{hypotheses}}} with
Furthermore, from \hyperref[{Z and Q}]{\tagform@{\ref*{Z and Q}}}, we see that we can use the same algebra for the model in \hyperref[{newmodel}]{\tagform@{\ref*{newmodel}}} as for the model in \hyperref[{model}]{\tagform@{\ref*{model}}} to define the test statistics, i.e.\ \hyperref[{Sigmac}]{\tagform@{\ref*{Sigmac}}}, \hyperref[{Sigmaf}]{\tagform@{\ref*{Sigmaf}}}, and so on, but now with empirical scores ${\bm{Q}}^\top \hat{\bm{s}}_g$ and ${\bm{Q}}^\top \hat{\bm{s}}_{gh}$ instead of $\hat{\bm{s}}_g$ and $\hat{\bm{s}}_{gh}$, respectively. This also applies to the bootstrap implementation in (ref). Of course, degrees-of-freedom corrections like the factor $m_c$ in \hyperref[{Sighat}]{\tagform@{\ref*{Sighat}}} need to reflect the total number of estimated coefficients.
Since very few regression models in economics contain just one regressor, the $\tau_\sigma$ test will almost always involve partialing out. It seems likely that the $\tau_\Sigma$ test will also involve partialing out in the vast majority of cases, so that the dimension of the vector $\hat{\bm{\theta}}$ upon which the $\tau_\Sigma$ test is based will be $k_1$ rather than $k$.
Most of the papers cited in the second paragraph of (ref) employ simulation experiments to study the finite\kern 0.04167em-sample properties of methods for cluster-robust inference. To our know\-ledge, all of these papers use some sort of random-effects, or single\kern 0.04167em-factor, model to generate the data. The key feature of these models is that all of the intra-cluster correlation for every cluster $g$ arises from a single random variable, say $\xi_g$, which affects every observation within that cluster equally. This yields disturbances that are equi-correlated within each cluster.
Although this type of DGP is convenient to work with and can readily generate any desired level of intra-cluster correlation, it cannot be used when a regression model has cluster fixed effects. Because the fixed effects completely explain the $\xi_g$, the residuals are always uncorrelated. Thus, for models with cluster fixed effects, it is always valid to use heteroskedasticity-robust (HR) standard errors whenever the intra-cluster correlation of the disturbances arises solely from a random-effects model. In such cases, the null hypothesis of our tests is satisfied, and they will have no (asymptotic) power. Of course, this is the desired outcome both in the statistical sense, because the null is satisfied, and in the practical sense, because cluster-robust (CR) standard errors are not needed.
In practice, HR and CR standard errors often differ greatly in models with cluster fixed effects; see, for example, BDM_2004, JGM-CJE, and (ref). Therefore, whatever processes are generating intra-cluster correlation in real-world data must be more complicated than simple random-effects models. Since we wish to investigate models with cluster fixed effects, we need to employ a DGP for which cluster fixed effects do not remove all of the intra-cluster correlation. To this end, we generate both the regressors and the disturbances in our experiments using factor models of the form
Here $\xi^1_g$ and $\xi^2_g$ are random effects, distributed as standard normal, which apply respectively to the odd-numbered and even-numbered observations within the \th{g} cluster. The $\zeta_{gi}$ are also distributed as standard normal. Under the DGP \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}}, the $z_{gi}$ have variance one, and the intra-cluster correlation of the odd (or even) observations is $\rho \ge 0$.
The DGP \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}} can be interpreted in a variety of ways, depending on the nature of the data. The idea is that there are two types of observations within each cluster, and all the intra-cluster correlation is within each type. For example, with clustering at the geographical level, there might be two sub\kern 0.04167em-regions. With clustering at the industry level, there might be two types of firm. The key assumption is that the researcher knows which cluster an observation belongs to, but not which type. Including cluster fixed effects explains some of the intra-cluster correlation by estimating an average of $\xi^1_g$ and $\xi^2_g$ for each cluster, but it does not explain all of it. Thus cluster-robust inference is still needed, and our tests should still have power.
In practice, of course, there might be more than than two types within each cluster, and the numbers of observations in each would almost certainly not be the same. It would be easy to make the DGP \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}} more complicated. However, our objective is not to mimic any actual dataset, but simply to generate data in a way that allows cluster fixed effects to be combined with cluster-robust standard errors.
The DGP \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}} makes no reference to fine and coarse clusters. It could be used to generate either finely or coarsely clustered data. The regressors ${\bm{X}}_1$ (that is, the ones whose coefficients are of interest; see (ref)) are generated using \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}}, and they are always coarsely clustered. This ensures that, if the disturbances are either independent ($\rho=0$), finely clustered, or coarsely clustered, the scores are also independent, finely clustered, or coarsely clustered, respectively.
In all experiments, each of the regressors in ${\bm{X}}_1$ is generated independently. This implies that there is no correlation among the coefficient estimates. It might seem that the extent of any such correlation would be important for the properties of the $\tau_\Sigma$ tests. However, that is not the case. We find numerically that the $\tau_\Sigma$ statistic is invariant to any transformation of ${\bm{X}}_1$ that does not change the subspace spanned by its columns. Thus there is no loss of generality in generating the columns of ${\bm{X}}_1$ independently.
Our first set of experiments is designed to investigate the rejection frequencies of asymptotic and bootstrap score\kern 0.04167em-variance tests under the null hypothesis. The model is
where the regressors $X^\ell_{ghi}$ are generated independently across $\ell$ by \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}} at the coarse level with $\rho=0.5$. The additional regressors in ${\bm{X}}^2_{gh}$ are either a constant term or a set of cluster fixed effects. When testing fine against coarse clustering, the fixed effects are at the fine level, and the disturbances are finely clustered with $\rho=0.1$. When testing independence against (coarse) clustering, the fixed effects are at the coarse level, and the disturbances are independent. The number of coarse clusters, which in this section we denote by $G_{\rm c}$, is allowed to vary. In the first set of experiments, there are always four fine clusters in each coarse cluster, so that $G_{\rm f} = 4G_{\rm c}$.
(ref) plots rejection frequencies at the 0.05 level for $\tau_\Sigma$ tests against $G_{\rm c}$, which varies from 6 to 36. We started at $G_{\rm c}=6$ to avoid singularities when $k_1=5$ and stopped at $G_{\rm c}=36$ because the results were hardly changing at that point. The values $k_1=1,\ldots,5$ imply that the number of degrees of freedom for the tests is 1, 3, 6, 10, or 15. Panels (a) and (c) concern tests of fine clustering against coarse clustering, and panels (b) and (d) concern tests of independence against clustering. The top two panels report rejection frequencies for asymptotic tests at the 0.05 level, and the bottom two report comparable ones for bootstrap tests. Notice that the vertical axes for the asymptotic tests are much longer than the ones for the bootstrap tests, because the latter work very much better.
One striking feature of (ref) is that, for the asymptotic tests, over-rejection increases sharply with $k_1$. This should not have been a surprise in view of the fact that, like the information matrix test White_1982, the $\tau_\Sigma$ test has degrees of freedom that are $O(k_1^2)$. DM_1992 found a similar tendency for the rejection rate of the information matrix test (in particular, the popular $N\kern -.08333em R^2$ form of it) to increase rapidly with the number of coefficients being tested.
When $G_{\rm c}$ is small, asymptotic tests of fine against coarse clustering, in Panel (a), over-reject more severely than tests of independence, in Panel (b). When $k_1=1$, there is almost no over-rejection for the tests of independence in Panel (b). For $k_1\ge2$, there is also more over-rejection in Panel (a) than in Panel (b) when $G_{\rm c}=6$, but the over-rejection diminishes much more rapidly as $G_{\rm c}$ increases in Panel (a) than in Panel (b).
The bootstrap versions of the tests perform very much better than the asymptotic ones. There is slight over-rejection in Panel (c) for smaller values of $G_{\rm c}$, which is really only noticeable for $k_1=1$ and $k_1=2$. In Panel (d), the bootstrap tests of independence work perfectly, except for experimental errors.
The bootstrap tests can be computationally demanding when the sample size is large, particularly for larger values of $k_1$. This is especially true for tests where the null hypothesis is no clustering, because the calculations in \hyperref[{Sigmac}]{\tagform@{\ref*{Sigmac}}}, \hyperref[{Sigmaf}]{\tagform@{\ref*{Sigmaf}}}, and \hyperref[{var2sided}]{\tagform@{\ref*{var2sided}}} involve score vectors of which the size is the number of clusters under the null hypothesis. This number is $N$ for tests of no clustering but only $G_{\rm f}$ for tests of fine clustering.
In (ref), we hold the number of coarse clusters constant at $G_{\rm c}=8$ and allow either the number of fine clusters per coarse cluster or the $N_g$ to vary. Results are shown for two specifications of \hyperref[{simmod}]{\tagform@{\ref*{simmod}}}. For the first of these, there are fine fixed effects when the null is fine clustering and cluster fixed effects when the null is independence, as in (ref). For the second, there is just a constant term. To make the figure readable, results are shown only for $k_1=1$, 3, and 5.
In Panels (a) and (c), the horizontal axis shows the number of fine clusters per coarse cluster, which varies between 3 and 12, so that the total number of fine clusters varies between 24 and 96. The rejection frequencies for asymptotic tests of fine against coarse clustering drop somewhat as $G_{\rm f}/G_{\rm c}$ increases. For the model with fixed effects, the asymptotic tests for $k_1=1$ work almost perfectly for $G_{\rm f}/G_{\rm c} \ge 8$, and all the bootstrap tests work almost perfectly for $G_{\rm f}/G_{\rm c} \ge 6$. The asymptotic tests always reject less often for the model with a constant term than for the model with fixed effects. For $k_1=1$, the former actually under-reject modestly for larger values of $G_{\rm f}/G_{\rm c}$.
In Panels (b) and (d), the horizontal axis shows the number of observations per coarse cluster, which varies between 25 and 400, on a log scale. It is evident that the asymptotic tests of independence perform better as the clusters become larger, although the curves are pretty flat at $N_g=400$. The asymptotic tests with just a constant over-reject much less than the tests with fixed effects. When $k_1=1$, these tests under-reject for all values of $N_g$. All the bootstrap tests work essentially perfectly.
(ref) shows rejection frequencies for both upper-tail and two\kern 0.04167em-sided $\tau_\sigma$ tests. The experimental design is essentially the same as for (ref), except that, since $k_1=1$, results for $G_{\rm c}=3$ and $G_{\rm c}=4$ are included. The asymptotic upper-tail tests over-reject noticeably more often than the asymptotic two\kern 0.04167em-sided tests. In contrast, the bootstrap upper-tail and two\kern 0.04167em-sided tests perform identically (and extremely well). Thus it seems to be valuable to bootstrap both types of $\tau_\sigma$ test, but particularly important to bootstrap upper-tail tests.
In the next set of experiments, we turn our attention to power, focusing on the special case of the $\tau_\sigma$ test for a single coefficient. The data are generated by \hyperref[{simmod}]{\tagform@{\ref*{simmod}}}, with one regressor and coarse fixed effects. As usual, the regressor is generated by \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}} with coarse clustering and $\rho=0.5$. The disturbances are generated by the same model, with $\rho$ varying between 0.00 and 0.10. We report results only for bootstrap tests with $B=999$. Using 999 instead of 399 reduces the, already quite small, power loss caused by using a finite number of bootstrap samples DM_2000.
(ref) shows the power of either two or three types of bootstrap $\tau_\sigma$ tests against coarse clustering as a function of the value of $\rho$ for the disturbances. The three types are upper-tail, symmetric, and equal-tail. As can be seen in both panels, all tests reject extremely close to 5% of the time when the null hypothesis is true. In Panel (a), the null hypothesis is fine clustering for three different values of $G_{\rm f}$. Power increases greatly when the number of fine clusters goes from 20 to 40. It increases further, but much more modestly, when $G_{\rm f}$ goes from 40 to 100. The upper-tail tests are more powerful than the symmetric ones, but only slightly more when $G_{\rm f}=100$. To avoid making the figure unreadable, Panel (a) omits the equal-tail tests, which have much less power than the other two types of tests.
In Panel (b), the null hypothesis is independence, and the alternative is clustering with 10 (coarse) clusters. There are three bootstrap tests for a model with just a constant term and three tests for a model with cluster fixed effects. As expected, the tests are more powerful when there is just a constant term, since the fixed effects explain some of the intra-cluster correlation. For each set of tests, the upper-tail test is slightly more powerful than the symmetric test, which in turn is substantially more powerful than the equal-tail test.
The results in Panel (b) illustrate the fact that it generally makes no sense to use equal-tail bootstrap SV tests. These tests are designed to reject equally often in each tail under the null hypothesis. Since the mean of the test statistics under the null is positive in our experiments, the equal-tail test implicitly uses asymmetric critical values, with the positive one being larger in absolute value than the negative one. This reduces its power against $\sigma^2_{\rm c} > \sigma^2_{\rm f}$, which is precisely the alternative we want SV tests to have power against.
Up to this point, all the simulations have involved equal-sized clusters. There are many ways in which coarse-cluster sizes, fine-cluster sizes, and the numbers of fine clusters per coarse cluster could vary. We next allow cluster sizes to vary for one level of clustering. (ref) considers $\tau_\sigma$ tests of no clustering and plots rejection frequencies against a measure of cluster size variation. The $N$ observations are allocated among $G$ clusters using the equation
where $\delta\ge0$, $[\cdot]$ denotes the integer part of its argument, and $N_G = N - \sum_{g=1}^{G-1} N_g$. This scheme has been used in MW-JAE, DMN_2019, and several other papers. In the experiments of (ref), $G=10$ and $N=1000$. When $\delta=0$, $N_g=100$ for all $g$. For $\delta=1$, the $N_g$ range from 61 to 155; for $\delta=2$, from 34 to 213; and for $\delta=4$, from 9 to 340. There is one regressor and 10 cluster fixed effects.
In Panel (a) of (ref), the null hypothesis of no clustering is true. The upper-tail asymptotic test over-rejects noticeably for small values of $\delta$, but rejection frequencies decline as $\delta$ increases, and they are less than 0.05 for $\delta=4$. In contrast, the upper-tail bootstrap test rejects almost exactly 5% of the time for all values of $\delta$. In Panel (b), the null hypothesis is false. Both tests have substantial power when $\delta$ is small, but it falls as $\delta$ increases. This makes sense, because the total number of off-diagonal elements in all the clusters increases with $\delta$, causing the number of terms in the variance \hyperref[{tauvar}]{\tagform@{\ref*{tauvar}}} to increase. The asymptotic test has noticeably more power than the bootstrap test for small values of $\delta$, but it has less power for the largest values, where it under-rejects under the null. The power differences almost certainly just reflect the size distortions of the asymptotic tests.
The results in (ref) suggest that the finite\kern 0.04167em-sample performance of SV tests inevitably depends on the pattern of cluster sizes, although probably much less for bootstrap tests than for asymptotic ones.
Our next set of experiments concerns the sequential testing procedure of (ref), using bootstrap tests. These experiments are quite similar to the ones in (ref), except that there are 8 coarse clusters, 48 fine clusters, and 2400 observations. As in (ref), there are 999 bootstrap samples. The model always contains coarse\kern 0.04167em-level fixed effects, and all of the tests are at the 0.05 level. The figure shows the outcomes of sequential, upper-tail $\tau_\sigma$ tests as $\rho$, the intra-cluster correlation for each set of disturbances generated by \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}}, varies within either coarse or fine clusters.
In Panel (a) of (ref), there is coarse clustering in the DGP, except when $\rho=0$. In that case, as expected, the procedure chooses no clustering (N) almost exactly 95% of the time, fine clustering (F) almost exactly 4.75% of the time, and coarse clustering (C) almost exactly 0.25% of the time. These results illustrate why the sequential testing algorithm does not inflate the Type I error. In this case, the true null is rejected almost exactly $\alpha$% of the time. Amongst the replications with false positives, the test concludes that fine clustering is appropriate about $(1-\alpha)$% of the time and that coarse clustering is appropriate the remaining $\alpha$% of the time.
As $\rho$ increases, the procedure chooses N or F less and less often. For very small values of $\rho$, it chooses N or F more often than C, but that changes quickly as $\rho$ increases. The gap between the solid red and dashed purple curves shows the fraction of the time that F is (incorrectly) chosen. This gap is always small, and it vanishes as $\rho$ becomes large.
The sequential procedure inevitably has less power than testing no clustering directly against coarse clustering. The outcome of testing N directly against C at the 0.05 level is shown by the blue dashed curve in Panel (a). The gap between this curve and the purple dashed curve that separates the F and C regions shows the power loss from using the sequential procedure. This power loss arises for two reasons. First, the test of N against F has less power than the test of N against C; see (ref). Second, even when N is correctly rejected against F, the latter is sometimes not rejected against C. When the investigator finds coarse clustering more plausible than fine clustering, it may therefore make sense to test no clustering directly against the former rather than to employ the sequential procedure.
In Panel (b) of (ref), there is fine clustering in the DGP, except when $\rho=0$. The sequential procedure again works very well. As $\rho$ increases, it incorrectly chooses no clustering a rapidly diminishing fraction of the time. For larger values of $\rho$, it incorrectly chooses coarse clustering about 5.2% of the time, because the bootstrap SV tests over-reject slightly with only 8 coarse and 48 fine clusters. Once again, the outcome of testing N directly against C is shown by the blue dashed line. This test works much less well than the sequential procedure, often failing to reject the false null hypothesis that the disturbances are not clustered. This is not surprising, since the alternative involves clustering at a coarser level than the DGP.
In (ref), we discussed several procedures for making inferences about a single regression coefficient when clustering may be either fine or coarse. We now investigate some of these procedures, notably pre\kern 0.04167em-test ones based on SV tests. There are four simulation experiments, each involving 12 coarse clusters. In two of them, we pre\kern 0.04167em-test the null of no clustering, and in the other two we pre\kern 0.04167em-test the null of fine clustering with 96 fine clusters.
The model is a variant of \hyperref[{simmod}]{\tagform@{\ref*{simmod}}}, with eight regressors plus coarse-level fixed effects, so that $k=G_{\rm c}+8=20$. The regressors are generated by \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}} with $\rho=0.5$. The disturbances $u_{ghi}$ are generated as a convex combination of two disturbances, $\epsilon^{\rm c}_{gi}$ and $\epsilon^{\rm f}_{ghi}$, with weights $\eta$ and $1-\eta$ respectively, rescaled so that the $u_{ghi}$ have unit variance. The $\epsilon^{\rm c}_{gi}$ are generated by \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}} with $\rho=0.25$. When the pre\kern 0.04167em-test null hypothesis is fine clustering, the $\epsilon^{\rm f}_{ghi}$ are generated in the same way as the $\epsilon^{\rm c}_{gi}$, but for 96 fine clusters instead of 12 coarse ones. When the pre\kern 0.04167em-test null hypothesis is no clustering, the $\epsilon^{\rm f}_{ghi}$ are i.i.d.\ normal.
The parameter $\eta$ determines the amount of correlation within coarse clusters. The pre\kern 0.04167em-test null hypotheses are true when $\eta=0$, so that there is either no intra-cluster correlation or only correlation within the fine clusters. The pre\kern 0.04167em-test null hypotheses are false when $\eta>0$, and the DGP moves further away from the pre\kern 0.04167em-test null as $\eta$ increases. In the experiments, we vary $\eta$ from 0 to 1.
There are several asymptotically valid standard errors for coarse clustering, fine clustering, and no clustering. The best-known variance matrix estimator with clustering, often referred to as CV$_1$, is the usual sandwich estimator \hyperref[{covbeta}]{\tagform@{\ref*{covbeta}}} with $\hat{\bm{\Sigma}}_{\rm c}$ given by \hyperref[{Sighat}]{\tagform@{\ref*{Sighat}}} or \hyperref[{Sigmac}]{\tagform@{\ref*{Sigmac}}}. However, recent work Hansen-jack,MNW-bootknife,MNW-influence suggests that the cluster jackknife, or CV$_3$, estimator usually performs better than CV$_1$, so we use the former for inference about the regression coefficient. For the case of no clustering, we use the HC$_3$ standard error of MW_1985, which is a jackknife estimator analogous to CV$_3$.
We focus on inference about $\beta_1$, one of the $\beta_\ell$ in \hyperref[{simmod}]{\tagform@{\ref*{simmod}}}. The pre\kern 0.04167em-test estimators that we study are based on upper-tail $\tau_\sigma$ tests. Upper-tail tests are more powerful than two\kern 0.04167em-sided tests, so that the former make fewer Type II errors; see (ref). Moreover, even when the difference between $\operatorname{Var}_{\rm c}(\hat\beta_1)$ and $\operatorname{Var}_{\rm f}(\hat\beta_1)$ is positive, $\widehat\operatorname{Var}_{\rm c}(\hat\beta_1)$ can be smaller than $\widehat\operatorname{Var}_{\rm f}(\hat\beta_1)$. This happens frequently in our experiments when $\eta$ is greater than 0 but small. Thus, investigators who do not wish to reject fine clustering in favor of coarse clustering when the coarse standard error is smaller than the fine one will choose to employ upper-tail pre\kern 0.04167em-tests.
The choice among various standard errors is an estimation problem. Thus, it seems reasonable to compare them on the basis of root mean squared error (RMSE). When the pre\kern 0.04167em-test null hypothesis is no clustering, the standard error is based on HC$_3$, CV$_3$, or the one chosen by pre\kern 0.04167em-tests at either the 0.05 or 0.20 level. When the pre\kern 0.04167em-test null is fine clustering, the standard error is based on fine CV$_3$, coarse CV$_3$, or the one chosen by pre\kern 0.04167em-tests at the same two levels. (ref) shows the RMSEs associated with each of these standard errors. In Panels (a) and (b), the pre\kern 0.04167em-test null hypothesis is no clustering. In Panels (c) and (d), it is fine clustering, with 96 clusters. There are 4800 observations in Panels (a) and (c) and 24,000 in Panels (b) and (d).
The HC$_3$ or fine CV$_3$ standard errors are the most accurate when $\eta=0$, and they continue to be the most accurate for small values of $\eta$. However, for larger values of $\eta$, they are by far the least accurate, because they are severely biased. In contrast, the coarse CV$_3$ standard errors are the least accurate when $\eta$ is small, but for moderate and larger values of $\eta$ they are the most accurate. The two pre\kern 0.04167em-test standard errors are substantially more accurate than the coarse CV$_3$ ones for small values of $\eta$ and almost identical to the latter for large values of $\eta$. In between, there is always a region where the pre\kern 0.04167em-test standard errors are slightly less accurate than the coarse CV$_3$ ones. This is barely noticeable for pre\kern 0.04167em-tests at the 0.20 level, but it is quite noticeable for pre\kern 0.08333em-tests at the 0.05 level, especially in Panel (c), where the SV tests have the least power.
In our view, the 0.20 pre\kern 0.04167em-test standard errors in (ref) perform substantially better than any of the others. They are much more accurate than coarse CV$_3$ standard errors for small values of $\eta$, slightly less accurate for some intermediate values, and essentially identical for larger values. Since using a more accurate standard error yields a confidence interval that provides a better sense of how reliable a coefficient estimate is, it seems reasonable to base confidence intervals on 0.20 pre\kern 0.04167em-test standard errors.
Of course, using a more accurate standard error does not guarantee better coverage. (ref) shows the coverage of confidence intervals using the four standard errors in (ref). The coarsely-clustered intervals always under-cover to some extent. With only 12 clusters, that is not surprising. If we had used CV$_1$ instead of CV$_3$ to construct the intervals, they would have under-covered to a somewhat greater extent. On the other hand, coverage would almost certainly have been closer to 95% if we had used the wild cluster bootstrap MNW-bootknife, but that would have been computationally very demanding to simulate. The coverage using HC$_3$ and the finely-clustered CV$_3$ is almost exactly 95% when $\eta=0$, but they always under-cover for $\eta>0$, and the under-coverage is very severe for most values of $\eta$. Indeed, their coverage always rapidly drops below 0.90, the lower limit of the vertical axis.
The pre\kern 0.04167em-test intervals over-cover slightly when $\eta=0$, which is a consequence of Type I errors in the pre\kern 0.04167em-tests. However, they under-cover more than the coarsely-clustered CV$_3$ intervals for intermediate values of $\eta$ because of Type II errors. The under-coverage is much more pronounced for pre\kern 0.04167em-tests at the 0.05 level than for pre\kern 0.04167em-tests at the 0.20 level. Because the sample size is five times larger in Panels (b) and (d) than in Panels (a) and (c), the pre\kern 0.04167em-tests are more powerful, and the pre\kern 0.04167em-test intervals converge more rapidly to the coarsely clustered CV$_3$ interval as $\eta$ increases.
To save computer time and programming effort, we use asymptotic SV tests in these experiments. In consequence, the levels of the pre\kern 0.04167em-tests are not exactly 0.05 and 0.20. In particular, the actual levels of tests at the 0.20 level are noticeably lower than 0.20, and the ones for tests at the 0.05 level are somewhat higher than 0.05. If we had used bootstrap pre\kern 0.04167em-tests, the under-coverage for moderate values of $\eta$ would have been a bit smaller for tests at the 0.20 level and a bit larger for tests at the 0.05 level. But all the curves for pre\kern 0.04167em-test confidence intervals would have looked very similar. They would also have looked very similar if we had used CV$_1$ and HC$_1$ instead of CV$_3$ and HC$_3$.
We now illustrate the use of our score\kern 0.04167em-variance tests in a realistic empirical setting. We employ the widely-used data from the Tennessee Student Teacher Achievement Ratio (STAR) experiment Finn_1990, Mosteller_1995. We use these data to estimate a cross\kern 0.04167em-sectional model similar to one in Krueger_1999. The STAR experiment randomly assigned students either to small-sized classes, regular-sized classes without a teacher's aide, or regular-sized classes with a teacher's aide. We are interested in the effect of being in a small class, or being in a class with an aide, on standardized test scores in reading.
We estimate the following cross\kern 0.04167em-sectional regression model:
The outcome variable $\text{read-one}_{sri} $ is the reading score in grade one of student $i$ in classroom $r$ in school $s$. We are interested in $\beta_s$ and $\beta_a$, which are the coefficients for the small-class and aide\kern 0.04167em-class dummies. Small-class equals 1 if a student attended a small class in grade one and equals 0 otherwise; aide\kern 0.04167em-class is constructed in the same way for classes with or without a teacher's aide. Additional control variables are collected in the vector of regressors ${\bm{x}}_{sri}$. These include dummy variables for whether the student was male, non-white, or received free lunches, as well as a dummy variable for whether the student's teacher was non-white. They also include the teacher's years of experience and the student's reading score in kindergarten. Finally, there are dummy variables for the student's quarter of birth, the student's year of birth, and the teacher's highest degree. There are thus 17 coefficients in total, not counting the constant term or the school fixed effects, if any.
OLS estimates for the model \hyperref[{eq:starcs}]{\tagform@{\ref*{eq:starcs}}} are presented in the top half of (ref). Two variants of the model are estimated. In the left panel, there is just a constant term. In the right panel, there are school fixed effects. It is impossible to use classroom fixed effects, because treatment was assigned at the classroom level. Three sets of standard errors and $t$-statistics are reported for each variant of the model. For each set, the first column reports results that are heteroskedasticity-robust (HR), using HC$_3$ standard errors. The next two columns report results that are cluster-robust (CR) at either the classroom (R) level or the school (S) level, using CV$_3$ standard errors. As in (ref), we employ HC$_3$ and CV$_3$, instead of the more commonly-used HC$_1$ and CV$_1$ estimators, because the former tend to yield more reliable inferences. The HR results would have been very similar if we had used HC$_1$ instead of HC$_3$. However, some of the CR results would have been noticeably different if we had used CV$_1$ instead of CV$_3$. The reason for this is interesting, and we discuss it below.
Because treatment was assigned at the classroom level, it seems plausible that clustering at that level would be appropriate. However, since there are multiple classrooms per school, and students from the same school probably have many common characteristics and peer effects, it might also seem natural to cluster at the school level instead of the classroom level; even more so if assignment was not entirely random.
Unfortunately, the dataset does not contain a classroom indicator. One was created by using the information on the school ID, teacher's race, teacher's experience, teacher's highest degree, teacher's career ladder stage, and treatment status. It is possible that this procedure occasionally grouped two classes into one class, when two teachers in the same school had exactly the same observable characteristics. However, since the largest observed class had only 29 students, this seems unlikely to have happened often. Moreover, it would not be a problem, because the true classes would always be nested within the larger, assumed class. What would be a problem is if classes were incorrectly partitioned, but this cannot happen.
For the model without school fixed effects, the estimated impact on test scores of being in a small class is $\hat\beta_s = 9.211$. Based on an HR standard error of 1.63, the $t$-statistic for the null hypothesis that $\beta_s=0$ is 5.64. When we instead use CR standard errors clustered at the classroom level, the standard error for $\beta_s$ increases to 3.23, and the $t$-statistic decreases to 2.81. Using CR standard errors clustered at the school level yields almost identical results; the standard error is 3.25, and the $t$-statistic is 2.83. In this case, the level at which we cluster makes no qualitative difference. For the model with school fixed effects, the estimate of the small-class effect is somewhat lower at $\hat\beta_s=8.095$. The HR $t$-statistic is now 5.20, the classroom-level CR $t$-statistic is 2.67, and the school-level CR $t$-statistic is 2.57. Once again, the level at which we cluster does not change the conclusions.
The estimated effect on test scores of being in a class with an aide is $\hat \beta_a = 6.245$ without school fixed effects and $\hat \beta_a = 4.170$ with them. Based on the HR $t$-statistics, there seems to be fairly strong evidence that $\beta_a \neq 0$ for both models. However, when we cluster at the classroom level, we cannot reject this null hypothesis at the 0.05 level for either specification. When we cluster at the school level, we can do so for the model without fixed effects ($P=0.031$), but not for the model with fixed effects.
The lower panel of (ref) shows the values of our SV test statistics, and the associated upper-tail asymptotic and bootstrap $P$ values, for the two coefficients of interest, both individually and jointly. It also shows results for the IM test for the model with school fixed effects, when that test can be calculated. For each specification, we consider three hypotheses: H$_{\rm N}$ is no clustering with possible heteroskedasticity, H$_{\rm R}$ is classroom-level clustering, and H$_{\rm S}$ is school-level clustering. These are nested as $\text{H}_{\rm N} \subseteq \text{H}_{\rm R} \subseteq \text{H}_{\rm S}$.
For testing H$_{\rm N}$ against H$_{\rm R}$, the SV tests, both asymptotic and bootstrap, very strongly reject the null in all cases. IM tests cannot be computed for this hypothesis, because the procedure requires the model to be estimated classroom by classroom, and the two treatment variables are invariant at that level. For testing H$_{\rm N}$ against H$_{\rm S}$, the SV tests also very strongly reject the null in all cases. This is not surprising. Since there is overwhelming evidence against H$_{\rm N}$ when tested against H$_{\rm R}$, and classrooms are nested within schools, there is inevitably also strong evidence against H$_{\rm N}$ when tested against H$_{\rm S}$.
IM tests can be computed when testing against H$_{\rm S}$, but only for the model with school fixed effects. For both coefficients, the IM tests suggest that H$_{\rm N}$ should not be rejected. This is inconsistent with the results of the score\kern 0.04167em-variance tests and surprising in view of the standard errors reported in the top part of the table; see below for further discussion.
The results for testing H$_{\rm R}$ against H$_{\rm S}$ differ depending on the model, the coefficient(s) of interest, and the testing procedure. Consider first the model with no fixed effects. Here, both $\tau_\sigma$ statistics are negative, so of course upper-tail tests do not reject the null. This reflects the fact that, for both coefficients, the CR standard errors for school clustering are smaller than those for classroom clustering. The $\tau_\Sigma$ test for both coefficients jointly is always two\kern 0.04167em-sided. With $P$ values of 0.157 (asymptotic) and 0.171 (bootstrap), it also fails to reject the null hypothesis. Thus we conclude that the classroom level is the right one at which to cluster for the model with just a constant term.
Consider next the model with school fixed effects. As we noted in (ref), the “correct” level of clustering may be different for different hypotheses. This is what we find here. For $\hat\beta_s$, all three SV tests reject the null hypotheses and consequently suggest that school clustering is appropriate. In contrast, for $\hat\beta_a$, the SV tests suggest quite clearly (at least when using bootstrap $P$-values) that classroom clustering is appropriate.
Closer examination reveals that, for the model with school fixed effects, the asymptotic and bootstrap tests for H$_{\rm R}$ against H$_{\rm S}$ always yield quite different $P$ values. This is easily seen for $\beta_a$, where the bootstrap $P$ value of 0.344 is more than ten times the asymptotic $P$ value of 0.031. But it is also true for the other two tests. For $\beta_{\rm s}$, the $\tau_\sigma$ test statistic of 4.366 has an asymptotic $P$ value of 0.000006 and a bootstrap $P$ value of 0.0044. For the joint test of both coefficients, the $\tau_\Sigma$ test statistic of 28.673 has an asymptotic $P$ value of 0.000003 and a bootstrap $P$ value of 0.0109. In the latter two cases, the bootstrap $P$ values are small, but they are many times larger than the asymptotic ones.
The differences between asymptotic and bootstrap $P$ values for SV tests of classroom against school clustering in the model with school fixed effects arise because there are only a few classrooms per school. The average is 4.4, and most schools have just 3 or 4 classrooms. Because the residuals are orthogonal to the school fixed effects, they must add to zero over all classrooms in each school. This mechanically creates negative correlation between the residuals across classrooms within each school, even if the disturbances are uncorrelated across classrooms. The negative correlation of the residuals leads to spurious correlation of the empirical scores whenever a regressor of interest, after being projected off the fixed effects and the other regressors, is correlated across classrooms within schools. Because student characteristics probably vary at the school level, this sort of correlation seems very likely.
In principle, the spurious correlation of the empirical scores could be either positive or negative. For the model \hyperref[{eq:starcs}]{\tagform@{\ref*{eq:starcs}}}, it is evidently positive and quite large. This explains why the bootstrap tests yield much larger $P$ values than the asymptotic tests. Equivalently, the bootstrap critical values are greater than the asymptotic ones. For example, the test statistic for H$_{\rm N}$ against H$_{\rm R}$ for $\beta_a$ is 7.625. The asymptotic critical value for an upper-tail test at the 0.05 level is 1.645, but the bootstrap critical value is 3.423.
Whenever there is a dummy variable that affects only a few clusters (in this case the classrooms within each school), OLS residuals will be negatively correlated across those clusters, even when the disturbances are uncorrelated. This distortion of the residuals can cause cluster-robust inference to be severely misleading; see, among others, MW-JAE,MW-EJ and Chaisemartin_2022. However, CV$_3$ standard errors are almost certainly much more reliable in such cases than CV$_1$ standard errors. As MNW-bootknife explains, the cluster jackknife implicitly involves transforming the empirical scores in a way that undoes at least part of the distortion induced by least squares. This is evidently happening here.
With 330 clusters, we would normally expect CV$_1$ and CV$_3$ standard errors to be almost identical. But this is not the case for the model with fixed effects and classroom clustering. The CV$_1$ standard errors with classroom clustering for $\hat\beta_s$ and $\hat\beta_a$ are 2.322 and 2.109, respectively. These are much smaller than the CV$_3$ standard errors of 3.028 and 2.814 reported in (ref). The latter are almost certainly much more reliable than the former. Note that the CV$_1$ standard error for $\hat\beta_a$ with school clustering is 2.422, which is almost identical to the CV$_3$ one in the table and greater than 2.109. Thus the ratio of the S and R standard errors is greater than one for CV$_1$ and less than one for CV$_3$. Because the former ratio is greater than one, the $\tau_\sigma$ statistic is positive.
In additional simulation experiments not reported here, we generated artificial samples using the actual regressors for the STAR model. When there are no school fixed effects, all the SV tests, both asymptotic and bootstrap, work very well. However, when there are fixed effects, the asymptotic tests over-reject severely (up to about 70% of the time). The bootstrap tests perform almost perfectly when testing H$_{\rm N}$ against either H$_{\rm R}$ or H$_{\rm S}$, but they reject between 7% and 9% of the time for the tests of H$_{\rm R}$ against H$_{\rm S}$. We also performed some experiments in which the number of classrooms per school was doubled. All tests performed very much better in this case. These results suggest that, when there are fixed effects at the coarse level with few fine clusters per coarse cluster, and the asymptotic and bootstrap $P$ values differ sharply, the former should not be believed, and the latter should be taken with a grain of salt.
The IM tests are undoubtedly also affected by the odd properties of OLS residuals with school fixed effects. However, many of the differences between the score\kern 0.04167em-variance tests and the IM tests in (ref) probably arise because calculating the latter for the model \hyperref[{eq:starcs}]{\tagform@{\ref*{eq:starcs}}} is tricky. The problem is that estimating all the coefficients for every one of the 75 schools is infeasible. For 34 schools, it is impossible to estimate at least one of $\beta_s$ and $\beta_a$ (17 schools in the case of $\beta_s$ and 21 schools in the case of $\beta_a$). This means that the IM tests have to be based on either 58 or 54 coarse clusters, instead of all 75. Additionally, the other regressors that are included vary across clusters, so that the coefficients $\beta_s$ and $\beta_a$ may have different interpretations for different clusters. The IM tests may effectively be testing different null hypotheses than the score\kern 0.04167em-variance tests, which are always based on estimates for the entire sample.
In summary, our score\kern 0.04167em-variance tests suggest that clustering at either the classroom or school level is essential, because the null hypothesis of no clustering is always strongly rejected against both alternatives. Which of these levels we should cluster at depends on the model and the coefficient(s) of interest. With just a constant term, the sequential testing procedure, using either asymptotic or bootstrap tests, suggests that we should choose H$_{\rm R}$ and cluster at the classroom level. However, with school fixed effects, we should apparently choose H$_{\rm R}$ if interest focuses on $\beta_a$ and H$_{\rm S}$ if it focuses on $\beta_s$ or on both coefficients. Both choices lead us to conclude that the effect of small classes is positive and significant at the 0.05 level, while the effect of a teacher's aide is also positive but not significant at that level.
The fact that we obtain different results for the three SV tests should not be surprising. The test statistics depend on empirical scores, and they are different for the three tests because the ${\bm{Z}}$ matrices in \hyperref[{newmodel}]{\tagform@{\ref*{newmodel}}}, which are vectors for the $\tau_\sigma$ tests, are different; see (ref). For the model with fixed effects, the residuals are clearly correlated at the school level. While part of this correlation is evidently spurious and caused by the fixed effects, the bootstrap results suggest that the disturbances are surely correlated at the school level, because the $\tau_\sigma$ test for $\beta_s$ and the $\tau_\Sigma$ test for the two coefficients both reject quite strongly. For $\beta_a$ by itself, however, the scores are apparently not correlated, leading the $\tau_\sigma$ test not to reject in that case.
Empirical research that uses cluster-robust inference typically assumes that the level of clustering is known. When it is unknown, the consequences can be serious. Clustering at too fine a level can result in tests that over-reject severely and confidence intervals that under-cover dramatically. However, clustering at too coarse a level can lead to loss of power and to confidence intervals that vary greatly in length across samples and are, on average, excessively long.
We have proposed two direct tests for the level of clustering in a linear regression model, which we call score\kern 0.04167em-variance (or SV) tests. Both tests are based on the variances of the scores for two nested levels of clustering, because it is these variances that appear in the “filling” of the sandwich covariance matrices that correspond to the two levels. Under the null hypothesis that the finer level is appropriate, many of these variances are zero. The test statistics are functions of the empirical counterparts of those variances. Tests based on them can be used either to test the null of no clustering against an alternative of clustering at a certain level or to test the null of “fine” clustering against an alternative of “coarser” clustering. We have also proposed a sequential procedure which can be used to determine the correct level of clustering without inflating the family-wise error rate; see (ref).
The simplest of our two tests is based on the statistic $\tau_\sigma$. It has the form of a $t$-statistic and tests whether the variance of a particular coefficient estimate is the same for two different levels of clustering. It will be attractive whenever interest focuses on a single coefficient, and it can be implemented as either a one\kern 0.04167em-sided, upper-tail test or as a two\kern 0.04167em-sided test. Since upper-tail $\tau_\sigma$ tests have more power than two\kern 0.04167em-sided ones ((ref)), we believe that they will usually be the procedure of choice. The second variant, based on the Wald-like statistic $\tau_\Sigma$, tests whether the covariance matrix of a vector of coefficient estimates is the same for two different levels of clustering. It is necessarily two\kern 0.04167em-sided.
Our tests can be implemented as either asymptotic tests or as wild bootstrap tests. In (ref) and (ref), we derive the asymptotic distribution of our tests, prove that they are consistent tests, and also prove the validity of the wild bootstrap implementations. In the simulation experiments of (ref), the asymptotic tests often work well for tests of a single coefficient, but they can be seriously over-sized for tests of several coefficients. The problem is most severe when testing a moderate number of fine clusters against a small number of coarse clusters. For the empirical example of (ref), where several regressors, including the key ones, vary only at the fine\kern 0.04167em-cluster level, the asymptotic tests seem to be quite over-sized when there are school fixed effects. When the asymptotic tests are seriously over-sized, the bootstrap tests always perform much better.
Our score\kern 0.04167em-variance tests are very different from the other tests for the correct level of clustering proposed in Ibragimov_2016 and Cai_2022; see (ref). All these tests may provide valuable information, although we believe that SV tests are particularly intuitive. As we discuss in (ref), SV tests can be used either as formal pre\kern 0.04167em-tests for choosing the level at which to cluster or simply as robustness checks.
Both our simulation results and the empirical example suggest that SV tests can have excellent power. In many cases, with both actual and simulated data, the value of the test statistic is so far beyond any reasonable critical value that we can reject the null hypothesis with something very close to certainty even without bothering to use the bootstrap. However, when our tests are used as pre\kern 0.04167em-tests to choose the level of clustering, they inevitably make some Type I errors when the true clustering level is fine, and they inevitably make some Type II errors when the true clustering level is coarse but the sample size and the extent of coarse clustering are not large enough for rejection to occur all the time; see (ref).
The score\kern 0.04167em-variance tests we have proposed are intended to provide guidance for applied researchers. In our view, it should be routine to report the results of SV tests whenever more than one level of clustering is plausible. This is especially important when investigators are considering the use of heteroskedasticity-robust standard errors or clustering at a very fine level, such as by individual or by family. In practice, however, it may be safest to report inferences based on more than one level of clustering, along with the outcomes of SV tests, as we did in (ref).