EconBase
← Back to paper

On the implementation of Approximate Randomization Tests in Linear Models with a Small Number of Clusters

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.

83,635 characters · 14 sections · 50 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.

On the Implementation of Approximate Randomization Tests in Linear Models with a Small Number of Clusters

spacing{1.1} \begin{abstract} This paper provides a user's guide to the general theory of approximate randomization tests developed in canay/etal:17 when specialized to linear regressions with clustered data. An important feature of the methodology is that it applies to settings in which the number of clusters is small -- even as small as five. We provide a step-by-step algorithmic description of how to implement the test and construct confidence intervals for the parameter of interest. In doing so, we additionally present three novel results concerning the methodology: we show that the method admits an equivalent implementation based on weighted scores; we show the test and confidence intervals are invariant to whether the test statistic is studentized or not; and we prove convexity of the confidence intervals for scalar parameters. We also articulate the main requirements underlying the test, emphasizing in particular common pitfalls that researchers may encounter. Finally, we illustrate the use of the methodology with two applications that further illuminate these points: one to a linear regression with clustered data based on mengQY2015 and a second to a linear regression with temporally dependent data based on munyo2015first. The companion {\tt R} and {\tt Stata} packages facilitate the implementation of the methodology and the replication of the empirical exercises. \end{abstract}

KEYWORDS: Randomization tests, linear regression, clustered data, time series

JEL classification codes: C12, C14

\thispagestyle{empty} \setcounter{page}{1}

Introduction

This paper provides a user's guide to the general theory of approximate randomization tests (ARTs) developed in canay/etal:17 when specialized to linear regressions with clustered data. Here, clustered data refers to data that may be grouped so that there may be dependence within each cluster, but distinct clusters are approximately independent in a way to be made precise below. Such data is remarkably common, including not only data that are naturally grouped into clusters, such as villages or repeated observations over time on individual units, but also data with weak temporal dependence, in which pseudo-clusters may be formed using blocks of consecutive observations. An important feature of the methodology is that it applies to commonly encountered settings in which the number of clusters is small -- even as small as five. In this respect, the proposed methodology contrasts sharply and meaningfully with many commonly employed methods for inference in such settings. We briefly elaborate on this point in our discussion of related literature below.

A principal goal of this paper is to make the general theory developed in canay/etal:17 more accessible by providing a step-by-step algorithmic description of how to implement the test and construct confidence intervals for the quantity of interest in these types of settings. In order to do so, we develop three novel results concerning the methodology in Section (ref). Our first result shows that what we view as the most natural implementation of the test, as described in Algorithm (ref), is numerically equivalent to an alternative implementation based on weighted scores (see Algorithm (ref)). Our second result shows that when the parameter of interest is a scalar parameter, studentizing or not the $t$-statistic entering the test does not affect the results of the test or the associated confidence intervals. We therefore focus on the unstudentized statistic in Algorithm (ref). Finally, our third result shows that the confidence sets for scalar parameters that are conceptually described by test inversion are indeed a closed interval of the real line. This further leads to a simple closed form expression for the lower and upper bound of the confidence intervals (see Algorithm (ref)). These results are new to this paper and play an important role in developing simple algorithms for the implementation of ARTs.

We additionally provide a discussion of the main requirements underlying the test in Section (ref). These requirements essentially demand that the quantity of interest is suitably estimable cluster-by-cluster. As discussed further in Section (ref), when this is not satisfied, a researcher need not conclude that it is not possible to exploit the results in canay/etal:17. Instead, several remedies are possible, including clustering more coarsely or changing the specification to ensure that this requirement is satisfied. We provide two applications that further elucidate these points: one to a linear regression with clustered data based on mengQY2015 and a second to a linear regression with temporally dependent data based on munyo2015first. The required software to replicate these empirical exercises and to aid researchers wishing to employ the methods elsewhere is provided in both {\tt R} and {\tt Stata}.\footnote{The {\tt Stata} and {\tt R} packages {\tt ARTs} can be downloaded from \hyperlink{http://sites.northwestern. edu/iac879/software/}{http://sites.northwestern. edu/iac879/software/}.}

The methodology described in this paper is part of a large and active literature on inference with clustered data. Following bertrand2004much, researchers are acutely aware of the need to adjust inferences appropriately to account for this sort of dependence. Many of the most commonly employed methods for doing so, however, are inadequate for the unusually common situation in which the number of clusters is small. Conventional wisdom suggests that the number of clusters is small when it is less than forty. For example, the method described in liang1986longitudinal, which has enjoyed considerable popularity due to its availability in software packages such as {\tt Stata}, is widely acknowledged to perform poorly when this rule-of-thumb is not satisfied. Similarly, the cluster wild bootstrap described in cameron2008bootstrap requires either a sufficiently large number of clusters or, as shown by canay/santos/shaikh:20, stringent homogeneity across clusters, to perform reliably. As explained further in Section (ref), the methods developed in canay/etal:17 and described in this paper, require neither a large number of clusters nor such homogeneity across clusters. We note that the methods by ibragimov2010t,ibragimov/muller:16, which are closely related to the ones described here, also do not require such restrictions, but are generally less powerful and permit testing a less rich variety of hypotheses. See canay/etal:17 for further discussion of these points as well as conley2018inference for an insightful and thorough review of the related literature more broadly.

The remainder of this paper is organized as follows. In Section (ref), we first formalize the setting and establish some notation. We then describe the implementation of approximate randomization tests (ARTs) in an algorithmic fashion, including how to use these tests to construct confidence intervals for the quantity of interest. In Section (ref) we present three results that play an important role in developing these algorithms. In Section (ref), we articulate the main requirements underlying the tests and discuss remedies for cases where these requirements are not satisfied. Our two empirical applications are contained in Section (ref). Finally, we provide some concluding remarks in Section (ref).

Review of ARTs in regression models

We start by reviewing the inference approach proposed by canay/etal:17 in the context of a linear regression model with clustered data. In order to do so, we index clusters by $j \in J \equiv \{1, \ldots, q\}$ and units in the $j$th cluster by $i \in I_{n,j} \equiv \{1, \ldots, n_j\}$. We also denote by $n=\sum_{j=1}^q n_j$ the total number of observations. The observed data consists of an outcome of interest, $Y_{i,j}$, and a vector of covariates, $Z_{i,j} \in \mathbf R^{d_z}$, that are related through the equation

equation[equation omitted — 77 chars of source]

where $\beta \in \mathbf{R}^{d_{z}}$ are unknown parameters and our requirements on $\epsilon_{i,j}$ are explained below in Section (ref). Our goal is to test

equation[equation omitted — 123 chars of source]

for given values of $c\in \mathbf R^{d_z}$ and $\lambda \in \mathbf R$, at level $\alpha \in (0,1)$. An important special case of this framework is a test of the null hypothesis that a particular component of $\beta$ equals a given value, i.e.,

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

for some $\ell \in \{1,\dots,d_{z}\}$, by simply setting $c$ to be a standard unit vector with a one in the $\ell$th component and zeros otherwise. More generally, the approach we describe below extends immediately to the case where the hypothesis of interest involves multiple elements of $\beta$, in which case the test becomes

equation[equation omitted — 107 chars of source]

for a given $p\times d_z$-dimensional matrix $R$ and $p$-dimensional vector $\Lambda$, at level $\alpha \in (0,1)$.

ARTs were developed more generally in canay/etal:17 and admit a variety of different applications that go beyond the linear model considered here. For example, the method accommodates non-linear models, non-linear hypotheses, or even applications that go beyond inference with a small number of clusters (e.g., canay/kamat:18 develop a variation that applies to inference in the regression discontinuity design). Here, we abstract away from the generality of the method and focus on the steps needed to use ARTs to test the null hypothesis in (ref) in the context of the model in (ref).

How to implement ARTs

The most straightforward way to test the hypotheses in (ref) via ARTs is by following the steps described in Algorithm (ref) below.

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

Algorithm (ref) involves five steps that are easy to implement from a computational standpoint, but some of the steps deserve some clarification. Step 1 involves $q$ within-cluster regressions that lead to $q$ estimates of $\beta$. This essentially demands that the parameter $\beta$ is identified cluster-by-cluster, and may fail to hold if some of the variables in the vector $Z_{i,j}$ are constant within cluster. We discuss possible remedies for this problem in Section (ref) and illustrate their use in one of the applications in Section (ref). An important feature of the method is that from Step 2 onwards, the original data is no longer needed as all the calculations only involve the $q$ estimators of the parameter $\beta$ obtained in Step 1.

Step 2 defines a type of unstudentized $t$-statistic that is appropriate for the null hypothesis in (ref). We discuss the connection to its studentized version in Section (ref) below. If the null hypothesis of interest is the one in (ref), then a Wald-type test statistic could be used instead, i.e.,

equation[equation omitted — 163 chars of source]

where

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

Step 3 does not require one to recompute the estimates of $\beta$. It rather uses the $q$ estimates from Step 1 and applies sign changes to the $q$-dimensional vector $\{S_{n,j}:j \in J\}$. Since the cardinality of $\mathbf G$ is $|\mathbf G|=2^q$, it exceeds $2000$ when $q>10$ and in such cases it may be convenient to use a stochastic approximation. This may be done while still controlling the rejection probability under the null hypothesis canay/etal:17. Formally, in this case we let

equation[equation omitted — 79 chars of source]

where $g^1=\iota\equiv (1,\dots,1)$ is the identity vector and $g^{b}=(g^{b}_1,\dots,g^{b}_q)$, for $b=2,\dots,B$, are i.i.d.\ Rademacher random variables; i.e., each $g^{b}_j$ equals $\pm1$ with equal probability. To retain validity of the test regardless of the value of $B$, we require that $g^1=\iota$. We note, however, that the power of the test may still depend on $B$. For this reason, we implement Algorithm (ref) with $\hat{\mathbf G}$ replacing $\mathbf G$ everywhere and set $B=1000$ (or any other reasonably large number chosen by the analyst).

Step 4 requires computing the $1-\alpha$ quantile of $\{T_n(g): g\in\mathbf G\}$, which can be typically obtained by sorting the values of $\{T_n(g): g\in\mathbf G\}$ and then taking the $\lceil |\mathbf G|(1-\alpha) \rceil^{\rm th}$ highest element in the ordered list. Thus, if we denote the ordered values of $\{T_n(g): g\in\mathbf G\}$ by $$ T^{(1)}_n\le T^{(2)}_n \le \cdots \le T^{(B)}_n~,$$ then we may define $\hat c_n(1-\alpha)$ in (ref) as $\hat c_n(1-\alpha)= T^{(\lceil |\mathbf G|(1-\alpha) \rceil)}$. This representation suggests that the test may have trivial power for very low values of $q$. For example, when $\alpha=10\%$, this problem arises if $q\le 4$. For $q=5$ the test already has non-trivial power and is only slightly conservative under the null. Similarly, when $\alpha=5\%$ the test has non-trivial power for any $q\ge 6$.

Step 5 is straightforward and it provides both the test $\phi_n$ and the $p$-value $\hat p_n$. Each of these correspond to the non-randomized version of ARTs as opposed to their randomized counterparts canay/etal:17 since practitioners often prefer tests that do not involve exogenous randomness. In any case, the differences between the randomized and non-randomized versions of the test have been found to be minimal in simulations canay/etal:17.

How to compute confidence intervals

We now discuss how to compute confidence intervals for the parameter $c^\prime \beta$ by developing a novel algorithm that exploits the properties derived in Section (ref). As before, a particularly important case is when $c$ selects the $\ell$th component of $\beta$ and then the confidence set is simply a confidence interval for $\beta_{\ell}$. Conceptually we can simply form the confidence set by collecting all values of $c^\prime \beta$ that cannot be rejected by our test at level $\alpha$. That is, for the test $\phi_n$ in (ref) we define

equation[equation omitted — 135 chars of source]

In an asymptotic framework where $n\to \infty$ while $q$ remains fixed, canay/etal:17 show that $\phi_n$ is asymptotically level $\alpha$ under $H_0$. It follows from that result that, by construction, $C_n$ covers $c^\prime \beta$ with probability at least equal to $1-\alpha$ asymptotically. In Section (ref) we show that $C_n$ is indeed a closed interval in $\mathbf{R}$ and so it takes the form

equation[equation omitted — 67 chars of source]

where $\lambda_{l}$ is the smallest value of $\lambda$ that cannot be rejected by $\phi_n$ and $\lambda_{u}$ is the largest value of $\lambda$ that cannot be rejected by $\phi_n$. The analysis in Section (ref) also reveals that $\lambda_{l}$ and $\lambda_{u}$ admit simple closed-form representations that we exploit to develop Algorithm (ref) below.

algorithm[algorithm omitted — 3,024 chars of source]

Algorithm (ref) requires three steps that are straightforward to compute and that exploit the results in Section (ref). We refer the reader to that section for the details on why $\lambda_l$ and $\lambda_u$ admit the expressions in (ref) and (ref), respectively.

Three results on implementation of ARTs

Before we review the main requirement underlying ARTs, we present three properties related to the implementation of ARTs that we believe practitioners should be aware of and that are novel to this paper. The first property establishes a connection between the implementation of ARTs as described in Algorithm (ref) and an alternative implementation based on weighted scores. The second property establishes the numerical equivalence of ARTs for the null in (ref) when the test statistics in (ref) is replaced by its studentized version. The third and final result shows that ARTs confidence set for $c'\beta$ is indeed a closed interval in $\mathbf{R}$ and provides a representation for the upper and lower bounds of the interval that lead to Algorithm (ref).

Equivalence with weighted scores

It turns out that ARTs can be implemented by an algorithm that does not involve estimating the parameter $\beta$ within each cluster. This alternative algorithm involves replacing Steps 1 and 2 in Algorithm (ref) by the two alternative steps described in Algorithm (ref) below, while keeping Steps 3 to 5 unaffected.

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

Note that Steps 3-5 remain unchanged given the alternative definition of $S_{n,j}$ in Step 2$^\prime$. When it comes to Steps 1 and 2, there are two differences worth discussing. The first difference is that Step 1$^\prime$ requires a single full-sample restricted least squares estimator of $\beta$ as opposed to the $q$ cluster-by-cluster estimators in Step 1 of Algorithm (ref). The second difference is that Step 2$^\prime$ is based on within-cluster weighted scores as opposed to the centered within-cluster estimates of $\beta$ in Step 2 of Algorithm (ref). Interestingly, these two implementations are numerically equivalent and so implementing ARTs via Algorithm (ref) or Algorithm (ref) leads to identical results. To see this formally, it is enough to show that $S_{n,j}$ as defined in (ref) and (ref) are the same using the following argument. For each $j\in J$,

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

where the fourth equality follows by adding and subtracting $\sqrt{n_j}c'\beta$ and the last equality holds because $c'\hat \beta_{n}^{\rm r}=c'\beta=\lambda$ under the null hypothesis in (ref). It thus follows that $S_{n,j}$ in (ref) and in (ref) are identical and so ARTs can be alternatively implemented via Algorithm (ref) or (ref). The following lemma summarizes our discussion above:

lemmaLet $\hat\Omega_{n,j}$ in (ref) be full rank for each $j\in J$. Denote by $C_n$ a confidence interval for $c'\beta$ computed using Algorithm (ref) and by $C_n^\prime$ a confidence interval for $c'\beta$ computed using Algorithm (ref). Then $C_n = C_n^\prime$.

Equivalence with studentized version of the t-statistic

The ART defined in (ref) of Algorithm (ref) is based on the unstudentized test statistic $T_n$ defined in (ref). It may perhaps appear more desirable to instead consider the studentized version of this test statistic as studentization commonly improves performance in a variety of other settings. Here, we prove that this is not the case for ARTs when the null hypothesis is the one in (ref) and that both versions of the test statistic lead to numerically identical results.

To see this, start by defining the studentized version of the test statistic in (ref) as $T_n^{\rm s} \equiv T_n^{\rm s}(\iota)$, where for each $g\in \mathbf G$,

equation[equation omitted — 291 chars of source]

Then note that

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

where $V_n \equiv \frac{1}{q} \sum_{j=1}^q g_j^2S^2_{n,j}$ does not depend on $g$ as $g_j^2=1$ for all $j\in J$. It follows that we can write the studentized test statistic as

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

Since the function $x \mapsto \frac{x}{\sqrt{1 - x^2}}$ is strictly increasing for $x \in [0,1)$, it follows that $T_n^{\rm s}(g)$ is a strictly monotonic transformation of $T_n(g)$ for each $g\in \mathbf G$. We conclude that $I\{T_n(g)\ge T_n(\iota)\}=I\{T_n^{\rm s}(g)\ge T_n^{\rm s}(\iota)\}$ for all $g\in \mathbf G$ and so the ART based on $T_n(g)$ and $T_n^{\rm s}(g)$ are identical. This discussion is summarized in the following lemma:

lemmaLet $\hat\Omega_{n,j}$ in (ref) be full rank for each $j\in J$. Denote by $C_n$ a confidence interval for $c'\beta$ computed using Algorithm (ref) and by $C_n^\prime$ a confidence interval for $c'\beta$ computed using Algorithm (ref) with $T_n^{\rm s}$ in place of $T_n$ and $T_n^{\rm s}(g)$ in place of $T_n(g)$. Here, $T_n^{\rm s}(g)$ is given by (ref) and $T_n^{\rm s}$ is understood to be $T_n^{\rm s}(\iota)$, where $\iota$ is the identity transformation. Then, $C_n = C_n^\prime$.

Convexity of the confidence intervals

The ART-based confidence intervals for $c'\beta$ defined in (ref) can be computed by test inversion. From a computational standpoint, however, computing confidence sets by test inversion may be cumbersome and the resulting set may not even be an interval. That is, it may not be closed and convex. In this section we prove that this is not a concern for ART-based confidence intervals for $c'\beta$ and so such confidence intervals could be easily computed by a standard bisection algorithm. In fact, our results go even further. We derive closed form expressions for the lower and upper bounds of the confidence interval that imply that computing ART-based confidence intervals for $c'\beta$ is straightforward from a computational standpoint. In order to derive these results, we slightly change our notation to make explicit the dependence on $\lambda$ of each of the elements entering the test in (ref). To this end, let

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

and note that $T_n = T_n(\iota, \lambda)$. Using this notation, we can re-write the confidence interval in (ref) as

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

which is simply the values of $\lambda$ for which the $p$-value of the test, as defined in (ref), is not below $\alpha$. In order to show that this confidence set is a closed interval, we claim that the $p$-value

equation[equation omitted — 166 chars of source]

is equal to $1$ for $\lambda_0 \equiv b(\iota)/a(\iota)$, monotonically increasing for any $\lambda<\lambda_0$, and monotonically decreasing for any $\lambda>\lambda_0$. The next lemma formalizes this result.

lemmaLet $\hat\Omega_{n,j}$ in (ref) be full rank for each $j\in J$. Let $a(g)$, $b(g)$, and $\lambda_0$ be defined as in (ref). The p-value in (ref) equals \begin{equation} \hat p_n(\lambda) = \begin{cases} \frac{1}{|\mathbf{G}|}\sum_{g \in \mathbf{G}} I \left\{\lambda\ge \lambda_{l}(g)\right\} & for \lambda<\lambda_0\\[3pt] 1 & for \lambda=\lambda_0\\[3pt] \frac{1}{|\mathbf{G}|}\sum_{g \in \mathbf{G}} I \left\{\lambda\le \lambda_{u}(g)\right\}& for \lambda>\lambda_0 \end{cases} , \end{equation} where $\{\lambda_u(g):g\in \mathbf G\}$ and $\{\lambda_l(g):g\in \mathbf G\}$ are defined in Algorithm (ref).
proofIt is useful to re-write $T_n(g,\lambda)$ in terms of $a(g)$ and $b(g)$. To this end, note that \begin{align} T_n(g,\lambda) &\equiv \Big| \frac{1}{q}\sum_{j=1}^q g_j S_{n,j}(\lambda) \Big|\notag = \Big| \frac{1}{q}\sum_{j=1}^q g_j \sqrt{n_j}c'\hat\beta_{n,j}-\lambda\frac{1}{q}\sum_{j=1}^q g_j \sqrt{n_j} \Big|\notag \\ &= |b(g) - \lambda a(g)| . \end{align} Given $g\in \mathbf G$ and $a(g)\ne 0$, $T_n(g,\lambda)$ is a “V-shaped” function of $\lambda$ taking the value $0$ at $\frac{b(g)}{a(g)}$ and with slope $-|a(g)|$ for all $\lambda<\frac{b(g)}{a(g)}$ and slope $|a(g)|\leq a(\iota) $ for all $\lambda>\frac{b(g)}{a(g)}$. Figure (ref) illustrates this for three values of $g$. First, note that $I \left\{T_n(g, \lambda_0)\ge T_n(\iota,\lambda_0) \right\}=I \{T_n(g, \lambda_0)\ge 0 \}=1$ for all $g\in\mathbf G$ and so it follows immediately that $\hat p_n(\lambda_0)=1$. Second, restrict attention to the set $\Lambda^{+}\equiv \{\lambda \in \mathbf R: \lambda> \lambda_0\}$ where $T_n(\iota,\lambda)$ is linearly increasing. In order to prove that $\hat p_n(\lambda)$ takes the form in (ref) we prove that $I \left\{T_n(g, \lambda)\ge T_n(\iota,\lambda) \right\}=I \left\{\lambda\le \lambda_{u}(g)\right\}$ for each $g\in\mathbf G$ by dividing the argument into three cases. Case 1: Consider $g\in \mathbf G$ such that $a(g)\ne 0$ and $|a(g)|\ne a(\iota)$. Since $|a(g)|< a(\iota)$, it follows that $T_n(g, \lambda)$ and $T_n(\iota, \lambda)$ intersect only once on $\Lambda^{+}$ and this holds regardless of whether $\frac{b(g)}{a(g)}< \lambda_0$ or $\frac{b(g)}{a(g)}\ge \lambda_0$ (see Figure (ref) for a graphical illustration of each of these cases). Denote the intersection point by $\lambda_{u}(g)$ and note that $T_n(g,\lambda)\ge T_n(\iota,\lambda)$ for all $\lambda_0< \lambda\le \lambda_{u}(g)$ and $T_n(g,\lambda)<T_n(\iota,\lambda)$ for all $\lambda>\lambda_u(g)$. Conclude that on $\Lambda^{+}$, \begin{equation} I\{T_n(g,\lambda)\ge T_n(\iota,\lambda)\} = I\{\lambda\le \lambda_{u}(g)\} . \end{equation} Simple algebra shows that the intersection point $\lambda_u(g)$ takes the form in (ref). Case 2: Consider $g\in \mathbf G$ such that $a(g) = 0$. Note that $b(\iota) - \lambda a(\iota)<0$ for $\lambda\in\Lambda^{+}$. It thus follows that for $\lambda\in \Lambda^{+}$, \begin{equation*} I\{T_n(g,\lambda)\ge T_n(\iota, \lambda) \} = I\{|b(g)|\ge |b(\iota) - \lambda a(\iota)| \} = I\left\{ \lambda \le \frac{b(\iota)}{a(\iota)} + \frac{|b(g)|}{a(\iota)} \right\} , \end{equation*} and so (ref) holds in this case with $\lambda_{u}(g) = \frac{b(\iota)}{a(\iota)} + \frac{|b(g)|}{a(\iota)} $, as defined in (ref). Case 3: Consider $g\in \mathbf G$ such that $|a(g)|=a(\iota)$ and so $g=\pm \iota$. If $g=\iota$, $I\{T_n(g,\lambda)\ge T_n(\iota, \lambda) \}=1$ for all $\lambda\in \mathbf R$. We conclude that (ref) holds with $\lambda_u(g)=\infty$. If $g=-\iota$, then we have that $a(-\iota)=-a(\iota)$ and $b(-\iota)=-b(\iota)$ so that $\frac{b(-\iota)}{a(-\iota)} =\lambda_0$ and again $I\left\{ T_n(-\iota, \lambda) \ge T_n(\iota, \lambda) \right\} = 1$ for all $\lambda\in \mathbf R$. We conclude that (ref) holds with $\lambda_u(g)=\infty$, as defined in (ref). This completes the proof of (ref) for the case $\lambda \in \Lambda^{+}$. Finally, the construction for $\lambda\in \Lambda^{-}\equiv \{\lambda \in \mathbf R: \lambda < \lambda_0\}$ parallels the one for $\lambda \in \Lambda^{+}$ so we omit the arguments here. Putting all the cases together, (ref) follows and this completes the proof.
figure[figure omitted — 4,321 chars of source]

Figure (ref) illustrates the $p$-value in (ref) as a function of $\lambda$ for the groups in Figure (ref). Since $ \hat p_n(\lambda)$ is right continuous and increasing for $\lambda<\lambda_0$, we can define $\lambda_l$ as the smallest value of $\lambda$ for which $ \hat p_n(\lambda)\ge \alpha$. Such value exists and is unique. Similar, since $ \hat p_n(\lambda)$ is left continuous and decreasing for $\lambda>\lambda_0$, we can define $\lambda_u$ as the largest value of $\lambda$ for which $ \hat p_n(\lambda)\ge \alpha$. Such value exists and is again unique. This argument leads to the representation of $C_n$ in (ref), showing that ART-based confidence intervals for $c'\beta$ are indeed intervals in $\mathbf R$. Furthermore, note that (ref) implies that the smallest value of $\lambda$ for which $ \hat p_n(\lambda)\ge \alpha$ can be defined as

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

which is just the definition of the $\alpha$ quantile of $\lambda_l(g)$, as defined in Algorithm (ref). A similar result holds for $\lambda_u$ and so $C_n$ can be computed in closed form by Algorithm (ref).

figure[figure omitted — 3,828 chars of source]

What we need for ARTs to work

The main requirement underlying ARTs is Assumption 3.1 in canay/etal:17. This assumption guarantees that the test delivers rejection probabilities under the null hypothesis that are close to the nominal level $\alpha$ in an asymptotic framework where $n\to \infty$ and $q$ remains fixed. In the context of the linear model in (ref), this translates into the following two conditions summarized in Assumption (ref) below.

assumptionLet $\{\hat\beta_{n,j}:j\in J\}$ be the cluster-by-cluster estimators of $\beta$ defined in Algorithm (ref). Assume that: \begin{enumerate}[(a)] • $\{\hat\beta_{n,j}:j\in J\}$ jointly converge in distribution at some (possibly unknown) rate; i.e., \begin{equation} \left( \begin{array}{c} a_{n,1}(\hat\beta_{n,1}-\beta)\\ \vdots\\ a_{n,q}(\hat\beta_{n,q}-\beta)\\ \end{array} \right)\overset{d}{\to} \left( \begin{array}{c} S_1\\ \vdots\\ S_q\\ \end{array}\right) \end{equation} for a sequences $a_{n,j}\to \infty$ and random variables $(S_1,\dots,S_q)'$. • The limiting random variables $(S_1,\dots,S_q)'$ are invariant to sign changes, i.e., \begin{equation} (g_1S_1,\dots,g_qS_q) \overset{d}{=} (S_1,\dots,S_q) , \end{equation} for any $g$ in $\mathbf G$, where $\mathbf G$ is defined in Step 4 of Algorithm (ref). \end{enumerate}

Condition (ref) holds, for example, when $Z_{i,j}$ and $\epsilon_{i,j}$ are uncorrelated and the analyst assumes some form of weak dependence within clusters that permits the application of an appropriate central limit theorem. In such a case, (ref) typically holds with $a_{n,j}=\sqrt{n_j}$ and each $S_j$ being a mean-zero normal random variable. In fact, under the commonly used assumption of independent clusters, it also follows that $S_j \perp \!\!\!\!\perp S_{j'}$ for any $j\ne j'$. In this case the normally distributed random variables may not be identically distributed but are indeed independent. Condition (ref), in turn, requires each $S_j$ to be symmetrically distributed around zero and independent of each other. This is immediately satisfied when each $S_j$ is a mean-zero normal random variable and clusters are independent. Importantly, these assumptions allow for the normally distributed random variables to have different variances across clusters; a type of heterogeneity not allowed by the cluster wild bootstrap approach popularized by cameron2008bootstrap and later studied formally by canay/santos/shaikh:20.

remark\rm The asymptotic normality in (ref) arises frequently in applications, but is not necessary for the validity of ARTs. All that is required is that the estimators $\{a_{n,j}(\hat\beta_{n,j}-\beta):j\in J\}$ have a limiting distribution that is the product of $q$ distributions that are symmetric about zero. This may even hold in cases where the estimators have infinite variances or are inconsistent. See canay/etal:17 for additional discussion on this point. \rule{2mm}{2mm}
remark\rm It is worthwhile to contrast the requirements of Assumption (ref) with those of “classical” methods, such as those described in liang1986longitudinal. These latter methods permit arbitrary dependence within each cluster, but require the size of the clusters to be small and the number of clusters to be large. As described above, Assumption (ref)(a), on the other hand, permits the number of clusters to be small, but requires the size of the clusters to be large and weak dependence within each cluster. We emphasize, however, that these restrictions are commonly employed in establishing the validity of other methods in settings with a small number of clusters, including, for example, the $t$-test approach by ibragimov2010t and the wild bootstrap canay/santos/shaikh:20. \rule{2mm}{2mm}
remark\rm We focus our exposition on the case where $Z_{i,j}$ is exogenous but we emphasize that the conditions in (ref) and (ref) typically hold in instrumental variable (IV) models. Accommodating IV to ARTs then only requires modifying Step 1 in Algorithm (ref) so that the least squares regression is replaced with the appropriate IV regression. Steps 2-6 remain unaffected. \rule{2mm}{2mm}

An implicit requirement behind ARTs that deserves further comments lies in Step 1 of Algorithm (ref), which requires that the analyst runs cluster-by-cluster regressions. This step implicitly assumes that the parameter $\beta$ is identified within each cluster. In practice, this means that the matrix $\hat \Omega_{n,j}$ in (ref) must be invertible for each $j\in J$ and hence the same requirement applies to Algorithm (ref). This restriction may be substantially important in some applications and so here we discuss common ways in which the problem may manifest and two alternative remedies.

One case in which running least squares cluster-by-cluster is not feasible is when the coefficient of interest is associated with a variable that only varies across clusters. For example, consider the model in (ref) and partition $Z_{i,j}$ into a constant term, a scalar variable that only varies across clusters, $Z^{(1)}_j$, and another variable that varies across and within clusters, $Z^{(2)}_{i,j}$. That is,

equation[equation omitted — 121 chars of source]

where the analysts' interest lies in the coefficient $\beta_1$, i.e., $c'\beta = \beta_1$. Clearly, the regression in Step 1 of Algorithm (ref) would not separately identify $\beta_0$ and $\beta_1$ as $Z^{(1)}_{j}$ is perfectly colinear with the constant term. The matrix $\hat \Omega_{n,j}$ in (ref) is simply singular. This situation arises, for example, in the empirical application considered by canay/etal:17supp where $j\in J$ indexes schools and the variable of interest is a treatment indicator at the school level. A natural remedy in a situation like this is clustering more coarsely (e.g., by combining clusters) to obtain variation within the re-defined clusters. This is possible for ARTs since the validity of the method does not rely on having a large number of clusters and thus it can afford to work with coarser clustering. In fact, in certain settings combining clusters may be quite natural. For example, canay/etal:17supp re-defined clusters as “pairs” of schools (as opposed to just schools) given that the treatment assignment mechanism of the experiment was a matched pairs design and so the pairs used at the randomization stage represented natural groupings. In other settings where it is less clear how to group clusters, any grouping that satisfies the requisite identification condition leads to a valid test, but it may be further desirable to combine such tests to limit concerns about “data snooping” across groupings. To this end, results in diciccio2020exact on combining tests may be relevant.

remark\rm A quick inspection of (ref) may lead the analyst to believe there is a workaround that does not involve combining clusters if one instead uses some estimator of $\beta_0$ from a full sample regression. For example, the full sample least squares estimator $\hat\beta_{n,0}$ from the regression in (ref). Then, assuming for simplicity that $Z_j^{(1)}\ne 0$ for all $j\in J$, one may consider modifying Step 1 in Algorithm (ref) by running a regression of $Y_{i,j}$ on an intercept and $Z^{(2)}_{i,j}$ (not including $Z_j^{(1)}$) and then redefining $\{\hat \beta_{n,j}:j \in J\}$ as the difference between the within cluster intercept estimates, $\hat\beta_{j,0}$ and the full sample estimate $\hat\beta_{n,0}$, i.e., $\hat \beta_{n,j}=\hat\beta_{j,0}-\hat\beta_{n,0}$. Such strategies unfortunately introduce dependence between the $q$ estimators of $\beta$ (as they all depend on $\hat\beta_{n,0}$) and thus end up violating one of the two main conditions needed for ARTs to be asymptotically valid; mainly condition (ref). \rule{2mm}{2mm}

Another case where the lack of identification within cluster may manifest is when the variable of interest actually varies within clusters but the model specification involves other variables that are collinear with some other variable (including the variable of interest or the constant term) within clusters. For example, consider the model in (ref) where instead of individuals indexed by $i\in I_{n,j}$, units within cluster are indexed over time $t\in T$. Partition $Z_{j,t}$ into the variable of interest, $Z^{(1)}_{j,t}$, and time fixed effects $\delta_t$. That is,

equation[equation omitted — 145 chars of source]

It then follows that, within each cluster $j\in J$, the time fixed effect $\delta_t$ absorbs all the variation in $Z^{(1)}_{j,t}$ and so $\beta_1$ is not identified. In cases like this the analyst could again combine clusters to obtain variation within the re-defined clusters. An alternative remedy is to change the specification by, for example, replacing the time fixed effect with a cluster-specific time trend. Such specification is more restrictive than the time fixed effect in the sense that it imposes a linear trend but, at the same time, is more general as it allows for heterogeneity across clusters in the linear trend. We illustrate this approach in the application we consider in Section (ref).

The need to identify $\beta$ within each cluster is in our view the main limitation of ARTs, but a limitation that needs to be dealt with in certain settings. One may then wonder why not simply use some other inference method that is valid when the number of clusters is small and that does not rely on estimating $\beta$ cluster-by-cluster. Perhaps the most popular approach in that category is the cluster wild bootstrap popularized by cameron2008bootstrap and recently studied formally by canay/santos/shaikh:20. While not having to estimate $\beta$ within each cluster represents an advantage over ARTs, this additional flexibility comes at a cost in terms of the degree of heterogeneity that the model can deal with. In particular, the results in canay/santos/shaikh:20 show that the cluster wild bootstrap is expected to work well in settings with a small number of clusters as long as the clusters are “homogeneous,” in a sense made precise in canay/santos/shaikh:20. Intuitively, it is required that the variance covariance matrix $\hat \Omega_{n,j}$ defined in (ref) is the same across clusters (up to scalar multiplication). Such stringent homogeneity condition is not required for ARTs to work well, as the method allows clusters to be arbitrarily heterogeneous as long as $\hat \Omega_{n,j}$ is invertible for $j\in J$.

remark\rm For ease of exposition, we have written the requirement in (ref) in terms of the differences $\hat \beta_{n,j} - \beta$, but it is possible to replace it with the differences $c’\hat \beta_{n,j} - c’\beta$ (or $R\hat \beta_{n,j} - R\beta$, depending on the null hypotheses of interest). In most cases, re-writing the condition in this way is not useful, but it is in cases where $c’\beta$ is identified within each cluster while $\beta$ is not. For example, consider the model in (ref) when the coefficient of interest is $\beta_2$ as opposed to $\beta_1$, i.e, $c'\beta = \beta_2$. In that case the entire term $\beta_0+Z_j^{(1)}\beta_1$ may be absorbed into a cluster-specific intercept without affecting the identification and estimation of $c'\beta = \beta_2$ within each cluster. \rule{2mm}{2mm}

Empirical applications

In this section we apply ARTs as described in Algorithm (ref) and ART-based confidence intervals as described in Algorithm (ref) in the context of two distinct empirical applications. The R and Stata packages and codes required to replicate the results in this section are available as part of the online supplemental material.

Meng, Qian and Yared (2015)

mengQY2015 argue that China's Great Famine, from 1959 to 1961, was the result of an inflexible food procurement policy by the central government. To make this point, they show that food production and mortality become positively correlated during the time of famine, when this coefficient is otherwise negative or not significantly different from 0 in normal times.

MQY consider the following regression,

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

where $j$ indexes provinces (ranging from 1 to 19) and $t$ indexes years (ranging from 1953 to 1982). Here,

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

In this application the level of clustering is a province, and so in order to apply ARTs as described in Section (ref), one needs to estimate $\beta=(\beta_1,\beta_2)'$ and $\delta_t$ province-by-province. This illustrates one of the situations where including time fixed effects province-by-province is infeasible for the implementation of ARTs, given that the only source of remaining variation within a province is indeed time. The second identification problem described in Section (ref) then arises. As we discussed in that section, one way to deal with this issue consists of replacing the time fixed effects with a cluster-specific time trend, i.e., in Step 1 of Algorithm (ref) estimate

equation[equation omitted — 130 chars of source]

We will refer to this as Analysis \#1. In addition, we also consider the following alternative specifications studied by MQY:

itemize• Analysis \#2: Repeating Analysis \#1 using only data between 1953 and 1965. • Analysis \#3: Repeating Analysis \#1 using four additional autonomous provinces. • Analysis \#4: Repeating Analysis \#2 using four additional autonomous provinces. • Analysis \#5: Repeating Analysis \#1 using actual rather than constructed grain production. • Analysis \#6: Repeating Analysis \#2 using actual rather than constructed grain production.

As with Analysis \#1, the above analyses differ from their MQY counterparts only in that a linear time trend $\gamma_j t$ replaces time fixed effects $\delta_t$. Table (ref) summarizes the number of clusters and the number of observations for each of these analyses. We caution, however, that in this application, in addition to the number of clusters being small, the number of observations within each cluster may also be small. See Remark (ref) for further discussion in relation to Assumption (ref)(a).

table[table omitted — 612 chars of source]

mengQY2015 consider the following two null hypotheses of interest,

equation[equation omitted — 119 chars of source]

In Table (ref) we replicate the main table in mengQY2015 using cluster robust standard errors (CCE) and also include the results associated with ARTs for both $H_0^{(1)}$ and $H_0^{(2)}$ in (ref). For $H_0^{(1)}$ we report $p$-values and 95% confidence intervals, while for $H_0^{(2)}$ we just report $p$-values following MQY. The authors note in footnote 33 that using the cluster wild bootstrap led to similar results as those presented in their main table so we do not include cluster wild bootstrap results here either.

table[table omitted — 2,113 chars of source]

We comment on the following main features of Table (ref):

enumerate• For the null hypothesis $H_0^{(1)}$ associated with the parameter $\beta_1$, the ART $p$-values are of comparable magnitude to traditional CCE $p$-values. Similarly, ART-based confidence intervals are of roughly the same length as those obtained based on CCE although the ART-based confidence intervals do not contain the LS estimates. This is because ART-based confidence intervals are centered around the mean of the province-by-province estimates, which may not necessarily be equal to the full sample LS estimate of $\beta_1$. • For the null hypothesis $H_0^{(2)}$ associated with the parameter $\beta_1 + \beta_2$, the ART $p$-value is sometimes higher and sometimes lower than the CCE $p$-value depending on the specification. Given the relatively small number of clusters in this application, the ART $p$-values are likely to be more reliable than those associated with CCE as CCE is known to perform poorly when the number of clusters is not sufficiently large.

Munyo and Rossi (2015)

munyo2015first study criminal recidivism of former prisoners by looking at the relationship between the number of inmates released from incarceration on a given day and the number of offenses committed on the same day. They claim that the liquidity constraints that inmates face on the day of release increase the likelihood of recidivism on the same day. Using data of 2631 days between January 1st 2004 and March 15 2011 collected from the criminal incidents reports in Montevideo in Uruguay, they estimate the following linear model by least squares

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

where $t$ indexes days and

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

We refer to this as Analysis \#1. munyo2015first additionally consider the following four analyses:

itemize• Analysis \#2: $Z_{t}$ includes a daily trend in place of a yearly trend. • Analysis \#3: $Z_{t}$ includes a monthly trend in place of a yearly trend. • Analysis \#4: $Z_{t}$ includes an intra-month daily trend, month- and year- level fixed effects and their interactions in place of a yearly trend. • Analysis \#5: $Z_{t}$ includes month- and year- level fixed effects and their interactions in place of a yearly trend.

Analysis \#5 is their preferred specification. munyo2015first report the results of these analyses in Table 2 in their paper. They report least squares estimates of $\beta$ with Newey-West heteroskedasticity-autocorrelation-consistent (HAC) standard errors. In addition, they report ART $p$-values as described in Algorithm (ref) for the null hypothesis that $H_0: c'\beta=0$ as in (ref), where $c$ selects the coefficient on the total number of inmates released on day $t$.

In this application the level of clustering is not naturally determined by the data, but pseudo-clusters may be formed using blocks of consecutive observations under the assumption of weak temporal dependence. In order to apply ARTs as described in Algorithm (ref) we then form $q$ pseudo-clusters by dividing the data into $q$ consecutive blocks of size $b_n = \lfloor n/q \rfloor$ where $n=2631$ is the number of total observations. More concretely, we define the $j$th pseudo-cluster as

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

and let the last $q$th pseudo-cluster contain all the remaining $n- b_n(q-1)$ observations. Note that in this application the number of pseudo-clusters $q$ is a tuning parameter that the analyst must specify. munyo2015first set $q=10$. We repeat their analyses with alternative values of $q$ and investigate how sensitive the results are to this choice. The relevant cluster information is given in Table (ref).

table[table omitted — 309 chars of source]
table[table omitted — 2,108 chars of source]

Table (ref) shows LS estimates of $\beta$, $p$-values for the hypothesis in (ref), and $95\%$ confidence intervals for each analysis. Following munyo2015first, we report results based on HAC standard errors. The table also shows ART $p$-values as described in Algorithm (ref) and ART-based 95% confidence intervals as described in Algorithm (ref) for $q=8$, $q=10$, and $q=16$.

We summarize the main findings of the results in Table (ref) as follows:

enumerate• The choice of $q$ is important for the results of ARTs but currently there is no theory developed to choose this tuning parameter according to some data dependent criteria. The smaller $q$ is, the more observations are available within each cluster. Having more observations per cluster is important for one of the requirements behind ARTs, mainly (ref). A small value of $q$, however, tends to affect the power of ARTs despite not really affecting the control of the rejection probability under the null hypothesis. This feature can be seen in Table (ref), where ARTs $p$-values are decreasing in $q$ across different specifications. In this application, where there are still over a hundred observations when $q=16$, a larger value of $q$ like $q=10$ or $q=16$ may be preferable to smaller values, like $q=8$, based on power considerations. Note, however, that except in Analyses \#4--5, where the choice of $q$ determines whether the null hypothesis is rejected at a given significance level, the results for Analyses \#1--3 are in all agreement at a $5\%$ level. • Overall, the test results based on standard $t$-test with HAC standard errors are consistent to those of ARTs when $q=16$. Both methods reject the null hypothesis $H_0:c'\beta=0$ at a $10\%$ nominal level across different specifications. The results support the authors' argument that the release of inmates from incarceration increase the chance of re-offenses on the day of release.

Computational gains of the new algorithm

Tables (ref) and (ref) report four alternative ways to compute ART-based confidence intervals in the two empirical applications we consider in this paper; mengQY2015 and munyo2015first. The first alternative is to compute the confidence intervals by a simple grid search algorithm. The second alternative involves a bi-section algorithm. We implement both of these methods using a studentized and an unstudentized test statistic to illustrate the result in Section (ref). The last alternative is to simply use Algorithm (ref), as reported in Sections (ref) and (ref). In each case, we also report computational times to illustrate the computational advantages of the algorithm we propose in this paper. The R and Stata codes required to replicate the results in this section are available as part of the online supplemental material.

Starting from Table (ref), we see that grid search take a significant amount of time to compute. Our convexity result (Lemma (ref)) facilitates the use of the bisection method, cutting implementation time by a factor of over 50. Moving from the bisection method to Algorithm (ref) further leads to a speed up of at least 2 times. A similar pattern emerges in Table (ref). Furthermore, comparing specification with $q=8$ that with $q=16$, the speed advantage of our method becomes far starker. For $q=16$, grid search takes almost 100 times as long as the bisection method. The bisection method, meanwhile, takes close to 10 times as long as Algorithm (ref).

table[table omitted — 1,849 chars of source]
table[table omitted — 3,867 chars of source]

Concluding remarks

The goal of this paper is to make the general theory developed in canay/etal:17 more accessible by providing a step-by-step algorithmic description of how to implement the test and construct confidence intervals in linear regression models with clustered data, as well as clarifying the main requirements and limitations of the approach. The main two takeaways are the following. First, ARTs-based confidence intervals for scalar parameters in linear regression models can be characterized in closed form and thus are straightforward to implement in practice. Algorithms (ref) and (ref) provide a clear explanation of how to apply ARTs in linear models, and the companion {\tt Stata} and {\tt R} packages available as part of the supplemental material are intended to facilitate doing so. Second, our discussion on the main requirements behind ARTs hopefully show that understanding the trade-offs between ARTs and other popular alternatives for inference with a small number of clusters, like the cluster wild bootstrap, is fundamental for practitioners to choose a method that aligns well with the features of their application. In particular, while ARTs essentially demand that the parameter of interest is suitably estimable cluster-by-cluster without imposing restrictions on the degree of heterogeneity across clusters, the cluster wild bootstrap requires the clusters to be sufficiently homogeneous canay/santos/shaikh:20 without demanding identification of the parameter of interest cluster-by-cluster.

\setcounter{lemma}{0} \setcounter{theorem}{0} \setcounter{corollary}{0} \setcounter{equation}{0} \setcounter{remark}{0}