The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
140,923 characters
Testing for the appropriate level of clustering in linear regression models
\title{Testing for the appropriate level of clustering\\in linear
regression models\thanks{We are grateful to the editor, Serena Ng, an
anonymous associate editor, three anonymous referees, Yevgeniy Feyman,
and participants at the 2018 Canadian Economics Association
conference, 2018 Canadian Econometric Study Group conference, 2020
Econometric Society World Congress, 2021 WEA conference, 2021 APPAM
conference, 2021 IAAE Annual Conference, Universit\'e de Montr\'eal,
University of Exeter, UCLA, UC Santa Barbara, Lakehead University,
Michigan State University, and Copenhagen Business School for
comments. MacKinnon and Webb thank the Social Sciences and Humanities
Research Council of Canada (SSHRC grants 435-2016-0871 and
435-2021-0396) for financial support. Nielsen thanks the Danish
National Research Foundation for financial support (DNRF Chair grant
number DNRF154). Computer code, including a \texttt{Stata} ado file,
for performing the testing procedures proposed here may be found at
\url{http://qed.econ.queensu.ca/pub/faculty/mackinnon/svtest/}.
}}
\author{James G. MacKinnon\thanks{Corresponding author. Address:
Department of Economics, 94 University Avenue, Queen's University,
Kingston, Ontario K7L 3N6, Canada. Email:\ \texttt{[email removed]}.
Tel.\ 613-533-2293. Fax 613-533-6668.}\\Queen's University\\
\texttt{[email removed]} \and
Morten \O rregaard Nielsen\\Aarhus University\\
\texttt{[email removed]} \and
Matthew D. Webb\\Carleton University\\
\texttt{[email removed]}}
\maketitle
\begin{abstract}
The overwhelming majority of empirical research that uses
cluster-robust inference assumes that the clustering structure
is known, even though there are often several possible ways in which a
dataset could be clustered. We propose two tests for the correct level
of clustering in regression models. One test focuses on inference about
a single coefficient, and the other on inference about two or more
coefficients. We provide both asymptotic and wild bootstrap
implementations. The proposed tests work for a null hypothesis of
either no clustering or ``fine'' clustering against alternatives of
``coarser'' clustering. We also propose a sequential testing procedure
to determine the appropriate level of clustering. Simulations suggest
that the bootstrap tests perform very well under the null hypothesis
and can have excellent power. An empirical example suggests that using
the tests leads to sensible inferences.
\medskip \noindent \textbf{Keywords:} CRVE, grouped data, clustered
data, cluster-robust variance estimator, robust inference, wild
bootstrap, wild cluster bootstrap.
\medskip \noindent \textbf{JEL Codes:} C12, C15, C21, C23.
\end{abstract}
\clearpage
\onehalfspacing
\section{Introduction}
\label{sec:intro}
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. \citet{Esarey_2019} and \citet{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,
\citet{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 \citep{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
\citep{MW-JAE,MW-EJ}, and can also reduce power \citep{MW-survey}.
In \Cref{sec:tests}, 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 \Cref{subsec:level}, 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 \citet{Ibragimov_2016} and recently by \citet{Cai_2022}. These
tests are very different from our tests and very different from each
other. We discuss them briefly in \Cref{subsec:other}.
The model of interest is discussed in \Cref{sec:model}. Our
score\kern 0.04167em-variance tests are described in \Cref{sec:tests}, 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. \Cref{sec:theory} provides asymptotic theory for the two
test statistics, the bootstrap tests, and the sequential testing
procedure. In \Cref{sec:partial}, 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
\Cref{sec:simulations}. An empirical example that deals with
clustering by classroom or school using the STAR dataset
\citep{Finn_1990,Mosteller_1995} is discussed in \Cref{sec:example}.
Finally, \Cref{sec:conclusion} concludes and offers some guidance for
empirical researchers. All mathematical proofs are given in
\Cref{sec:proofs}.
\section{The Regression Model with Clustering}
\label{sec:model}
We focus on the linear regression model
\begin{equation}
\label{model}
{\bm{y}} = {\bm{X}} {\bm\beta} + {\bm{u}},
\end{equation}
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
\begin{equation}
\label{betahat}
\hat{\bm\beta} = ({\bm{X}}^\top{\bm{X}})^{-1}{\bm{X}}^\top{\bm{y}}
= {\bm\beta}_0 + ({\bm{X}}^\top{\bm{X}})^{-1}{\bm{X}}^\top{\bm{u}}
= {\bm\beta}_0 + ({\bm{X}}^\top{\bm{X}})^{-1}\sum_{g=1}^G {\bm{X}}_g^\top{\bm{u}}_g,
\end{equation}
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
\begin{equation}
\label{def Sigma}
{\rm E}({\bm{s}}_g{\bm{s}}_{g'}^\top)= {\mathbb I} (g=g') {\bm{\Sigma}}_g, \quad g,g'=1,\ldots,G,
\end{equation}
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}]{\textup{\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
\begin{equation}
\label{covbeta}
\widehat\operatorname{Var} (\hat{\bm\beta} ) =
({\bm{X}}^\top{\bm{X}})^{-1} \hat{\bm{\Sigma}} ({\bm{X}}^\top{\bm{X}})^{-1}\kern -.08333em,
\end{equation}
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}]{\textup{\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}]{\textup{\tagform@{\ref*{def Sigma}}}} is
\begin{equation}
\label{Sighat}
\hat{\bm{\Sigma}}_{\rm c} = m_{\rm c} \sum_{g=1}^G {\bm{X}}_g^\top\hat{\bm{u}}_g
\hat{\bm{u}}_g^\top{\bm{X}}_g = m_{\rm c} \sum_{g=1}^G \hat{\bm{s}}_g\hat{\bm{s}}^\top_g,
\end{equation}
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}]{\textup{\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 \citet{DMN_2019} and
\citet{HansenLee_2019}.
\begin{remark}
\label{rem:hetmodel}
In the special case in which each cluster has $N_g =1$ observation, we
can use
\begin{equation}
\label{HCmid}
\hat{\bm{\Sigma}}_{\rm het} = \sum_{i=1}^N \hat{u}^2_i\kern 0.08333em {\bm{X}}_i^\top{\bm{X}}_i
= {\bm{X}}^\top\kern -.08333em \operatorname{diag} (\hat{u}^2_1,\dots,\hat{u}^2_N) {\bm{X}}\kern -.08333em,
\end{equation}
where ${\bm{X}}_i$ is the \th{i} row of the ${\bm{X}}$ matrix and $\hat{u}_i$
is the \th{i} residual. The variance matrix obtained by setting
$\hat{\bm{\Sigma}} = \hat{\bm{\Sigma}}_{\rm het}$ in \hyperref[{covbeta}]{\textup{\tagform@{\ref*{covbeta}}}} is the famous
heteroskedasticity-consistent variance matrix estimator (HCCME) of
\citet{Eicker_1963} and \citet{White_1980}. Of course, the matrix
$\hat{\bm{\Sigma}}_{\rm het}$ can be modified in various ways to improve its
finite\kern 0.04167em-sample properties \citep{MW_1985,JGM_2013}. The simplest is
to multiply it by \mbox{$m_{\rm het}=N/(N-k)$}, so that
\hyperref[{covbeta}]{\textup{\tagform@{\ref*{covbeta}}}} becomes what is usually called HC$_1$.
\end{remark}
\begin{remark}
\label{rem:AAIW} As \citet*{AAIW_2023} point out, when the object of
interest is the average treatment effect in a finite population,
cluster-robust standard errors based on \hyperref[{Sighat}]{\textup{\tagform@{\ref*{Sighat}}}} can be
``unnecessarily conservative.'' Consequently, they develop an approach
to inference that depends both on how the data were sampled and on how
treatment was assigned. In this paper, however, we follow most of the
literature on cluster-robust inference and rely on the traditional
approach in which every sample is treated as a random outcome from a
data-generating process (DGP). The objective is to draw inferences
about the parameters of the DGP, which may be interpreted as features
of an infinitely large population; see \citet{MNW-guide} for
additional details.
\end{remark}
\section{The Testing Procedure}
\label{sec:tests}
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}]{\textup{\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 \citet{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}]{\textup{\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$,
\begin{equation}
\label{def Sigma gh}
{\rm E} ({\bm{s}}_{gh}{\bm{s}}_{gh'}^\top ) = {\mathbb I} (h=h') {\bm{\Sigma}}_{gh},
\quad h,h' = 1,\ldots ,M_g ,
\end{equation}
where each of the ${\bm{\Sigma}}_{gh}$ is a $k \times k$ matrix. Thus,
\hyperref[{def Sigma}]{\textup{\tagform@{\ref*{def Sigma}}}} and \hyperref[{def Sigma gh}]{\textup{\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 \Cref{rem:two-way}.
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}]{\textup{\tagform@{\ref*{def Sigma}}}} and \hyperref[{def Sigma gh}]{\textup{\tagform@{\ref*{def Sigma gh}}}}, these
matrices are
\begin{equation}
\label{Sigmas}
{\bm{\Sigma}}_{\rm c} = \sum_{g=1}^G {\bm{\Sigma}}_g
\quad\textrm{and}\quad
{\bm{\Sigma}}_{\rm f} = \sum_{g=1}^G \sum_{h=1}^{M_g}{\bm{\Sigma}}_{gh}.
\end{equation}
We consider the null and alternative hypotheses
\begin{equation}
\label{hypotheses}
\H{0}\!: \underset{N \to \infty}{\lim} \kern 0.08333em {\bm{\Sigma}}_{\rm f} \kern 0.08333em {\bm{\Sigma}}^{-1}_{\rm c} = {\bf I}
\quad \textrm{and} \quad
\H{1}\!: \underset{N \to \infty}{\lim} \kern 0.08333em {\bm{\Sigma}}_{\rm f} \kern 0.08333em {\bm{\Sigma}}^{-1}_{\rm c} \neq {\bf I}.
\end{equation}
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 \citet{DMN_2019}.
\begin{remark}
\label{rem:hyp}
In \hyperref[{hypotheses}]{\textup{\tagform@{\ref*{hypotheses}}}}, we are not directly testing the fine clustering
condition in \hyperref[{def Sigma gh}]{\textup{\tagform@{\ref*{def Sigma gh}}}}. Instead, we are testing an
important implication of the clustering structure. Specifically, we
test whether ${\bm{\Sigma}}_{\rm c}={\bm{\Sigma}}_{\rm f}$, which implies that
a valid CRVE for $\hat{\bm\beta}$ is given by \hyperref[{covbeta}]{\textup{\tagform@{\ref*{covbeta}}}} with
$\hat{\bm{\Sigma}} = \hat{\bm{\Sigma}}_{\rm f}$.
\end{remark}
\begin{remark}
\label{rem:hetnull}
An important null hypothesis is that no CRVE is needed because the
HCCME considered in \Cref{rem:hetmodel}, obtained by combining
\hyperref[{covbeta}]{\textup{\tagform@{\ref*{covbeta}}}} and \hyperref[{HCmid}]{\textup{\tagform@{\ref*{HCmid}}}}, is valid. In this case, each fine
cluster has just one observation, so that $M_g=N_g$ and $N_{gh}=1$ for
all $g$ and~$h$.
\end{remark}
\begin{remark}
\label{rem:partial}
In practical applications, the number of coefficients in regression
models, and hence the size of the CRVE matrices, is often large, so
that comparing these matrices directly can be impractical.
Furthermore, it is usually only one coefficient, or a small subset of
them, that is actually of interest. Many coefficients typically
correspond to fixed effects and other conditioning variables that are
not of primary interest. By partialing out the latter, it is possible
to reduce the dimensionality of the test and focus on the parameter(s)
of interest. We discuss this issue in \Cref{sec:partial}.
\end{remark}
\subsection{Test Statistics}
\label{sec:stats}
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 (\Cref{sec:partial}), 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}]{\textup{\tagform@{\ref*{Sighat}}}} is
\begin{equation}
\label{Sigmac}
\hat{\bm{\Sigma}}_{\rm c}
= m_{\rm c} \sum_{g=1}^G \hat{\bm{s}}_g\hat{\bm{s}}^\top_g
= m_{\rm c} \sum_{g=1}^G \left( \sum_{h=1}^{M_g}
\hat{\bm{s}}_{gh} \right)\!
\left( \sum_{h=1}^{M_g} \hat{\bm{s}}_{gh} \right)^{\!\!\!\top}\!.
\end{equation}
Similarly, we can write, c.f.\ \hyperref[{def Sigma gh}]{\textup{\tagform@{\ref*{def Sigma gh}}}} and \hyperref[{Sigmas}]{\textup{\tagform@{\ref*{Sigmas}}}},
\begin{equation}
\label{Sigmaf}
\hat{\bm{\Sigma}}_{\rm f} = m_{\rm f} \sum_{g=1}^G \sum_{h=1}^{M_g}
\hat{\bm{s}}_{gh} \hat{\bm{s}}_{gh}^\top ,
\end{equation}
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}]{\textup{\tagform@{\ref*{Sigmac}}}} and \hyperref[{Sigmaf}]{\textup{\tagform@{\ref*{Sigmaf}}}}
reduce to the scalars
\begin{equation}
\label{Sigmascalar}
\hat\sigma^2_{\rm c} = m_{\rm c} \sum_{g=1}^G
\left( \sum_{h=1}^{M_g} \hat s_{gh} \right)^{\!\!2}
\quad\textrm{and}\quad
\hat\sigma^2_{\rm f} = m_{\rm f} \sum_{g=1}^G \sum_{h=1}^{M_g}
\hat s_{gh}^2 .
\end{equation}
The quantities given in \hyperref[{Sigmac}]{\textup{\tagform@{\ref*{Sigmac}}}}, \hyperref[{Sigmaf}]{\textup{\tagform@{\ref*{Sigmaf}}}}, and
\hyperref[{Sigmascalar}]{\textup{\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
\begin{equation}
\label{sigmahet}
\hat\sigma^2_{\rm het} =
\sum_{g=1}^G\sum_{h=1}^{M_g}\sum_{i=1}^{N_{gh}} \hat s^2_{ghi},
\end{equation}
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}]{\textup{\tagform@{\ref*{Sigmascalar}}}}, namely,
\begin{equation}
\label{thetascalar}
\hat\theta = \hat\sigma^2_{\rm c} - \hat\sigma^2_{\rm f}.
\end{equation}
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,
\begin{equation}
\label{thetaSigma}
\hat{\bm{\theta}} = \operatorname{vech} (\hat{\bm{\Sigma}}_{\rm c} - \hat{\bm{\Sigma}}_{\rm f}),
\end{equation}
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}]{\textup{\tagform@{\ref*{thetascalar}}}} and
\hyperref[{thetaSigma}]{\textup{\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
\begin{equation}
\label{thetaknown}
\theta = \sum_{g=1}^G\sum_{h_1=1}^{M_g}\sum_{h_2\neq h_1}^{M_g}
s_{gh_1} s_{gh_2}.
\end{equation}
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}]{\textup{\tagform@{\ref*{def Sigma gh}}}}.
The variance of $\theta$ in \hyperref[{thetaknown}]{\textup{\tagform@{\ref*{thetaknown}}}} is, under the null
hypothesis,
\begin{equation}
\label{var1}
\operatorname{Var} (\theta ) = \sum_{g=1}^G\sum_{h_1=1}^{M_g} \sum_{\ell_1=1}^{M_g}
\sum_{h_2\neq h_1}^{M_g}\sum_{\ell_2\neq \ell_1}^{M_g}
{\rm E} ( s_{gh_1}s_{g\ell_1}s_{gh_2}s_{g\ell_2}).
\end{equation}
The expectation of any product of scores can only be nonzero, under
the null and \hyperref[{def Sigma gh}]{\textup{\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}]{\textup{\tagform@{\ref*{var1}}}} simplifies to
\begin{equation}
\label{tauvar}
\operatorname{Var}(\theta) = 2 \sum_{g=1}^G \sum_{h_1=1}^{M_g}\sum_{h_2\neq h_1}^{M_g}
\sigma^2_{gh_1}\sigma^2_{gh_2} ,
\end{equation}
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}]{\textup{\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}]{\textup{\tagform@{\ref*{def Sigma gh}}}}. This suggests the variance
estimator
\begin{equation}
\label{varfast}
\widehat\operatorname{Var} (\hat\theta) =
2 \sum_{g=1}^G\left(\sum_{h=1}^{M_g} \hat s^2_{gh}\!\right)^{\!\!2}
- 2 \sum_{g=1}^G \sum_{h=1}^{M_g} \hat s^4_{gh}.
\end{equation}
This equation avoids the triple summation in \hyperref[{tauvar}]{\textup{\tagform@{\ref*{tauvar}}}} by
squaring the sums of squared empirical scores, which then requires
that the second term be subtracted. In deriving \hyperref[{varfast}]{\textup{\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}]{\textup{\tagform@{\ref*{thetascalar}}}} and \hyperref[{varfast}]{\textup{\tagform@{\ref*{varfast}}}} yields the studentized
test statistic
\begin{equation}
\label{eq:taus}
\tau_\sigma = \frac{\hat\theta}
{\sqrt{\widehat\operatorname{Var} (\hat\theta})}.
\end{equation}
In \Cref{sec:theory}, we show that $\tau_\sigma$ is asymptotically
distributed as ${\rm N}(0,1)$.
\begin{remark}
\label{rem:onesided}
The statistic defined in \hyperref[{eq:taus}]{\textup{\tagform@{\ref*{eq:taus}}}} yields either a one\kern 0.04167em-sided
or a two\kern 0.04167em-sided test. Upper-tail tests may often be of primary
interest, because we expect the diagonal elements of ${\bm{\Sigma}}_{\rm c}$
to exceed the corresponding elements of ${\bm{\Sigma}}_{\rm f}$ when there
is positive correlation within clusters under the alternative.
However, since this is not necessarily the case, two\kern 0.04167em-sided tests
based on $\tau_\sigma^2$ may also be of interest. The asymptotic
theory in \Cref{sec:theory} handles both cases.
\end{remark}
\begin{remark}
\label{rem:Moulton}
Consider again the special case in which the null is
heteroskedasticity with no clustering. When the elements of ${\bm{x}}$
display little intra-cluster correlation, the contrast $\hat\theta$,
and hence the absolute value of $\tau_\sigma$, will tend to be small,
even if the residuals display a great deal of intra-cluster
correlation. This is what we should expect, because in that case the
so\kern 0.04167em-called Moulton factor, the ratio of clustered to non-clustered
standard errors \citep{Moulton_1986}, will be relatively small. Of
course, the opposite will be true when the elements of ${\bm{x}}$ display
a lot of intra-cluster correlation. Thus, all else equal, SV tests may
well yield different results for different choices of~${\bm{x}}$.
\end{remark}
When $k>1$, so that $\hat{\bm{\theta}}$ is a vector, the variance estimator
analogous to \hyperref[{varfast}]{\textup{\tagform@{\ref*{varfast}}}} is
\begin{equation}
\label{var2sided}
\widehat\operatorname{Var}( \hat{\bm{\theta}} ) =
2\sum_{g=1}^G {\bm{H}}_k \bigg(\sum_{h=1}^{M_g} \hat{\bm{s}}_{gh}\hat{\bm{s}}_{gh}^\top
\otimes\sum_{h=1}^{M_g}\hat{\bm{s}}_{gh}\hat{\bm{s}}_{gh}^\top\!\bigg) {\bm{H}}_k^\top
-2\sum_{g=1}^G \sum_{h=1}^{M_g}{\bm{H}}_k \big(\hat{\bm{s}}_{gh}\hat{\bm{s}}_{gh}^\top
\otimes \hat{\bm{s}}_{gh}\hat{\bm{s}}_{gh}^\top\big) {\bm{H}}_k^\top .
\end{equation}
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}$ \citep[p.~354]{Harville_1997}, and $\otimes$ denotes
the Kronecker product. A studentized (Wald) statistic is then given by
\begin{equation}
\label{eq:tauGs}
\tau_\Sigma =\hat{\bm{\theta}}^\top
\widehat\operatorname{Var}( \hat{\bm{\theta}} )^{-1} \hat{\bm{\theta}} .
\end{equation}
In \Cref{sec:theory}, we show that $\tau_\Sigma$ is asymptotically
distributed as $\chi^2(k(k+1)/2)$.
\begin{remark}
\label{rem:Phillips}
As pointed out by a referee, \citet{CP_2018} develop simple measures
of the discrepancy between two positive definite symmetric matrices. When
bootstrapped, these measures can be used as alternative test
statistics for cases where $k>1$. Preliminary simulations suggest that
these bootstrap tests can work well, although not (in general) better
than our proposed bootstrap tests based on \hyperref[{eq:tauGs}]{\textup{\tagform@{\ref*{eq:tauGs}}}}. A full
analysis is beyond the scope of this paper and is therefore left for
future work.
\end{remark}
\begin{remark}
\label{rem:two-way} It is possible to use the test statistics
\hyperref[{eq:taus}]{\textup{\tagform@{\ref*{eq:taus}}}} and \hyperref[{eq:tauGs}]{\textup{\tagform@{\ref*{eq:tauGs}}}} for testing one\kern 0.04167em-way against
two\kern 0.04167em-way clustering. Suppose there are two alternative clustering
dimensions, labeled A and~B, and their intersection is labeled~I.
These could correspond to, say, state~(A) and year~(B), with the
intersection denoting observations that correspond to the same year in
the same state. Let $\hat{\bm{\Sigma}}_j$ denote the (one\kern 0.04167em-way) CRVE in
\hyperref[{Sighat}]{\textup{\tagform@{\ref*{Sighat}}}} under clustering dimension $j \in \{
\rm{A},\rm{B},\rm{I} \}$. Then the two\kern 0.04167em-way CRVE \citep*{CGM_2011}
is given by \hyperref[{covbeta}]{\textup{\tagform@{\ref*{covbeta}}}} with
\begin{equation}
\label{eq:twowayCRVE}
\hat{\bm{\Sigma}}_{\kern 0.04167em\rm TW} = \hat{\bm{\Sigma}}_{\rm A} + \hat{\bm{\Sigma}}_{\rm B}
- \hat{\bm{\Sigma}}_{\kern 0.04167em\rm I}.
\end{equation}
If we test the null of one\kern 0.04167em-way clustering by~A against the
alternative of clustering by both A and~B, then $\hat{\bm{\Sigma}}_{\rm c}=
\hat{\bm{\Sigma}}_{\kern 0.04167em\rm TW}$ and $\hat{\bm{\Sigma}}_{\rm f}=\hat{\bm{\Sigma}}_{\rm A}$.
Therefore, the vector of contrasts in \hyperref[{thetaSigma}]{\textup{\tagform@{\ref*{thetaSigma}}}} becomes
\begin{equation}
\label{eq:twowaytheta}
\hat{\bm{\theta}} = \operatorname{vech} ( \hat{\bm{\Sigma}}_{\kern 0.04167em\rm TW} - \hat{\bm{\Sigma}}_{\rm A} )
= \operatorname{vech} ( \hat{\bm{\Sigma}}_{\rm B} - \hat{\bm{\Sigma}}_{\kern 0.04167em\rm I} ).
\end{equation}
The result in \hyperref[{eq:twowaytheta}]{\textup{\tagform@{\ref*{eq:twowaytheta}}}} shows that testing the null of
one\kern 0.04167em-way clustering by~A against the alternative of two\kern 0.04167em-way
clustering by~A and~B must lead to the same test statistic as testing
the null of one\kern 0.04167em-way clustering by~I against the alternative of
one\kern 0.04167em-way clustering by~B.
Although it is straightforward to derive a test statistic based on
\hyperref[{eq:twowaytheta}]{\textup{\tagform@{\ref*{eq:twowaytheta}}}}, the asymptotic analysis of this statistic
would be different from the analysis for testing nested one\kern 0.04167em-way
clustering in \Cref{sec:theory} below. For example, to derive the
asymptotic null distribution of a statistic based on
\hyperref[{eq:twowaytheta}]{\textup{\tagform@{\ref*{eq:twowaytheta}}}} for testing clustering by~I against clustering
by~B would mean analyzing it under the DGP that clustering is in fact
by~A. Therefore, we leave this analysis for future work.
Similarly, if there is two\kern 0.04167em-way clustering under both the null and
alternative hypotheses, it may be feasible to calculate
score\kern 0.04167em-variance statistics similar to \hyperref[{eq:taus}]{\textup{\tagform@{\ref*{eq:taus}}}} and
\hyperref[{eq:tauGs}]{\textup{\tagform@{\ref*{eq:tauGs}}}}. However, the analysis of the asymptotic null
distribution would require completely different, and technically
nontrivial, techniques
\citep*{DDG_2021,MNW_multi,Menzel_2021,Chiang_2022}.
\end{remark}
In this section, we have proposed two score\kern 0.04167em-variance tests of
\hyperref[{hypotheses}]{\textup{\tagform@{\ref*{hypotheses}}}}. They both involve comparing different variance
estimates of the empirical scores, namely, the two scalars in
\hyperref[{Sigmascalar}]{\textup{\tagform@{\ref*{Sigmascalar}}}} for the $\tau_\sigma$ test and the matrices in
\hyperref[{Sigmac}]{\textup{\tagform@{\ref*{Sigmac}}}} and \hyperref[{Sigmaf}]{\textup{\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 \Cref{sec:partial}. 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 \Cref{sec:simulations}, 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}]{\textup{\tagform@{\ref*{eq:taus}}}} and~\hyperref[{eq:tauGs}]{\textup{\tagform@{\ref*{eq:tauGs}}}},
which often perform much better in finite samples, as we also show in
\Cref{sec:simulations}. These bootstrap implementations are described
next.
\subsection{Bootstrap Implementation}
\label{subsec:bootstrap}
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 \citep{Wu_1986, Liu_1988} or the wild
cluster bootstrap \citep{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
\Cref{subsec:boot}.
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 \Cref{rem:hetnull}, 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
\citet{DF_2008}, \citet{DMN_2019}, and \citet{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$.
\begin{algorithm}[Bootstrap test implementation]
\label{alg:BS}
Let $B >\!> 1$ denote the number of bootstrap replications, and let
$\tau$ denote the chosen test statistic.
\begin{enumerate}
\item Estimate model \hyperref[{model}]{\textup{\tagform@{\ref*{model}}}} by OLS to obtain the
residuals~$\hat{\bm{u}}$.
\item Compute the empirical score vector $\hat{\bm{s}}$ and use it to
compute~$\tau$.
\item For $b=1,\ldots,B$,
\begin{enumerate}
\item generate the vector of bootstrap dependent variables ${\bm{y}}^{*b}
= {\bm{u}}^{*b}$ from the residual vector $\hat{\bm{u}}$ using the wild
cluster bootstrap corresponding to the null hypothesis, or the
ordinary wild bootstrap if the null does not involve clustering.
\item Regress ${\bm{y}}^{*b}$ on ${\bm{X}}$ to obtain the bootstrap residuals
$\hat{\bm{u}}^{*b}$, and use these, together with ${\bm{X}}$, to compute
$\tau^{*b}$, the bootstrap analog of~$\tau$.
\end{enumerate}
\item Compute the bootstrap $P$~value $\hat{P}^\ast =B^{-1}\sum_{b=1}^B
{\mathbb I}(\tau^{*b} > \tau)$.
\end{enumerate}
\end{algorithm}
As usual, if $\alpha$ is the level of the test, then $B$ should be
chosen so that $(1-\alpha)(B+1)$ is an integer \citep{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 \citep{DM_2000}.
\begin{remark}
\label{rem:BSonesided}
When $\tau$ is defined as $\tau_\sigma$, \Cref{alg:BS} yields a
one\kern 0.04167em-sided upper-tail test. When $\tau$ is defined as $|\tau_\sigma|$,
$\tau_\sigma^2$, or $\tau_\Sigma$, it yields a two\kern 0.04167em-sided test; see
\Cref{rem:onesided}.
\end{remark}
\begin{remark}
\label{rem:critvals}
If desired, bootstrap critical values can be calculated as quantiles
of the $\tau^{*b}$. For example, when $B=999$ and the $\tau^{*b}$ are
sorted from smallest to largest, the 0.05 critical value for a
one\kern 0.04167em-sided upper-tail test is $\tau^{*b'}$ for $b' = (1-0.05)(B+1) =
950$.
\end{remark}
\begin{remark}
\label{rem:ordinary}
We could use the ordinary wild bootstrap instead of the wild cluster
bootstrap in \Cref{alg:BS}, even when the null hypothesis involves
clustering. The same intuition as in \citet{DMN_2019} applies, whereby
the ordinary wild bootstrap would lead to asymptotically valid tests
because the statistics are asymptotically pivotal. There may be cases,
like the ones considered in \citet{MW-EJ} and/or ones in which the
number of fine clusters is small, in which the wild bootstrap would
perform better than the wild cluster bootstrap. However, we believe
that such cases are likely to be rare.
\end{remark}
\subsection{Choosing the Level of Clustering by Sequential Testing}
\label{subsec:level}
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
\citep[e.g.,][]{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.
\citet{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}]{\textup{\tagform@{\ref*{covbeta}}}} and \hyperref[{Sigmas}]{\textup{\tagform@{\ref*{Sigmas}}}}.
Here we assume that ${\bm{\Sigma}}_0$ corresponds to no clustering, c.f.\
\Cref{rem:hetmodel,rem:hetnull}, 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.
\begin{algorithm}[Nested sequential testing procedure]
\label{alg:seq}
Let $m=0$. Then:
\begin{enumerate}
\item Test $\H{0}\!: \underset{N \to \infty}{\lim} \kern 0.08333em {\bm{\Sigma}}_m {\bm{\Sigma}}_{m+1}^{-1} = {\bf I}$
against $\H{1}\!: \underset{N \to \infty}{\lim} \kern 0.08333em {\bm{\Sigma}}_m {\bm{\Sigma}}_{m+1}^{-1} \neq {\bf I}$.
\item If the test in step~1 does not reject, choose $\hat{m}=m$ and
stop.
\item If $m=p-1$ and the test in step~1 rejects, choose $\hat{m}=p$
and stop.
\item If $m \leq p-2$ and the test in step~1 rejects, increment
$m$ by $1$ and go to step~1.
\end{enumerate}
\end{algorithm}
We can equivalently state the sequential testing problem in
\Cref{alg:seq} as a type of estimation problem. Specifically,
\begin{equation}
\label{mhat}
\hat{m} = \min \{ m \in \{ 0,1,\dots,p \} \textrm{ such that }
\H{0}\!: \underset{N \to \infty}{\textrm{plim}} \kern 0.08333em {\bm{\Sigma}}_m {\bm{\Sigma}}_{m+1}^{-1} = {\bf I}
\textrm{ is not rejected} \}.
\end{equation}
Of course, the $\hat{m}$ resulting from \Cref{alg:seq} and from
\hyperref[{mhat}]{\textup{\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 \Cref{alg:seq} or
\hyperref[{mhat}]{\textup{\tagform@{\ref*{mhat}}}} is (nearly) consistent. Precise asymptotic
properties of the proposed sequential procedure are established in
\Cref{subsec:seq}, and finite\kern 0.04167em-sample performance is investigated
by Monte Carlo simulation methods in \Cref{subsec:seqsims}.
\subsection{Other Tests for the Level of Clustering}
\label{subsec:other}
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 \Cref{subsec:asy} requires that $G\to\infty$.
However, see the discussion in \Cref{subsec:fixedG}.
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 \Cref{sec:example}. \citet{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 \citet{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.
\citet{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
\citet{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.
\subsection{Inference about Regression Coefficients}
\label{subsec:infreg}
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}]{\textup{\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 \citep[e.g.,][]{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 \Cref{subsec:pretest}, we report the results of
some simulation experiments that compare confidence intervals based on
several alternative procedures.
As \citet{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 \Cref{subsec:other}, 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 \citep{Cai_2022}.
\section{Asymptotic Theory}
\label{sec:theory}
In \Cref{subsec:asy}, 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 (\Cref{subsec:boot}) and
prove asymptotic results for the sequential testing procedure
(\Cref{subsec:seq}). We first state and discuss the assumptions needed
for our proofs, which may be found in \Cref{sec:proofs}.
\begin{assumption}
\label{as:cluster}
The sequence ${\bm{s}}_{gh}=\sum_{i=1}^{N_{gh}}{\bm{X}}_{ghi}^\top u_{ghi}$ is
independent across both $g$ and~$h$.
\end{assumption}
\begin{assumption}
\label{as:moments}
For all $g,h$, it holds that ${\rm E} ({\bm{s}}_{gh}) = {\bm{0}}$ and $\operatorname{Var}
({\bm{s}}_{gh}) ={\bm{\Sigma}}_{gh}$. Furthermore, $\sup_{g,h,i}{\rm E} \Vert{\bm{s}}_{ghi}
\Vert^{2\lambda}<\infty$ for some $\lambda >1$.
\end{assumption}
\begin{assumption}
\label{as:X}
The regressor matrix ${\bm{X}}$ satisfies $\sup_{g,h,i} {\rm E} \Vert {\bm{X}}_{ghi}
\Vert^2 <\infty$ and $N^{-1}{\bm{X}}^\top{\bm{X}} \overset{P} \longrightarrow {\bm{\Xi}}$, where ${\bm{\Xi}}$ is
finite and positive definite.
\end{assumption}
\begin{assumption}
\label{as:eigen}
Let $\omega_{\min}(\cdot)$ and $\omega_{\max}(\cdot)$ denote the
minimum and maximum eigenvalues of the argument. Then
$\inf_{g,h}N_{gh}^{-1}\omega_{\min} ({\bm{\Sigma}}_{gh} )>0$ and
$\sup_{g,h}\omega_{\max} \big({\bm{\Sigma}}_{gh} (\sum_{h=1}^{M_g}
{\bm{\Sigma}}_{gh})^{-1} \big)<1$.
\end{assumption}
\begin{assumption}
\label{as:size}
For $\lambda$ defined in \Cref{as:moments}, the cluster sizes satisfy
\par \vspace{\abovedisplayskip}
\hfill $\displaystyle \frac{\sup_g N_g^2 \sup_{g,h} N_{gh}^2}
{\sum_{g=1}^G \omega_{\min}\big(\sum_{h=1}^{M_g}{\bm{\Sigma}}_{gh}\big)^2}
\longrightarrow 0
\quad \textrm{and} \quad
\frac{N^{1/\lambda}\sup_g N_g \sup_{g,h}N_{gh}^{3-1/\lambda}}
{\sum_{g=1}^G \omega_{\min}\big(\sum_{h=1}^{M_g}{\bm{\Sigma}}_{gh}\big)^2}
\longrightarrow 0.$
\end{assumption}
\Cref{as:cluster} is the assumption of (at most) ``fine'' clustering,
which implies that the null hypothesis in \hyperref[{hypotheses}]{\textup{\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
\Cref{as:moments} and the multicollinearity condition in \Cref{as:X}
are standard in linear regression models.
Next, the conditions in \Cref{as:eigen} 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 \Cref{as:size} 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 \Cref{as:size} 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 \Cref{as:size}, which
allows larger clusters without dominating the limit; a similar
tradeoff was found in \citet{DMN_2019}. Furthermore, more homogeneity
in cluster sizes allows for fewer and larger clusters. We illustrate
these tradeoffs in the following remarks.
\begin{remark}
\label{rem:tradeoff}
Under \Cref{as:cluster,as:eigen}, both denominators in \Cref{as:size}
are bounded from below by $\sum_{g=1}^G N_g^2 \geq c N \inf_g N_g$, and
a sufficient condition for \Cref{as:size} is
\begin{equation}
\label{suff1a}
\sup_{g,h} N_{gh}^2 \left( \frac{\sup_g N_g}{\inf_g N_g} \right)
\left( \frac{\sup_g N_g}{N} \right) \longrightarrow 0
\quad\textrm{and}\quad
\left( \frac{\sup_g N_g}{\inf_g N_g} \right)^{\!\lambda}
\left( \frac{\sup_{g,h}N_{gh}^{3\lambda-1}}{N^{\lambda-1}} \right)
\longrightarrow 0.
\end{equation}
If the cluster sizes are bounded under the alternative, i.e.\ $\sup_g
N_g < \infty$, then \hyperref[{suff1a}]{\textup{\tagform@{\ref*{suff1a}}}} is easily satisfied. Note that
$\sup_g N_g /N \to 0$, and hence $G\to\infty$, is implied by
\Cref{as:size}, and it is therefore not stated explicitly. Suppose,
on the other hand, that \Cref{as:eigen} were strengthened to assume that
$\inf_{g,h} N_{gh}^{-2} \omega_{\min}({\bm{\Sigma}}_{gh})>0$, as would be the
case if a random-effects or factor-type model were assumed under the
null. In that case, \Cref{as:size} and \hyperref[{suff1a}]{\textup{\tagform@{\ref*{suff1a}}}} could be weakened
substantially.
\end{remark}
\begin{remark}
\label{rem:sizes}
It is interesting to consider a setup for clusters that are relatively
homogeneous, but possibly unbounded, in size. To make this concrete,
suppose the coarse clusters have size $N_g =O( N^\alpha )$ for
$g=1,\ldots,G$, where $\alpha \in [0,1)$ and `$O (\cdot )$' is to be
understood as an exact rate subject to $N_g$ being an integer. Because
$\sum_{g=1}^G N_g = N$, it then holds that $G = O ( N^{1-\alpha} )$.
Similarly, for each $g$, the fine clusters have size $N_{gh} = O (
N_g^\gamma )$ for $h=1,\ldots,M_g =O ( N_g^{1-\gamma})$. That is, when
$\alpha$ is large (small), there are few large (many small) coarse
clusters. Similarly, when $\gamma$ is large (small), there are few large
(many small) fine clusters per coarse cluster. Under this setup,
\hyperref[{suff1a}]{\textup{\tagform@{\ref*{suff1a}}}} is satisfied if $\alpha (2\gamma+1)<1$ and $\alpha\gamma
<(\lambda-1)/(3\lambda-1)$.
The important implication of this setup is that the implied
restrictions on the cluster sizes in \Cref{as:size} are very weak. In
fact, if we assume that the fine cluster sizes are bounded (i.e.,
$\gamma = 0$), which applies, for example, in the important special
case in which the scores are independent but heteroskedastic under the
null, then we can allow $G=O(N^{1-\alpha})$ for any $\alpha <1$. That
is, the number of coarse clusters can be arbitrarily close to $O(1)$.
For example, we allow $G=O(N^{0.1})$ and $N_g = O(N^{0.9})$, which
corresponds to very few and very large coarse clusters. In this sense,
our asymptotic framework can nearly accommodate the fixed-$G$ setup;
see \Cref{subsec:fixedG}.
\end{remark}
\subsection{Theory for Asymptotic Tests}
\label{subsec:asy}
\begin{theorem}
\label{thm:asy}
Let \Cref{as:cluster,as:moments,as:eigen,as:X,as:size} be satisfied.
Then, as $N \to \infty$, it holds that
\begin{align*}
\operatorname{Var} ({\bm{\theta}} )^{-1/2} \hat{\bm{\theta}} &\overset{d} \longrightarrow {\rm N} (0,{\bf I}), &
\operatorname{Var} ({\bm{\theta}} )^{-1}\widehat\operatorname{Var} (\hat{\bm{\theta}} ) &\overset{P} \longrightarrow {\bf I},
\quad\textrm{and}\\
\frac{\hat\theta}{\sqrt{\operatorname{Var}(\theta)}} &\overset{d} \longrightarrow {\rm N} (0,1), &
\frac{\widehat\operatorname{Var} (\hat\theta)}{\operatorname{Var} (\theta)} &\overset{P} \longrightarrow 1.
\end{align*}
\end{theorem}
\begin{remark}
\label{rem:self-norm}
Observe that the statement of the asymptotic distributions in
\Cref{thm:asy} only concerns quantities that are self-normalized. For
example, in the scalar case, these are either $\hat\theta$ divided
by its true standard error or the estimated variance of $\hat\theta$
divided by the true variance. This is because the appropriate rates of
convergence are not known in general; see the discussion
below~\hyperref[{hypotheses}]{\textup{\tagform@{\ref*{hypotheses}}}}.
\end{remark}
The asymptotic distributions of the test statistics follow
immediately from \Cref{thm:asy}.
\begin{corollary}
\label{cor:asy}
Let \Cref{as:cluster,as:moments,as:eigen,as:X,as:size} be satisfied.
Then, as $N \to \infty$, it holds that
\begin{equation*}
\tau_\Sigma \overset{d} \longrightarrow \chi^2 \big(k(k+1)/2\big)
\quad\textrm{and}\quad
\tau_\sigma \overset{d} \longrightarrow {\rm N}(0,1).
\end{equation*}
\end{corollary}
We next consider the asymptotic behavior of the test statistics under
the alternative. Because \Cref{as:cluster} implies that \H{0} is true,
we do not make that assumption. Instead, we impose the following
conditions:
\begin{assumption}
\label{as:clusteralt}
The sequence ${\bm{s}}_g = {\bm{X}}_g^\top {\bm{u}}_g = \sum_{h=1}^{N_g}{\bm{s}}_{gh}$
is independent across~$g$.
\end{assumption}
\goodbreak
\begin{assumption}
\label{as:sizealt}
The cluster sizes satisfy
\par \vspace{\abovedisplayskip}
\hfill $\displaystyle \frac{\sup_g N_g^{3/2}N^{1/2}}
{\sum_{g=1}^G \omega_{\min}({\bm{\Sigma}}_g)} \longrightarrow 0.$
\end{assumption}
\goodbreak
\Cref{as:clusteralt} 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
\Cref{as:sizealt}.
\begin{remark}
\label{rem:tradeoffalt}
As in \Cref{rem:tradeoff}, there is a tradeoff between cluster size
heterogeneity and intra-cluster correlation, in this case correlation
within coarse clusters. Specifically, under \Cref{as:eigen}, the
denominator in \Cref{as:sizealt} is bounded from below by $\sum_{g=1}^G
N_g = N$\kern -.08333em, and hence a sufficient condition for \Cref{as:sizealt} is
\begin{equation}
\label{sizealt}
\frac{\sup_g N_g^3}{N}\longrightarrow 0.
\end{equation}
Suppose instead that \Cref{as:eigen} were strengthened to assume that
$\inf_g N_g^{-2}\omega_{\min}({\bm{\Sigma}}_g)>0$ (as in
\Cref{rem:tradeoff}, this could be due to a random-effects model or a
factor-type model). That is, more correlation is assumed within the
coarse clusters, so that there is a stronger departure from the null
hypothesis. In this case, the denominator in \Cref{as:sizealt} is
bounded from below by $\sum_{g=1}^G N_g^2 \geq \inf_g N_g N$.
Therefore, a sufficient condition for \Cref{as:sizealt} is
\begin{equation}
\label{sizealt2}
\frac{\sup_g N_g^3}{\inf_g N_g^2 N}\longrightarrow 0.
\end{equation}
With relatively homogeneous coarse clusters as in \Cref{rem:sizes},
i.e.\ coarse clusters where $\sup_g N_g$ and $\inf_g N_g$ are of the
same order of magnitude, the condition \hyperref[{sizealt2}]{\textup{\tagform@{\ref*{sizealt2}}}} reduces to
$\sup_g N_g /N \to 0$, which is clearly minimal and implied by
\Cref{as:size}.
\end{remark}
\begin{theorem}
\label{thm:cons}
Let \Cref{as:moments,as:X,as:eigen,as:size,as:clusteralt,as:sizealt} be
satisfied, and suppose \H{0} in \hyperref[{hypotheses}]{\textup{\tagform@{\ref*{hypotheses}}}} is not true. Then,
as $N \to \infty$, it holds that
\begin{equation*}
\tau_\Sigma \overset{P} \longrightarrow +\infty \quad\textrm{and}\quad |\tau_\sigma| \overset{P} \longrightarrow +\infty.
\end{equation*}
\end{theorem}
It follows immediately from \Cref{thm:cons} 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}]{\textup{\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.
\subsection{Theory for Bootstrap Tests}
\label{subsec:boot}
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.
\begin{theorem}
\label{thm:boot}
Let \Cref{as:moments,as:X,as:eigen,as:size,as:clusteralt} be satisfied
with $\lambda \geq 2$, and assume that ${\rm E}^\ast |v^\ast|^{2\lambda}
<\infty$. Then, as $N \to \infty$, it holds for any $\epsilon >0$ that
\begin{equation*}
P \big( \sup_{x \in \mathbb R} \big| P^\ast (\tau^\ast \leq x) -
P_0 (\tau \leq x) \big| > \epsilon \big) \longrightarrow 0 .
\end{equation*}
\end{theorem}
First, note that the bootstrap theory requires a slight strengthening of
the moment condition since at least four moments are now required.
Second, \Cref{thm:boot} shows that the bootstrap $P$~values in
\Cref{alg:BS} are asymptotically valid under \Cref{as:cluster} and \H{0}.
Third, note that neither the null hypothesis nor \Cref{as:cluster} is
imposed in \Cref{thm:boot}. Thus \Cref{thm:asy,thm:cons,thm:boot}
together show immediately that the bootstrap tests are consistent. We
summarize these results in the following corollary.
\begin{corollary}
\label{cor:boot}
Let \Cref{as:moments,as:X,as:eigen,as:size} be satisfied with $\lambda
\geq 2$, and assume that ${\rm E}^\ast |v^\ast|^{2\lambda} < \infty$. As
$N \to \infty$, it holds that:
\begin{itemize}
\item[(i)] If \Cref{as:cluster} is satisfied and \H{0} is true, then
$\hat P^\ast \overset{d} \longrightarrow {\rm U}(0,1)$, where ${\rm U}(0,1)$ is a uniform random
variable on~$[0,1]$.
\item[(ii)] If \Cref{as:clusteralt,as:sizealt} are satisfied and \H{0}
is not true, then $\hat P^\ast \overset{P} \longrightarrow 0$.
\end{itemize}
\end{corollary}
\subsection{Theory for Sequential Testing Procedure}
\label{subsec:seq}
The next theorem provides theoretical justification for the sequential
testing procedure given in \Cref{alg:seq}.
\begin{theorem}
\label{thm:seq}
Let $\hat m$ be defined in \Cref{alg:seq} or \hyperref[{mhat}]{\textup{\tagform@{\ref*{mhat}}}}. Suppose
\Cref{as:cluster} is satisfied when the ``fine'' clustering level in
\hyperref[{hypotheses}]{\textup{\tagform@{\ref*{hypotheses}}}} is $m = m_0 \in \{ 0,1,\ldots ,p \}$ (and hence also
when $m>m_0$), and suppose \H{0} in \hyperref[{hypotheses}]{\textup{\tagform@{\ref*{hypotheses}}}} is not true for
clustering levels $m<m_0$. Suppose also that
\Cref{as:moments,as:X,as:eigen,as:size,as:sizealt} are satisfied, and
let $\alpha$ denote the nominal level of the tests. As $N \to \infty$, it holds
that
\begin{itemize}
\item[(i)] if $m_0 \leq p-1$, then $P(\hat m \leq m_0-1) \to 0$,
$P(\hat m =m_0 ) \to 1-\alpha$, and $P(\hat m \geq m_0+1) \to \alpha$;
\item[(ii)] if $m_0 = p$, then $P ( \hat m \leq m_0-1 ) \to 0$ and
$P(\hat m = m_0 ) \to 1$.
\end{itemize}
\end{theorem}
The results in \Cref{thm:seq} show that $\hat m$ defined in
\Cref{alg:seq} or \hyperref[{mhat}]{\textup{\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.
\subsection{Fixed-$G$ Asymptotic Theory}
\label{subsec:fixedG}
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
\citet{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 \citet{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}]{\textup{\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
\begin{equation}
\label{chi square}
\frac{\theta}{\sqrt{\operatorname{Var} (\theta)}}\overset{d} \longrightarrow
\sum_{g=1}^G w_g^2 (\chi_{1,g}^2 -1),
\end{equation}
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}]{\textup{\tagform@{\ref*{chi square}}}}.
The limiting distribution in \hyperref[{chi square}]{\textup{\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}]{\textup{\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 \Cref{thm:asy} and \Cref{cor:asy} 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.
\section{Dimension Reduction by Partialing Out}
\label{sec:partial}
As discussed in \Cref{rem:partial}, 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 \Cref{sec:tests} 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}]{\textup{\tagform@{\ref*{model}}}} can then be rewritten as
\begin{equation}
\label{newmodel}
{\bm{y}} = {\bm{Z}}{\bm\beta}_1 + {\bm{X}}_2{\bm{\delta}} + {\bm{u}}, \qquad
{\bm{Z}} = {\bm{M}}_{{\bm{X}}_2}{\bm{X}}_1,
\end{equation}
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}]{\textup{\tagform@{\ref*{model}}}} and \hyperref[{newmodel}]{\textup{\tagform@{\ref*{newmodel}}}} have exactly the same explanatory
power and the same disturbances, ${\bm{u}}$. The coefficient ${\bm\beta}_1$ in
\hyperref[{newmodel}]{\textup{\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}]{\textup{\tagform@{\ref*{betahat}}}},
\begin{equation}
\label{beta1hat}
\hat{\bm\beta}_1 = ({\bm{Z}}^\top{\bm{Z}})^{-1} {\bm{Z}}^\top{\bm{y}}
= {\bm\beta}_{1,0} + ({\bm{Z}}^\top{\bm{Z}})^{-1}\sum_{g=1}^G {\bm{Z}}_g^\top{\bm{u}}_g,
\end{equation}
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
\begin{equation}
\label{Z and Q}
{\bm{Z}} = {\bm{X}} {\bm{Q}} \quad\textrm{with}\quad
{\bm{Q}} =[ {\bf I}_{k_1} , \; -{\bm{X}}_1^\top{\bm{X}}_2 ({\bm{X}}_2^\top{\bm{X}}_2 )^{-1}]^\top .
\end{equation}
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}]{\textup{\tagform@{\ref*{beta1hat}}}} and
\hyperref[{Z and Q}]{\textup{\tagform@{\ref*{Z and Q}}}}, we obtain the following sandwich formula, c.f.\
\hyperref[{covbeta}]{\textup{\tagform@{\ref*{covbeta}}}},
\begin{equation}
\label{covbeta1}
\widehat\operatorname{Var} (\hat{\bm\beta}_1 ) =
({\bm{Z}}^\top{\bm{Z}})^{-1} {\bm{Q}}^\top \hat{\bm{\Sigma}}\kern 0.08333em {\bm{Q}} ({\bm{Z}}^\top{\bm{Z}})^{-1}\kern -.08333em .
\end{equation}
Under \Cref{as:X}, ${\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}]{\textup{\tagform@{\ref*{covbeta1}}}} is
clearly an estimator of ${\bm{A}}^\top {\bm{\Sigma}} {\bm{A}}$.
The matrix ${\bm{Q}}$ in \hyperref[{Z and Q}]{\textup{\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}]{\textup{\tagform@{\ref*{newmodel}}}}, we consequently replace the hypotheses in
\hyperref[{hypotheses}]{\textup{\tagform@{\ref*{hypotheses}}}} with
\begin{equation}
\label{newhypotheses}
\H{0}\!: \underset{N \to \infty}{\lim} \kern 0.08333em ({\bm{A}}^\top {\bm{\Sigma}}_{\rm f}{\bm{A}} ) \kern 0.08333em ({\bm{A}}^\top
{\bm{\Sigma}}_{\rm c} {\bm{A}} )^{-1} = {\bf I}
\quad \textrm{and} \quad
\H{1}\!: \underset{N \to \infty}{\lim} \kern 0.08333em ( {\bm{A}}^\top {\bm{\Sigma}}_{\rm f}{\bm{A}} ) \kern 0.08333em ({\bm{A}}^\top
{\bm{\Sigma}}_{\rm c} {\bm{A}} )^{-1} \neq {\bf I}.
\end{equation}
Furthermore, from \hyperref[{Z and Q}]{\textup{\tagform@{\ref*{Z and Q}}}}, we see that we can use the same
algebra for the model in \hyperref[{newmodel}]{\textup{\tagform@{\ref*{newmodel}}}} as for the model in
\hyperref[{model}]{\textup{\tagform@{\ref*{model}}}} to define the test statistics, i.e.\ \hyperref[{Sigmac}]{\textup{\tagform@{\ref*{Sigmac}}}},
\hyperref[{Sigmaf}]{\textup{\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 \Cref{alg:BS}. Of course, degrees-of-freedom
corrections like the factor $m_c$ in \hyperref[{Sighat}]{\textup{\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$.
\begin{remark}
\label{rem:empscores}
The empirical scores ${\bm{Q}}^\top\hat s_{gh}$ and
${\bm{Q}}^\top\hat{\bm{s}}_{gh}$ depend on the matrix ${\bm{Z}} =
{\bm{M}}_{{\bm{X}}_2}{\bm{X}}_1$, which is the residual matrix from regressing
${\bm{X}}_1$ on ${\bm{X}}_2$. Therefore, different choices for ${\bm{X}}_1$ will
yield different empirical scores, and hence different test statistics;
see \Cref{rem:Moulton}. This is also reflected in the hypotheses in
\hyperref[{newhypotheses}]{\textup{\tagform@{\ref*{newhypotheses}}}}, where different choices for ${\bm{X}}_1$ will yield
a different ${\bm{A}}$ matrix and hence different null and alternative
hypotheses.
\end{remark}
\begin{remark}
\Cref{thm:asy,thm:cons,thm:boot,thm:seq} continue to hold with the new
definitions given in this section, with ${\bm\beta}_1$ replacing ${\bm\beta}$,
${\bm{Z}}$ replacing ${\bm{X}}$, and $k_1$ replacing $k$. Because the matrix
${\bm{Q}} \overset{P} \longrightarrow {\bm{A}}$ under \Cref{as:X}, it acts only as a fixed constant
in all asymptotic arguments; that is, ${\bm{Q}}^\top{\bm{s}}_{gh} = {\bm{A}}^\top
{\bm{s}}_{gh}(1+o_P(1))$. Thus, the same proofs apply with ${\bm{s}}_{gh}$
replaced by ${\bm{A}}^\top{\bm{s}}_{gh}$.
\end{remark}
\begin{remark}
\label{rem:sufficient}
Careful inspection of the proofs shows that, in the setup of this section,
we can replace ${\bm{\Sigma}}_{gh}$ with ${\bm{A}}^\top{\bm{\Sigma}}_{gh}{\bm{A}}$ in
\Cref{as:size}. This could be attractive in some cases. Suppose, for
example, that ${\bm{X}}_1$ and ${\bm{X}}_2$ are (asymptotically) orthogonal,
such that ${\bm{A}}^\top{\bm{\Sigma}}{\bm{A}}$ is equal to the diagonal block of
${\bm{\Sigma}}$ corresponding to ${\bm{X}}_1^\top{\bm{u}}$. Suppose also that ${\bm{X}}_1$
and ${\bm{u}}$ are both finely clustered, but the ${\bm{X}}_2$ are independent.
Then ${\bm{A}}^\top{\bm{\Sigma}}_{gh}{\bm{A}}$ satisfies the condition in
\Cref{rem:tradeoff}, while ${\bm{\Sigma}}_{gh}$ only satisfies the
corresponding condition in \Cref{as:eigen}, and hence using
${\bm{A}}^\top{\bm{\Sigma}}_{gh}{\bm{A}}$ in \Cref{as:size} would lead to a weaker
condition.
\end{remark}
\section{Simulation Experiments}
\label{sec:simulations}
Most of the papers cited in the second paragraph of \Cref{sec:intro}
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 \emph{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, \citet{BDM_2004},
\citet{JGM-CJE}, and \Cref{sec:example}. 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
\begin{equation}
\label{facDGP}
\begin{aligned}
z_{gi} &= \rho^{1/2}\kern 0.08333em\xi^1_g + (1-\rho)^{1/2}\kern 0.04167em\zeta_{gi}
\;\;\mbox{ if $i$ is odd}\cr
z_{gi} &= \rho^{1/2}\kern 0.08333em\xi^2_g + (1-\rho)^{1/2}\kern 0.04167em\zeta_{gi}
\;\;\mbox{ if $i$ is even.}\cr
\end{aligned}
\end{equation}
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}]{\textup{\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}]{\textup{\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}]{\textup{\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}]{\textup{\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 \Cref{sec:partial}) are generated using \hyperref[{facDGP}]{\textup{\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.
\subsection{Performance under the Null Hypothesis}
\label{subsec:null}
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
\begin{equation}
\label{simmod}
y_{ghi} = \sum_{\ell=1}^{k_1} \beta_\ell X^\ell_{ghi} + {\bm{X}}^2_{gh}{\bm{\delta}}
+ u_{ghi},
\end{equation}
where the regressors $X^\ell_{ghi}$ are generated independently across
$\ell$ by \hyperref[{facDGP}]{\textup{\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}$.
\begin{figure}[tb]
\begin{center}
\caption{Rejection frequencies for $\tau_\Sigma$ tests at
0.05 level, $G_{\rm c}$ varying}
\label{fig:1}
\includegraphics[width=\textwidth]{ctfigv2A.pdf}
\end{center}
{\footnotesize
\textbf{Notes:} The regressors are generated by \hyperref[{facDGP}]{\textup{\tagform@{\ref*{facDGP}}}} with
$\rho=0.5$ and $1 \leq k_1 \leq 5$. The regressand is generated
by~\hyperref[{simmod}]{\textup{\tagform@{\ref*{simmod}}}}. The disturbances are independent standard normals in
Panels~(b) and~(d) and finely clustered with $\rho=0.1$ in Panels~(a)
and~(c). $G_{\rm c}$ denotes the number of coarse clusters. Each
coarse cluster contains 400 observations, so that $N=400G_{\rm c}$. In
Panels~(a) and~(c), there are $G_{\rm f} = 4G_{\rm c}$ fine clusters,
each containing 100 observations. Bootstrap tests employ $B=399$.
Panel~(c) uses the wild cluster bootstrap, and Panel~(d) uses the
ordinary wild bootstrap. There are 400,000 replications.}
\end{figure}
\Cref{fig:1} 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 \Cref{fig:1} 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 \citep{White_1982}, the $\tau_\Sigma$ test has degrees of
freedom that are $O(k_1^2)$. \citet{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}]{\textup{\tagform@{\ref*{Sigmac}}}}, \hyperref[{Sigmaf}]{\textup{\tagform@{\ref*{Sigmaf}}}}, and
\hyperref[{var2sided}]{\textup{\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.
\begin{figure}[tb]
\begin{center}
\caption{Rejection frequencies for $\tau_\Sigma$ tests at
0.05 level, $G_{\rm f}$ or $N_g$ varying}
\label{fig:2}
\includegraphics[width=\textwidth]{ctfigv2B.pdf}
\end{center}
{\footnotesize
\textbf{Notes:}
The regressors are generated as in \Cref{fig:1}, but only for $k_1=1$,
3, and~5. There are $G_{\rm c}=8$ coarse clusters. In Panels~(a)
and~(c), there are between 3 and 12 fine clusters per coarse cluster,
each with 100 observations. In Panels~(b) and~(d), there are just coarse
clusters, with between 25 and 400 observations per coarse cluster.
Bootstrap tests employ $B=399$. Panel~(c) uses the wild cluster
bootstrap, and Panel~(d) uses the ordinary wild bootstrap. There are
400,000 replications.}
\end{figure}
In \Cref{fig:2}, 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}]{\textup{\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 \Cref{fig:1}. 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.
\begin{figure}[tb]
\begin{center}
\caption{Rejection frequencies for $\tau_\sigma$ tests at
0.05 level, $G_{\rm c}$ varying}
\label{fig:3}
\includegraphics[width=\textwidth]{ctfigv2C.pdf}
\end{center}
{\footnotesize
\textbf{Notes:} There is one regressor, which is generated by
\hyperref[{facDGP}]{\textup{\tagform@{\ref*{facDGP}}}} with $\rho=0.5$. In Panels~(a) and~(c), the disturbances
are finely clustered with~$\rho=0.1$ and $G_{\rm f} = 4G_{\rm c}$ fine
clusters, each with 100 observations. In Panels~(b) and~(d), they are
independent standard normals. $G_{\rm c}$ denotes the number of coarse
clusters, each of which contains 400 observations, so that
$N=400G_{\rm c}$. Bootstrap tests employ $B=399$. Panel~(c) uses the
wild cluster bootstrap, and Panel~(d) uses the ordinary wild bootstrap.
There are 400,000 replications.}
\end{figure}
\Cref{fig:3} 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 \Cref{fig:1}, 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.
\subsection{The Power of Bootstrap Tests}
\label{subsec:alt}
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}]{\textup{\tagform@{\ref*{simmod}}}}, with one
regressor and coarse fixed effects. As usual, the regressor is
generated by \hyperref[{facDGP}]{\textup{\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
\citep{DM_2000}.
\begin{figure}[tb]
\begin{center}
\caption{Power of bootstrap $\tau_\sigma$ tests at 0.05 level when there
is coarse clustering}
\label{fig:4}
\includegraphics[width=\textwidth]{ctfigv2D.pdf}
\end{center}
{\footnotesize
\textbf{Notes:} The data are generated by \hyperref[{simmod}]{\textup{\tagform@{\ref*{simmod}}}} with coarse
fixed effects and coarse (or no) clustering. There are 5000
observations, 10 coarse clusters, 20, 40, or 100 fine clusters,
400,000 replications, and 999 bootstraps.}
\end{figure}
\Cref{fig:4} 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. \Cref{fig:5} 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
\begin{equation}
N_g = \left[\frac{N \exp(\delta g/G)}{\sum_{j=1}^G \exp(\delta
j/G)}\right]
\enspace
\hbox{for } g=1,\ldots,G-1,
\label{gamma-eq}
\end{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 \citet{MW-JAE}, \citet{DMN_2019}, and several other papers. In
the experiments of \Cref{fig:5}, $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 \Cref{fig:5}, 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}]{\textup{\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 \Cref{fig:5} 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.
\begin{figure}[tb]
\begin{center}
\caption{Cluster size variation and the performance of upper-tail
$\tau_\sigma$ tests}
\label{fig:5}
\includegraphics[width=\textwidth]{ctfigv2E.pdf}
\end{center}
{\footnotesize
\textbf{Notes:} The regressor is generated by \hyperref[{facDGP}]{\textup{\tagform@{\ref*{facDGP}}}} with
$\rho=0.5$. The disturbances are generated by \hyperref[{facDGP}]{\textup{\tagform@{\ref*{facDGP}}}} with
$\rho=0.0$ in Panel~(a) and $\rho=0.1$ in Panel~b). There are 10
clusters, 1000 observations, and cluster fixed effects. Cluster sizes
vary according to~\hyperref[{gamma-eq}]{\textup{\tagform@{\ref*{gamma-eq}}}}. All tests are at the nominal 0.05
level. There are 400,000 replications and 399 bootstraps.}
\end{figure}
\subsection{The Sequential Testing Procedure}
\label{subsec:seqsims}
Our next set of experiments concerns the sequential testing procedure
of \Cref{subsec:seq}, using bootstrap tests. These experiments are
quite similar to the ones in \Cref{fig:3,fig:4}, except that there are
8 coarse clusters, 48 fine clusters, and 2400 observations. As in
\Cref{subsec:alt}, 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}]{\textup{\tagform@{\ref*{facDGP}}}},
varies within either coarse or fine clusters.
\begin{figure}[tb]
\begin{center}
\caption{Outcomes for sequential upper-tail bootstrap tests at 0.05
level}
\label{fig:6}
\includegraphics[width=\textwidth]{ctfigv2F.pdf}
\end{center}
{\footnotesize
\textbf{Notes:} There is one regressor, generated by \hyperref[{facDGP}]{\textup{\tagform@{\ref*{facDGP}}}}
with $\rho=0.5$, plus cluster fixed effects. The regressand is
generated by \hyperref[{simmod}]{\textup{\tagform@{\ref*{simmod}}}} with clustered disturbances at either the
coarse level (left panel) or the fine level (right panel), for $\rho$
between 0.0 and~0.64. There are 8 coarse clusters, 48 fine clusters,
and 2400 observations. Bootstrap tests use $B=999$, and there are
400,000 replications. The solid red and dashed purple curves separate
the three outcomes of the sequential procedure; the red curve
separates N from~F, and the purple curve F from~C. The dashed blue
curve shows the outcome of a direct test of N against~C.}
\end{figure}
In Panel~(a) of \Cref{fig:6}, 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 \Cref{fig:4}. 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 \Cref{fig:6}, 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.
\subsection{Making Inferences about a Regression Coefficient}
\label{subsec:pretest}
In \Cref{subsec:infreg}, 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}]{\textup{\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}]{\textup{\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}]{\textup{\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}]{\textup{\tagform@{\ref*{covbeta}}}} with
$\hat{\bm{\Sigma}}_{\rm c}$ given by \hyperref[{Sighat}]{\textup{\tagform@{\ref*{Sighat}}}} or \hyperref[{Sigmac}]{\textup{\tagform@{\ref*{Sigmac}}}}.
However, recent work \citep{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 \citet{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}]{\textup{\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 \Cref{fig:4}. 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.
\begin{figure}[tb]
\begin{center}
\caption{Root mean squared errors of four standard error estimates}
\label{fig:7}
\includegraphics[width=\textwidth]{ctfigv2H.pdf}
\end{center}
{\footnotesize \textbf{Notes:} The regressors are generated by
\hyperref[{facDGP}]{\textup{\tagform@{\ref*{facDGP}}}} with coarse clustering and $\rho=0.5$, and the
disturbances are generated as discussed in the second paragraph of
this subsection. When $\eta=0$, there is either no clustering (top
panels) or fine clustering (bottom panels), depending on the
pre\kern 0.04167em-test null hypothesis. When $\eta>0$, there is coarse clustering.
The pre\kern 0.04167em-test estimators are based on upper-tail $\tau_\sigma$
tests. There are $400,\kern -.08333em000$ replications.}
\end{figure}
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. \Cref{fig:7} 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 \Cref{fig:7}
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.
\begin{figure}[tb]
\begin{center}
\caption{Coverage of four 95\% confidence intervals}
\label{fig:8}
\includegraphics[width=\textwidth]{ctfigv2G.pdf}
\end{center}
{\footnotesize \textbf{Notes:} These results are for the same experiments
as in \Cref{fig:7}.}
\end{figure}
Of course, using a more accurate standard error does not guarantee
better coverage. \Cref{fig:8} shows the coverage of confidence intervals
using the four standard errors in \Cref{fig:7}. 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
\citep{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$.
\section{Empirical Example}
\label{sec:example}
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
\citep{Finn_1990, Mosteller_1995}. We use these data to estimate a
cross\kern 0.04167em-sectional model similar to one in \citet{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:
\begin{equation}
\label{eq:starcs}
\text{read-one}_{sri} = \alpha + \beta_s\kern 0.08333em \text{small-class}_{sr}
+\beta_a\kern 0.08333em \text{aide\kern 0.04167em-class}_{sr} + {\bm{x}}^\top_{sri}{\bm{\delta}} + u_{sri}.
\end{equation}
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.
\begin{table}[tp]
\caption{STAR Example}
\label{tab:starcs}
\vspace*{-0.5em}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}
llcd{2.3}d{1.3}d{1.3}cd{2.3}d{1.3}d{1.3}d{1.3}}
\toprule
& & & \multicolumn{3}{c}{Without School FE}
& & & \multicolumn{3}{c}{With School FE} \\
\cmidrule{4-6}\cmidrule{9-11}
\multicolumn{2}{l}{Estimates} & & \multicolumn{1}{c}{HC$_3$(N)}
& \multicolumn{1}{c}{CV$_3$(R)}
& \multicolumn{1}{c}{CV$_3$(S)} &&
& \multicolumn{1}{c}{HC$_3$(N)} & \multicolumn{1}{c}{CV$_3$(R)} &
\multicolumn{1}{c}{CV$_3$(S)} \\
\midrule
small & $\hat\beta_s$ & & 9.211 & 9.211 & 9.211 & &&8.095 &8.095 & 8.095 \\
& s.e. & & 1.633 & 3.273 & 3.253 && &1.556 &3.028 & 3.120 \\
& $t$-stat. & & 5.640 & 2.814 & 2.831 && &5.203 &2.673 & 2.595 \\
aide & $\hat\beta_a$ & & 6.245 & 6.245 & 6.245 & &&4.170 & 4.170 &4.170 \\
& s.e. & & 1.664 & 3.343 & 2.847 & &&1.587 & 2.814 & 2.429 \\
& $t$-stat. & & 3.752 & 1.868 & 2.194 && &2.628 & 1.482 & 1.717 \\
\midrule
& & & \multicolumn{3}{c}{Without School FE}
& & \multicolumn{4}{c}{With School FE} \\
\cmidrule{4-6}\cmidrule{8-11}
\multicolumn{2}{l}{Cluster tests} &
& \multicolumn{1}{c}{SV stat.} &
\multicolumn{1}{c}{asy.\ $P$} & \multicolumn{1}{c}{boot $P$} & &
\multicolumn{1}{c}{SV stat.} & \multicolumn{1}{c}{asy.\ $P$} &
\multicolumn{1}{c}{boot $P$} & \multicolumn{1}{c}{IM $P$} \\
\midrule
small & H$_{\rm N}$ vs H$_{\rm R}$ & &28.388 &0.000 &0.000
& &12.757 &0.000 &0.000 & \multicolumn{1}{c}{---} \\
& H$_{\rm N}$ vs H$_{\rm S}$ & &16.409 &0.000 &0.000
& &18.308 &0.000 &0.000 &0.251 \\
& H$_{\rm R}$ vs H$_{\rm S}$ & & -0.101 &0.540 &0.543
& &4.366 &0.000 &0.004 &0.000 \\
aide & H$_{\rm N}$ vs H$_{\rm R}$ & &25.693 &0.000 &0.000
& &7.625 &0.000 &0.000 & \multicolumn{1}{c}{---} \\
& H$_{\rm N}$ vs H$_{\rm S}$ & &10.102 &0.000 &0.000
& &7.696 &0.000 &0.000 &0.438 \\
& H$_{\rm R}$ vs H$_{\rm S}$ & & -1.765 &0.961 &0.973
& &1.871 &0.031 &0.344 & 0.000 \\
both &H$_{\rm N}$ vs H$_{\rm R}$ & &1075.469 &0.000 &0.000
& &180.448 &0.000 &0.000 & \multicolumn{1}{c}{---} \\
&H$_{\rm N}$ vs H$_{\rm S}$ & &322.367 &0.000 &0.000
& &385.950 &0.000 &0.000 & \multicolumn{1}{c}{---} \\
&H$_{\rm R}$ vs H$_{\rm S}$ & &5.215 &0.157 &0.171
& &28.673 &0.000 &0.011 & \multicolumn{1}{c}{---} \\
\bottomrule
\end{tabular*}
\vskip 6pt {\footnotesize \textbf{Notes:} There are 3,989 observations
and either 330 classroom clusters (denoted R for ``room'') or 75 school
clusters (denoted~S). The null hypotheses of no clustering, classroom
clustering, and school clustering are called H$_{\rm N}$, H$_{\rm
R}$, and H$_{\rm S}$, respectively. Values of the $\tau_\sigma$
statistic (for ``small'' and ``aide'') or the $\tau_\Sigma$ statistic
(for ``both'') are shown under ``SV~stat.'' All other numbers in the
lower panel are $P$~values. For the $\tau_\sigma$ tests, asymptotic
$P$~values are upper-tail and based on the ${\rm N}(0,1)$ distribution. For
the $\tau_\Sigma$ tests, they are based on the $\chi^2(3)$
distribution. Upper-tail bootstrap $P$~values use $B=99,\kern -.08333em999$. IM
tests use $S=9,\kern -.08333em999$. Data and \texttt{Stata} files may be found at
\url{http://qed.econ.queensu.ca/pub/faculty/mackinnon/svtest/}.}
\end{table}
OLS estimates for the model \hyperref[{eq:starcs}]{\textup{\tagform@{\ref*{eq:starcs}}}} are presented in the top
half of \Cref{tab:starcs}. 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 \Cref{subsec:pretest}, 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 \Cref{tab:starcs} 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
\Cref{rem:Moulton,rem:empscores}, 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}]{\textup{\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, \citet{MW-JAE,MW-EJ} and \citet{Chaisemartin_2022}. However,
CV$_3$ standard errors are almost certainly much more reliable in such
cases than CV$_1$ standard errors. As \citet{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
\Cref{tab:starcs}. 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
\Cref{tab:starcs} probably arise because calculating the latter for
the model \hyperref[{eq:starcs}]{\textup{\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}]{\textup{\tagform@{\ref*{newmodel}}}}, which are vectors for the $\tau_\sigma$
tests, are different; see \Cref{rem:empscores}. 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.
\section{Conclusion}
\label{sec:conclusion}
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 \Cref{subsec:level}.
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
(\Cref{subsec:alt}), 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 \Cref{sec:theory} and \Cref{sec:proofs}, 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
\Cref{sec:simulations}, 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 \Cref{sec:example}, 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 \citet{Ibragimov_2016}
and \citet{Cai_2022}; see \Cref{subsec:other}. All these tests may
provide valuable information, although we believe that SV tests are
particularly intuitive. As we discuss in \Cref{subsec:infreg}, 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 \Cref{subsec:pretest}.
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 \Cref{sec:example}.