EconBase
← Back to paper

A Modified Randomization Test for the Level of Clustering

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

66,639 characters · 14 sections · 44 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

A Modified Randomization Test for the Level of Clustering

\def\spacingset#1{ {#1}} \spacingset{1}

\if00 \fi

\if10 {

center[center omitted — 96 chars of source]

} \fi

abstractSuppose a researcher observes individuals within a county within a state. Given concerns about correlation across individuals, it is common to group observations into clusters and conduct inference treating observations across clusters as roughly independent. However, a researcher that has chosen to cluster at the county level may be unsure of their decision, given knowledge that observations are independent across states. This paper proposes a modified randomization test as a robustness check for the chosen level of clustering in a linear regression setting. Existing tests require either the number of states or number of counties to be large. Our method is designed for settings with few states and few counties. While the method is conservative, it has competitive power in settings that may be relevant to empirical work.

{\it Keywords:} Linear Regression, Clustered Standard Errors, Small-Cluster Asymptotics

\spacingset{1.45}

Introduction

Consider the following regression:

equation*[equation* omitted — 73 chars of source]

where a researcher wants to perform inference on $\beta$. If the researcher is concerned about correlation between $U_{i}$ and $U_{i'}$, it is frequently helpful to group observations into independent clusters. These independent clusters can then be used to construct cluster-robust covariance estimators (CCE) as in lz1986, or for approximate randomization tests as in crs2017 and ccks2021.

However, these procedures require the assignment of units to clusters be known ex ante. In practice, researchers often have some freedom in choosing the level at which to cluster their standard errors. For example, those working with the American Community Survey (ACS) can cluster their data either at the individual, county or state level. Alternatively, those working with firm data from COMPUSTAT have the option to cluster firms at either the 4-digit, 3-digit or 2-digit Standard Industrial Classification (SIC) level.

Clustering at the correct level is important for valid inference. A large body of simulation evidence shows that ignoring cluster dependence -- in other words, clustering at too fine a level -- leads to type I errors that exceed the nominal error by as much as 10 times (bdm2004; cgm2008). On the other hand, clustering at excessively coarse levels can also lead to problems. For one, coarse clusters tend to be few in number. It is well-known that confidence intervals based on the cluster-robust standard errors tend to under-cover when the number of clusters is small (see mhe2008 for instance), leading to poor size control. In the absence of under-coverage issues, unnecessarily coarse levels of clustering can also lead to tests with poor power since the researcher assumes less information than they actually have. aaiw17 demonstrate via simulations, in a many-cluster setting, that CCEs based on coarse clusters can be too large. They also provide theoretical results in this vein, though they do so in the context of their “design-based" asymptotics that differ from those traditionally used to analyze clustered standard errors. Nonetheless, the problems with tests based on excessively coarse-clustering arise even with few clusters -- the setting of interest for our paper. We present a simple simulation to demonstrate these issues in Appendix (ref).

Given the above considerations, a researcher may choose to cluster at a fine level (e.g. individual or county) even when a coarse level of clustering (e.g. state), which is known to be valid, is also available. Nonetheless, they may be unsure if the fine level is appropriate. That is, whether observations across the fine clusters are approximately independent.

To help researchers assess the validity of their chosen clusters, we propose a modified randomization test that can be used as a robustness check for a given clustering specification. Our test requires large (fine) sub-clusters, but is justified under asymptotics that take the number of (coarse) clusters and (fine) sub-clusters as fixed. Inference is difficult in this setting because scores are not independent across sub-clusters even asymptotically, as we will explain in Section (ref). Randomization tests, which typically require some type of asymptotic independence, thus cannot be directly applied. We get around this problem by searching for worst-case values of the unobserved parameters to guard against over-rejection. We describe a simple method to search for this value, so that the computational complexity of the test is of the same order as the number of sub-clusters. This is reasonable since our test is targeted towards applications with few sub-clusters. Our test has no power against negative correlation. However, ignoring negative correlation leads to variance estimators that are too large, and is thus less of an issue if the researcher is concerned about size control when performing inference on $\beta$.

To our knowledge, there are two other tests for the level of clustering. mnw2020 proposes a test based on having large number of coarse clusters, relying on the wild bootstrap to improve finite sample performance. Meanwhile, im2016 proposes a test for the case when there are many sub-clusters. Our test, which takes the number of clusters and sub-clusters to be fixed, handles a more challenging situation, though this comes at the cost of being conservative, especially in settings with homogeneous clusters. However, as our simulations in Section (ref) show, it has competitive power given heterogeneous clusters -- a setting that could be relevant for empirical work. Indeed, our test detects correlation in the clusters chosen by gllqsx2019, demonstrating its potential usefulness in applied work (see Section (ref)). Finally, we note that the test of im2016 also has no power against negative correlation, although that of mnw2020 does not share this limitation.

aaiw17 takes a different approach to this issue. They argue for a “design-based" perspective on clustering, requiring researchers to determine ex ante the uncertainty that they face in either sampling or treatment assignment. For example, if the researcher believes that in their specific context, treatment assignment occurs at the sub-cluster level, then sub-clusters should be used for computing standard errors, regardless of whether or not residuals are correlated across the sub-clusters. While insightful, this approach requires researchers to answer an alternative question on which there is equally little theoretically guidance. We therefore develop our method under the “model-based" framework, in which the researcher has in mind some data-generating process that entails dependent clusters.

The remainder of this paper is organized as follows. Section (ref) describes our proposed test. Section (ref) presents Monte Carlo simulations. Section (ref) demonstrates an application to gllqsx2019. Section (ref) concludes. Proofs are collected in Appendix (ref).

The Proposed Test

Model and Assumptions

In the following, we assume that the researcher has conducted inference on $\beta \in \mathbf{R}$, and seeks a robustness check for the level of clustering used for said inference. As will become clear in Section (ref), using a scalar $\beta$ yields computational advantages, though the test can be feasibly computed for moderate dimensions of $\beta$. For this reason and for ease of exposition we limit our discussion to the scalar case.

Consider the linear regression:

equation[equation omitted — 130 chars of source]

where $\beta \in \mathbf{R}$ is the parameter of interest and $\gamma \in \mathbf{R}^d$ is a nuisance parameter. Suppose there are $r$ clusters, indexed by $k \in \mathcal{K}$. Within each cluster $k$, there are $q_k$ sub-clusters, indexed by $j \in \mathcal{J}_k$. Within each sub-cluster $j$, there are $n_j$ individuals indexed by $i \in \mathcal{I}_j$. Let $\mathcal{J} = \bigcup_{k \in \mathcal{K}} \mathcal{J}_k$ and $\mathcal{I} = \bigcup_{j \in \mathcal{J}} \mathcal{I}_j$. Further, let $n = \sum_{j \in \mathcal{J}} n_j$ and $q = \sum_{k \in \mathcal{K}} q_k = |\mathcal{J}|$. We also write $i \in \mathcal{I}_k$ when $i \in \mathcal{I}_j$ and $j \in \mathcal{J}_k$. In the following, we suppress dependence on $j$ and $k$ whenever this does not cause confusion.

assumptionSuppose that for every cluster $j$, there exists a vector $\Pi_j$, with a consistent estimator $\hat{\Pi}_j$, such that for all $i \in \mathcal{I}_j$: \begin{equation} X_i = W_i'\Pi_j + \varepsilon_i, \quad E[W_i\varepsilon_i] = 0 . \end{equation}

Suppose that within a sub-cluster, $W$ has full rank. Then $\hat{\Pi}_j$ can be chosen as the sub-cluster level OLS estimator of $X$ on $W$. Otherwise, we can just drop variables until we obtain a linearly independent subset $\tilde{W}$. The entries of $\hat{\Pi}_j$ corresponding to the dropped variables can then be set to $0$ while the remaining entries are chosen to be the corresponding coefficients from the sub-cluster level regression of $X$ on $\tilde{W}$. Alternatively, if the researcher is willing to assume that $\Pi_j$ is identical across clusters, $\hat{\Pi}_j$ can also be obtained from the full sample regression of $X$ on $W$. Now define

equation[equation omitted — 136 chars of source]

where $\hat{U}_i$ is the full-sample OLS residual using equation ((ref)). Suppose we know that clusters are independent, so that $E[Z_{i}Z_{i'}] = 0$ when $i \in \mathcal{I}_k, i' \in \mathcal{I}_{k'}$ and $k \neq k'$. Under this assumption, we test the null hypothesis that sub-clusters are uncorrelated:

equation[equation omitted — 146 chars of source]

against the alternative hypothesis that there exists sub-clusters within at least one cluster that exhibit correlation:

equation*[equation* omitted — 162 chars of source]

Note that changing the choice of $X_i$ and $W_i$ corresponds to testing different null hypotheses and could lead to differing outcomes. If a researcher wants to test the level of clustering used for inference on $\beta$, $X_i$ should be projected onto $W_i$. Similarly, if inference was conducted on $\gamma$, then $W_i$ should take the place of $X_i$ in equation ((ref)).

remarkA researcher interested in inference on $\beta$ only has to test the residualized hypothesis of equation ((ref)). This is because to the first order, the asymptotic distribution of \begin{equation*} \sqrt{n}\left(\hat{\beta} - \beta\right) \quad and \quad \sqrt{n_j}\left(\hat{\beta}_j - \beta\right) \end{equation*} depends only on \begin{equation*} \frac{1}{\sqrt{n}} \sum_{i \in \mathcal{I}} Z_i \quad and \quad \frac{1}{\sqrt{n_j}} \sum_{i \in \mathcal{I}_j} Z_i . \end{equation*} respectively. Hence, if the $Z_i$'s exhibit no correlation across clusters, then conducting inference using the sub-clusters is appropriate. We flesh out this argument in Appendix (ref). The fact that tests for different coefficients require different adjustments for clustering is unsurprising. A similar phenomenon arises in methods employing degrees of freedom correction for inference with a small number of clusters. Here, each slope parameter in a regression may require a test with different degrees of freedom (see ik2016 and bm2002).
remarkAs with im2016 and mnw2020, we require the researcher to specify independent clusters which nest the potentially correlated sub-clusters. While this is not always feasible, researchers seeking to test the level of clustering typically have a few choices available to them. Coarser clusters are also frequently considered to be more believably independent than the finer ones. It is natural to apply our tests in these instances.

We further assume the following:

assumptionSuppose $q$ and $r$ are fixed, but $n_j \to \infty$ for all $j \in \mathcal{J}$. Let ${Z}_{i}$ be defined as in equation ((ref)). Suppose there exists $\Omega \in \mathbf{R}^{q \times q}$ such that the $q$-vector $S_n \overset{d}{\to} S$, where \begin{equation} S_n := \begin{pmatrix} \frac{1}{\sqrt{n_1}} \sum_{i \in \mathcal{I}_1} {Z}_{i} \\ \vdots \\ \frac{1}{\sqrt{n_{q}}} \sum_{i \in \mathcal{I}_{q}} {Z}_{i} \end{pmatrix} \quad and \quad S := N \left( \mathbf{0}, \Omega \right). \end{equation} Further, let $\hat{\beta}$ and $\hat{\gamma}$ be the (joint) respective OLS estimators of $\beta$ and $\gamma$, as defined in equation ((ref)) and $\hat{\Pi}_j$ be the estimator of $\Pi_j$ as defined in equation ((ref)). Suppose: \begin{equation*} \hat{\beta} \overset{p}{\to} \beta, \quad \hat{\gamma} \overset{p}{\to} \gamma, \quad \sqrt{n_j}\left(\hat{\Pi}_j - \Pi_j\right) = O_p(1) \quad for all j \in \mathcal{J} . \end{equation*}

In other words, we assume that the errors are weakly correlated within each sub-cluster $j$. Imposing weak dependence within a (sub-)cluster is not an uncommon assumption (see for instance the discussion in crs2017 and bch2011). We note that under $H_0$, $\Omega$ is a diagonal matrix. On the other hand, under the alternative, it has a block diagonal structure due to correlation between sub-clusters.

remarkOur assumption that $S_n \to S$ does not implicitly assume that (sub-)clusters have similar sizes. Intuitively, this is because our randomization test assigns “equal weight" to each sub-cluster: each sub-cluster is normalized by its own $n_j$, and the sign of each $S_{n,j}$ contribute equally to the sign mismatch within its parent cluster. As such, heterogeneous sub-cluster sizes pose no issue for our test. Nonetheless, the quality of the asymptotic approximation is determined by $\min_{j \in [q]} n_j$, so the smallest cluster has to be large. We expand on this point in Appendix (ref) and explain how the restricted heterogeneity assumptions that are required for inference with clustered data are not needed in our case.
remarkWithout further assumptions on $Z_i$, our test requires large sub-clusters. This rules out testing the null of no clustering where there is only one observation in each sub-cluster. However, the test is valid for the null of no clustering if we are willing to assume that each $Z_i$ is symmetrically distributed around $0$. Such assumptions can be found in the econometrics literature. For example, df2008 use it to justify a wild-bootstrapped based $F$-test for the linear regression model. Nonetheless, we consider this assumption to be highly restrictive and hence justify our test via large sub-cluster asymptotics.

Test Statistic and Critical Value

In this subsection, we define the test statistic and explain the need to search over the worst case critical value. Before doing so, we first consider the infeasible test in which the true parameters -- $\beta$, $\gamma$ and $\Pi$ as defined in equations ((ref)) and ((ref)) -- are observed. Readers who are only interested in the details of implementation can skip to the end of Section (ref).

Infeasible Test

Suppose we know $\beta$, $\gamma$ and $\Pi$. Given $Y_i$ and $X_i$, we can back out $U_i$ and construct the vector $S_n^*$, whose $j^\text{th}$ entry is

equation[equation omitted — 195 chars of source]

Given $S^*_{n}$, we can then define the infeasible test statistic:

equation[equation omitted — 223 chars of source]

The inner sum is the net number of positive ${S}^*_{n,j}$ within each cluster $k$. Intuitively, if the sub-clusters are independent, the net number of positive $S^*_{n,j}$ should be close to $0$. Conversely, if they are positively correlated, this number will be large in absolute value, since many sub-clusters will have $S_{n,j}^*$ of the same sign. On the other hand, if they are negatively correlated, this number will be more concentrated around $0$ than in the independent case. As will become clear below, our test interprets large absolute values of $T(S_n^*)$ as violation of the null. For this reason, we it will not have power against negative correlation.

remarkThere are two advantages to having a test statistic that depends only on the sign of the $S_{n,j}^*$'s. Firstly, large and small realizations of $S_{n,j}^*$ contribute the same amount to $T(S_n^*)$. As such, the performance of our test is not affected even if sub-clusters have wildly differing variances, a source of heterogeneity that may be important in applied work. We demonstrate this robustness property via simulations in Section (ref). Secondly, the feasible version of this test requires searching over the worst case values of the test statistic. As will become clear in Section (ref), this search is simplified by our choice of test statistic.

Next, denote by $\mathbf{G}$ the set of sign changes. $\mathbf{G}$ can be identified with the set of $g \in \{-1, 1\}^{q}$ so that

equation*[equation* omitted — 119 chars of source]

Now let $p^*(S_n^*)$ be the proportion of $T\left(g{S}^*_n\right)$ that are no smaller than $T\left({S}^*_n\right)$:

equation[equation omitted — 186 chars of source]

The test rejects the null hypothesis when $p(S_n^*)$ is small -- that is, when $T\left({S}^*_n\right)$ is extreme relative to $T\left(g{S}^*_n\right)$:

align[align omitted — 136 chars of source]

The intuition for the randomization test is as follows. Since $S^*_{n,j}$ involves only units within the same sub-cluster, under the null hypothesis, $S_n^*$ converges to a mean-zero normal distribution with independent components. Independence, together with symmetry of normal random variables about their means, implies that for any $g \in \mathbf{G}$, $gS^*_n$, has the same distribution as $S_n^*$. Hence, the randomization distribution $\{T(gS_n^*)\}_{g\in \mathbf{G}}$ is in fact the distribution of $T(S_n^*)$ conditional on the values of $|S_n^*|$, where $|\cdot|$ is applied component-wise. Rejecting the null hypothesis when we observe values of $T(S_n^*)$ that are extreme relative to $\{T(gS_n^*)\}_{g\in \mathbf{G}}$ therefore leads to a test with the correct size.

Note that the randomization test defined above is non-randomized. Randomization tests can also employ a randomized rejection rule for the situation when

equation*[equation* omitted — 329 chars of source]

Using a randomized rejection rule, we have that provided the necessary symmetry properties hold in finite sample, the randomization test will have size equal to $\alpha$ exactly. The test defined in equation ((ref)) is conservative since it never rejects when the above situation occurs. However, we present the deterministic version since the test that we propose is based on it.

Na\"ive Test

Tests based on $S_n^*$ are infeasible since $\beta$, $\gamma$ and the $\Pi_j$'s are unknown. Suppose we simply replaced $Z_i$ with $\hat{Z}_i$ and performed the randomization test with the estimated scores. It turns out that this procedure is incorrect. To see this, let $\tilde{S}_n$ be $S_n^*$ but with $\hat{Z}_i$ replacing $Z_i$. Then we can write each component of $\tilde{S}_n$ as:

align[align omitted — 525 chars of source]

In the above equation, $S_{n,j}^*$ is the part that is informative about cluster structure. However, each component now has an additional nuisance term $A_j$ that does not go away under asymptotics that take the number of sub-clusters to be fixed. Because $\hat{\beta} - \beta$ is common across the $A_j$'s, it induces correlation across $\tilde{S}_{n,j}$ even when the $S_{n,j}^*$'s are independent, leading potentially to over-rejection. Addressing this complication which does not arise in frameworks taking $q \to \infty$ results in the conservativeness of our test.

Feasible Test

If we knew $\hat{\beta} - \beta$, we could back out ${S}^*_{n,j}$ for the randomization test using equation ((ref)). Since that is not possible, we propose to search over values of $\hat{\beta} - \beta$ to ensure that the test controls size when the unobserved term takes on extreme values.

For a given $\lambda \in \mathbf{R}$, let $\hat{S}_n({\lambda})$ be $q\times 1$ vector whose $j^\text{th}$ entry is the following term:

equation*[equation* omitted — 267 chars of source]

Note that $\hat{S}_{n,j}(\hat{\beta} - \beta) = S^*_{n,j} + o_p(1)$. Define:

equation*[equation* omitted — 228 chars of source]

For a given $\lambda$, this is just the test statistic in equation ((ref)) but with $\hat{S}_{n,j}(\lambda)$ taking the place of $S_{n,j}^*$. As before, we denote by $\mathbf{G}$ the set of sign changes and write:

equation*[equation* omitted — 152 chars of source]

Now let $p(\hat{S}_n(\lambda))$ be the proportion of $T\left(g\hat{S}_n(\lambda)\right)$ that takes on extreme values relative to $T\left(\hat{S}_n(\lambda)\right)$:

equation[equation omitted — 220 chars of source]

We can then define the randomization test as:

align[align omitted — 175 chars of source]

We can then prove the following result:

theoremUnder assumptions (ref) and (ref), $\underset{n \to \infty}{\lim\sup} \,\,\mathbb{E}[\phi_n] \leq \alpha$.

The test is a two-stage process. In the first stage, it searches for the value of $\lambda$ that leads to the largest $p$-value. In the second stage, the test rejects if this worst-case $p$-value is still smaller than the desired level of significance $\alpha$. Since the worst-case $p$-value bounds the true $p$-value from above, the rejection rule based on the worst-case $p$-value must be conservative.

As the Monte Carlo simulations in Section (ref) shows, the test has size that could be much smaller than $\alpha$ under the null hypothesis. However, the same simulations also show that the test has reasonable power under the alternative hypothesis, particularly in settings where clusters are heterogeneous in their variances. The potential usefulness of our test is further seen in the empirical application (Section (ref)), where it detects dependence in the clusters chosen by gllqsx2019.

remarkThe worst-case test has no power if $r = 1$ since $\lambda = \text{median}(\{\tilde{S}_{n,j}\})$ will set exactly half the signs of $\hat{S}_{n,j}(\lambda)$ to be positive and half to be negative, so that the signs are completely balanced. However, this is no longer true with $r > 1$ since only a single value can be chosen to balance signs across multiple clusters. The implementation procedure provides further intuition for power in this test. See the next subsection.
remarkAs with standard randomization tests, $|\mathbf{G}|$ may sometimes be too large so that computation of $p(\hat{S}_n(\lambda))$ becomes onerous. In these instances, it is possible to replace $p(\hat{S}_n(\lambda))$ with a stochastic approximation. Formally, let $ \hat{\mathbf{G}} = \left\{ g^1, ... ,g^B \right\}~, $ where $g^1$ is the identity transformation and $g^2, ..., g^B$ are i.i.d. Uniform($\mathbf{G}$). Using $\hat{\mathbf{G}}$ instead of ${\mathbf{G}}$ in equation ((ref)) does not affect validity of theorem (ref). For implementation, we follow crs2017 in evaluating the $p(\hat{S}_n(\lambda))$ completely when $q \leq 10$ and approximating it with $B = 1000$ when $q > 10$.
remarkWe advocate the use of our test as a robustness check, after a researcher has chosen a level of clustering for inference, in the same spirit that manipulation tests are routinely used in studies with regression discontinuity designs, or in tests for pre-trends in studies involving difference-in-differences. In particular, the original inference results should be presented with results of the current test, regardless of the outcome. Conceptually, this is different from using the test as a pre-test to select the level of clustering prior to inference. The distinction is important as pre-testing is known to induce uniformity issues, where inference in the second stage (on $\beta$) suffers from distortion due to mistakes in the pre-test (that happen with positive probability). These same concerns are articulated by im2016, who argue that their test “merely provides empirical evidence on the plausibility of one particular clustering assumption”. We take exactly the same view of our test.

Implementation

In this subsection we describe an efficient way of searching for $\lambda \in \mathbf{R}$. This search is simplified by the fact that $p(\hat{S}_n(\lambda))$ depends only on the sign of $\hat{S}_{n,j}(\lambda)$'s. As such, to find $\sup_{\lambda \in \mathbf{R}}$, we only need to search over sign combinations of $\hat{S}_{n,j}$. When $\beta$ is scalar, the search can be completed in $O(q)$ time. This is reasonable since the test is designed for use when $q$ is small.

Suppose for now that $\sum_{i \in \mathcal{I}_j} \left(X_i - W_i'\hat{\Pi}_j\right)^2 > 0$ for all $j \in \mathcal{J}$. Define:

align*[align* omitted — 170 chars of source]

Then, $ \hat{S}_{n,j}({\lambda}) \geq 0 \Leftrightarrow R_j + \lambda \geq 0~. $ Sort the values of $R_j$'s so that $ R^{(1)} \geq R^{(2)} \geq ... \geq R^{(q)}~. $ We must have that $R^{(j)} +\lambda \geq 0 \Rightarrow R^{(j')} +\lambda \geq 0$ for all $j' \leq j$. Let $\hat{S}_{n, (1)}(\lambda), ..., \hat{S}_{n, (q)}(\lambda)$ denote the values of $\hat{S}_{n,j}(\lambda)$ corresponding to $R^{(1)}, ..., R^{(q)}$. Therefore, we only need to consider sequences of the form

equation*[equation* omitted — 123 chars of source]

for some cut-off $j$. Since the $p$-value, as defined in equation ((ref)), depends only on the sign of $\hat{S}_n$, we can compute it using $\check{S}_n$ in the place of $\hat{S}_{n,j}(\lambda)$:

equation*[equation* omitted — 128 chars of source]

Here, we see that even when we are searching over the worst case $\lambda$, we are only allowed to choose the cut-off point at which the signs change. We can therefore complete the search with no more than $q$ randomization tests. Assuming that the time it takes for each test is $O(1)$, the procedure takes $O(q)$ time. The restriction that $\check{S}_{n, (j)} \geq \check{S}_{n, (j')}$ for all $j \leq j'$ also gives the test power. If all combinations of signs for the $S_{n,j}$'s were allowed, the test will always return a $p$-value of 1 and will have no power.

Finally, suppose there are sub-clusters such that $\sum_{i \in \mathcal{I}_j} \left(X_i - W_i'\hat{\Pi}_j\right)^2 = 0$. We can repeat the above procedure excluding these sub-clusters. In the final step, we set $\check{S}_{n,j}$ corresponding to these clusters to 0. Hence,

remarkWe can further reduce computation time by the following. Let $$R^+_k =\min_{j' \in \mathcal{J}_k} \left\{R_{j'} \text{ greater than or equal to $>0.5$ of } \{R_j, j \in \mathcal{J}_k\} \right\}$$ be the “upward-conservative" median. Also define the “downward-conservative" median: $$R^-_k = \max_{j' \in \mathcal{J}_k} \left\{R_{j'} \text{ less than or equal to $>0.5$ of } \{R_j, j \in \mathcal{J}_k\} \right\}~.$$ Now let $R^+ = \max_{k \in \mathcal{K}} R_k^+ \quad \mbox{and} \quad R^- = \min_{k \in \mathcal{K}} R_k^-~.$ We only need to consider cutoffs below $R^+$. Setting the sign cutoff at the argmax of $R^+$ results in situation in which all clusters have at least half of their entries being $-1$. If we now set the extreme $S_{n,j}$'s to $-1$, this will increase the net number of $-1$'s in all clusters. Since our test is based on sign imbalance within clusters, such sequences will lead to a strictly larger test statistic and smaller $p$-values than if they were set to 1. For the same reason, we only need to consider cutoffs above $R^-$.

We summarise the implementation procedure in Algorithm (ref).

algorithm[algorithm omitted — 1,262 chars of source]

Comparison with Existing Tests

To our knowledge, two other tests have been proposed for the level of clustering. They take either the number of sub-clusters in each cluster to infinity or the number of clusters to infinity. We assume both to be fixed. For ease of exposition, we restrict our discussion of these tests to the univariate case.

im2016 (IM hereafter) adopts an asymptotic framework that takes $q_k \to \infty$ for all $k \in \mathcal{K}$. Consider estimating a regression coefficient cluster-by-cluster. Let $\hat{\beta}_k$ denote coefficients estimated using only cluster $k$. The IM test is based on the asymptotic distribution of an estimator for the variance of $\frac{1}{r} \sum_{k = 1}^r\hat{\beta}_k$. Let this variance be denoted by $V$ and let $\hat{\Omega}^\text{CCE}_k$ be the cluster-robust variance estimator for $\hat{\beta}_k$, where the clustering is done at the sub-cluster level using $j \in \mathcal{J}_k$. Under the null hypothesis, $\hat{\Omega}^\text{CCE}_k$ consistently estimates the variance of each $\hat{\beta}_k$.

Under either the null or the alternative, but maintaining the assumption that coarse clusters are independent, consider estimating $V$ by:

equation*[equation* omitted — 152 chars of source]

IM show that under the null, $\hat{V} \overset{d}{\to} V^W$, where $V^W = \frac{1}{r-1} \sum_{k = 1}^r ({W}_k - \bar{W})^2$ and

equation*[equation* omitted — 118 chars of source]

The IM test constructs a reference distribution $\hat{V}^W$ by drawing $W$ from $$N(0, \text{diag}(\hat{\Omega}^\text{CCE}_1, \hat{\Omega}^\text{CCE}_2, ..., \hat{\Omega}^\text{CCE}_r))$$ and seeing if $\hat{V}$ is larger than the $\left( 1-{\alpha}\right)^\text{th}$ quantile of $\hat{V}^W$.

There are two limitations to the IM test that our test does not share. Firstly, they require the regression to be estimated cluster-by-cluster. This would be infeasible in, for example, differences-in-differences set ups where treatment varies at the cluster level. Secondly, since their asymptotics take $q_k \to \infty$, we expect the test to have poor properties when $q_k$ is small. Instead, our test is expected to have good properties even when $q_k$ is small as long as $n_j$ is large. These benefits come at a cost. We expect our test to perform worse if observations within sub-clusters are highly correlated, whereas the IM test allows unrestricted covariance within sub-clusters. Our test is also conservative under the null hypothesis. We note also that neither test has power against negative correlations. This is because both tests use test statistics that take on large value relative to their reference distributions only when there is positive correlation.

mnw2020 (MNW hereafter) considers an asymptotic framework that takes $r \to \infty$. In the same spirit as IM, the MNW test is a Hausman-type test based on the variance of regression coefficients. Consider the full sample regression coefficient $\hat{\beta}$. Under the null hypothesis, the (full-sample) cluster-robust covariance estimator at the sub-cluster level, denoted, $\hat{\Omega}^\text{CCE}_J$, is consistent for the asymptotic variance-covariance matrix.

Under either the null or the alternative, but maintaining the assumption that coarse clusters are independent, the (full-sample) cluster-robust covariance estimator at the cluster level, denoted, $\hat{\Omega}^\text{CCE}_K$, is consistent for the asymptotic variance-covariance matrix. Under the null hypothesis, the authors show that their test statistic converges to a standard normal distribution: $ \frac{\hat{\Omega}^\text{CCE}_K - \hat{\Omega}^\text{CCE}_J}{\hat{V}^{MNW}} \overset{d}{\to} N(0,1) $ for an appropriately defined $\hat{V}^{MNW}$.

It is well known that the cluster-robust covariance estimator can be severely biased when $r$ is small. In order to deal with such situations, the authors propose to conduct the test using wild (sub-)cluster bootstrap. They prove the consistency of this approach in their large-$r$ framework, showing power even against alternatives with negative correlations.

Compared to the MNW test, our test is theoretically justified when both $r$ and $q$ are small, provided that $n_j$'s are large. Our test could therefore be preferable in such applications since it is presently not known if the MNW test remains valid once we take $r$ and $q$ to be fixed. However, as with the IM test, the MNW test allows unrestricted covariance within sub-clusters, whereas our test is expected to have poor performance if observations within sub-clusters are highly correlated. Our test is also conservative relative to the MNW test. On the other hand, simulation evidence suggests that it has comparable performance with the MNW test when clusters have differing variances (see Section (ref)).

Monte Carlo Simulations

In this section, we examine the finite sample performance of our worst-case randomization test (WCR) together with the IM and bootstrap version of the MNW tests via Monte Carlo simulations. We also study the performance of the na\"ive randomization test (NR) as described in Section (ref). We consider two data generating processes described below.

Model 1: Model 1 is defined by the following:

gather*[gather* omitted — 301 chars of source]

In particular, we set $X_{t,j,k} = \beta = 1$ and $\phi = 0.25$. Errors are correlated within a sub-cluster, according to an $AR(1)$ process, with autocorrelation coefficient $\phi$. $\rho$ captures the importance of cluster level shock. Since $\frac{1}{\sqrt{1 - \phi^2}} U_{t,j,k}$ has unit variance, $\rho$ is exactly the relative variance of cluster- to sub-cluster-level shocks. $\sigma_{j,k}$ controls the variance of the unobserved term in each cluster $k$. Here in Section (ref), we set $\sigma_{j,k} = 1$ for all $j \in \mathcal{J}, k \in \mathcal{K}$. In Section (ref), we explore the consequences of cluster heterogeneity by varying $\sigma_{j,k}$.

Model 2: This is the model used in the simulations of mnw2020, with the constant omitted. Let $m_k = \sum_{j \in \mathcal{J}_k} n_j$ be the total number of observations in cluster $k$. Let $U_k$ be the $m_k \times 1$ vector of $U_{t,j,k}$ for all observations in cluster $k$. Then

align*[align* omitted — 116 chars of source]

where $\xi_k$ is a $10 \times 1$ vector distributed as:

align*[align* omitted — 167 chars of source]

and $W_\xi$ is the $m_k \times 10$ loading matrix with the $(i,j)^\text{th}$ entry $\mathbf{1}\left\{ j = \lfloor (i-1)10/m_k \rfloor + 1 \right\}$. Under this model, $\frac{1}{10}$ of the observations in each cluster are correlated because they depend directly on the same $\xi_{k, l}$. In addition, there is correlation between the $\xi_{k,l}$'s since it is generated according to an AR(1) process. Observations are then ordered so that every sub-cluster contains the same number of observations that depend on each $\xi_{k,l}$. Finally, $\beta = (1,1)'$ and the two covariates are independent and generated in the same way as $U$. This model features more complex correlations between and within the sub-clusters. Clusters are independent and identically distributed. As in Section 5.2 of mnw2020, we set $\phi = 0.5$. $\rho$ here is directly comparable to $w_\xi$ in their simulations.

For our simulations, we perform the test at the 5% level. 1,000 Monte Carlo simulations were drawn for each combination of the parameters. The non-standard reference distribution in IM is evaluated using 1,000 Monte Carlo draws. Wild bootstrap in MNW is evaluated using $399$ draws as in their simulations.

Performance over values of $r$, $q_k$ and $n_j$

To understand the size and power of each of our tests in scenarios with few clusters and few sub-clusters, we consider equal-sized clusters and sub-clusters, with $r \in \{4, 8, 12\}$, $q_k \in \{4, 8, 12\}$ and $n_j \in \{25, 50, 100\}$. We consider $\rho \in \{0, 0.5\}$.

Table (ref) presents results under the null hypothesis ($\rho = 0$). Across the two models, we see that regardless of $r$, the IM test performs poorly when $q_k$ is small. With $q_k = 4$, type I error is between 15% and 20%. By $q_k = 12$, however, the size is between 6-7%. Comparatively, our test, which is highly conservative, has type I error less than 2% across all values of $q_k$. The MNW and NR tests perform well across the board. Table (ref) presents results under the alternative $\rho = 0.5$. Relative to the IM and MNW tests, our test has power that is consistently lower. In particular, our test does poorly when $q_k$ is small. This is the weakness of the worst-case approach.

table[table omitted — 4,267 chars of source]
table[table omitted — 4,273 chars of source]

Figure (ref) presents power of the tests for $r = 8$, $q_k = 8$, $n_j = 100$ as we vary $\rho$ from $0$ to $2$ in model 1 and $0$ to $1$ in model 2. Across the two models, we see that the IM and MNW tests have greater power than our test. However, as $\rho$ increases, our test quickly catches up in power.

figure[figure omitted — 265 chars of source]

Effect of Cluster-Level Heterogeneity

The previous section suggests that our test has poor performance compared to all other tests, including NR. However, a different picture emerges once we allow clusters and sub-clusters to be heterogeneous in their variances.

We first consider what happens when clusters are heterogeneous. Specifically, we return to model 1 but with $\sigma_{j,1} \in \{5,10,15\}$. That is, when all sub-clusters in cluster 1 are much noisier than the rest. Figure (ref) plots power curves with $r = 8, q_k = 12, n_j = 100$ for $\sigma_{j,1} \in \{5,10,15\}$. These curves are directly comparable with Figure (ref). Starting from the within test comparison, we see that the performance of our test is unaffected by $\sigma_{j,1}$. However, power of IM and MNW quickly degrade as $\sigma_{j,1}$ increases. Turning to the across test comparison, we see that the tests perform similarly when $\sigma_{j,1} = 5$. As $\sigma_{j,1}$ increases to 10, our test starts to have more power than the IM and MNW tests for $\rho \geq 1$. The across test comparison also shows how the NR test fails to control size. In particular, when $\sigma_{j,1}$, an NR test with nominal size 5% could wrongly reject over 40% of the time.

We see the same patterns when sub-clusters are heterogeneous. Consider again model 1 but with $\sigma_{1,k} \in \{5,10,15\}$. That is, when the first sub-cluster in each cluster is much noisier than the rest. Figure (ref) presents the results. Again, our test is not affected by changing $\sigma_{1,k}$. The power of the IM test falls by a large extent as $\sigma_{1,k}$ increases. The MNW test is also negatively affected by $\sigma_{1,k}$, though less so than the IM test.

figure[figure omitted — 391 chars of source]
figure[figure omitted — 325 chars of source]

All in all, the simulation evidence suggests that our test manages to maintain type I error below $\alpha$ when $q$ is small, whereas the IM and NR tests may see size distortion in such a setting. The cost of size control in a fixed $q$ setting is that the procedure is very conservative. This conservativeness limits the power of our test. However, the performance of our test is less sensitive to heterogeneous variances within and across clusters, such that it could be more powerful than the IM and MNW tests when some clusters or sub-clusters are much noisier than others. Hence, our test is suited for applications with small $q$ and heterogenous clusters. Indeed, as we will see in the next section, our test detects dependence in the clusters of gllqsx2019, demonstrating its potential relevance for empirical work.

Application: gllqsx2019

In recent years, the poor performance of American students in assessment tests such as the Programme for International Student Assessment (PISA) has raised concerns among policymakers. gllqsx2019 argues that the testing gap reflects, among other things, the low effort that American students put in on tests, especially when compared to their higher scoring counterparts in other countries.

The authors test their hypothesis by a randomized controlled experiment in which students were rewarded with cash for correct answers in a 25-question test. Those assigned to the treatment group were offered roughly \$1 USD per correct answer, while the control group received no payment. Students were informed right before the test started to prevent them from changing their effort in test preparation. The experiments were conducted at 4 schools in Shanghai and 2 schools in the US. Due to logistical reasons, the authors randomized treatment at the class level for some schools and individual level for others.

Various regression analyses were conducted to study the effect of treatment on test-taking effort and test performance. Panel A in Table 3 examines whether monetary incentive increased the probability that students attempt a given question -- a proxy for effort. It does so by estimating the following equation:

equation*[equation* omitted — 67 chars of source]

Here, the unit of analysis is a question and $Y_{qi}$ is an indicator for whether student $i$ attempted question $q$. $Z_i$ is the treatment indicator and $W_i$ is a vector of control variables, which include terms such as gender, ethnicity as well as question number fixed effects. We focus on Column 1 in Panel A, which looks at US students' responses to all 25 questions in the test, and Column 4, which looks at Shanghai students' responses to the same test.

The authors present their linear regression estimate of $\beta$, together with standard errors clustered at the level of randomization. However, other levels of clustering are plausible:

itemize• G: Group Level, that is, the level of randomization. • S: School Level. • SY: Experiments in Shanghai schools were conducted in 2016 and then 2018. We could plausibly interact school and year of experiment. • ST: Schools in the US separate students into tracks (Honors, Regular, Others). We could plausibly interact school and track.

We will refer to these levels of clustering by their initials hereafter. More information on the sizes of clusters can be found in appendix (ref).

While the authors chose to cluster their standard errors by $G$, it seems reasonable to be concerned about correlation across individuals within the same school or among those who took the test in the same year. If these clusters were not independent, $t$-tests using the presented standard errors could lead to the wrong conclusions.

table[table omitted — 1,117 chars of source]

Table (ref) presents the OLS estimates from gllqsx2019 as well as the $p$-values that would be obtained from testing the null hypothesis that $\beta = 0$ using several methods. Specifically, we consider the wild cluster bootstrap (cgm2008), approximate randomization tests (crs2017) and the t-distribution based procedure of ibragimov2010t, denoted IM2010. We perform these tests using the various plausible levels of clustering. For Column 1, we consider the increasingly coarse levels of clustering $G$, $ST$ and $S$. For the US, there are no schools sampled over multiple years, so $SY$ is the same as $S$. For Column 4, we consider the increasingly coarse levels of clustering $G$, $SY$ and $S$. In Shanghai schools, students are not separated by track, so $ST$ is the same as $S$.

remarkgllqsx2019 present clustered standard errors but do not use them for inference. Instead, they conduct randomization inference by permuting treatment status as in young2019channeling. This procedure tests the null hypothesis that the distribution of the $Y_{qi}$'s are the same with and without treatment. This is a stronger null hypothesis than the null of $0$ average treatment effect ($\beta = 0$). We believe that the latter hypothesis is typically the one of interest and test it in our Table (ref).

Turning to the results, for column 1, we see that CCE SE's decrease as we move to increasingly coarse levels of clustering. Correspondingly, $p$-values from CCE-based $t$-tests decrease as we coarsen the clusters. Such a pattern is typically interpreted as arising from the downward bias of CCEs with few clusters (mhe2008), so that these $p$-values would be considered unreliable. Faced with downward bias, practitioners commonly turn to the wild cluster bootstrap. With this method, the $p$-values increase as we coarsen the clusters. While clustering at $G$ and $ST$ may lead one to conclude that there is strong evidence that $\beta \neq 0$, the $p$-value at $S$ suggests the absence of strong evidence. The same phenomenon arises with approximate randomization tests: at $ST$ there appears to be strong evidence that $\beta \neq 0$. At $S$, this is no longer true. With IM2010, the test does not reject in either case. We note that ART and IM2010 cannot be applied with $G$ as the chosen level of clustering, since both methods require $\beta$ to be estimated cluster-by-cluster. The results for column 4 are qualitatively similar. At $G$, CCE-based $t$-test and the wild cluster bootstrap find strong evidence that $\beta \neq 0$. This conclusion is overturned once we cluster at either $SY$ or $S$.

table[table omitted — 817 chars of source]

To assess the validity of the above specifications, we apply our WCR test, the IM test and the MNW tests. Table (ref) presents the resulting $p$-values. The notation $G \to S$ means that the null hypothesis involves sub-clusters $G$ and coarse clusters $S$. For Column 1, clustering at $G$ appears to be appropriate, as all 3 tests fail to reject the null hypotheses $G \to ST$ and $G \to S$. For Column 4, all 3 tests find strong evidence that sub-clusters $G$ are inappropriate. The WCR test has higher $p$-values than the IM and MNW tests, likely due to its lower power. Nonetheless, they are close to 5%. The WCR and MNW tests do not reject the null hypothesis for $SY \to S$, whereas the IM test does. Given the that there are at most 2 school$\times$year per school, the IM test is likely to over-reject. As such, we consider the conclusion of the WCR and MNW test to be more reliable in this instance. Thus, results based on clustering at $SY$ are plausible.

All in all, we see that settings with varying numbers of clusters and sub-clusters arise in empirical work. Our test, designed for applications with few clusters and sub-clusters is relevant and appears to work well in practical settings.

Conclusion

We propose to test for the level of clustering in a regression by means of a modified randomization test. We show that the test controls size even when the number of clusters and sub-clusters are small, provided that the size of sub-clusters are relatively large. This is a challenging situation not accommodated by existing tests. To ensure size control, our procedure may be conservative when clusters are homogeneous. However, in settings with heterogeneous clusters, it has power that is comparable with other tests. As such, our test can be useful when the researcher faces an application with few sub-clusters, particularly when these clusters are likely to be heterogeneous. Finally, we note that the test is easy to implement and could serve as a helpful robustness check to researchers working with clustered data. An R package is available from the author's website.

center[center omitted — 48 chars of source]
description• Technical details including proof of theorem 1, details concerning the application and additional Monte Carlo simulations.