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.
74,025 characters · 20 sections · 52 citation commands
Combining Clusters for the Approximate Randomization Test
\fontfamily{ppl}
It is common to have clustered data in economics where observations within each cluster are correlated with each other. The recent approximate randomization test proposed by \citet*[CRS]{canayetal2017ecta} imposes weak assumptions on the dependence structure and can be used to conduct inference when there is a small number of clusters. However, a main requirement in CRS is that the target parameter has to be identified within each cluster. If this requirement is not satisfied, researchers can combine clusters in a way that ensures the target parameter can be identified within each combined cluster in order to apply CRS. Any method of combining clusters that solves the identification issue leads to a valid test. However, this creates ambiguity on how clusters should be combined. This paper addresses this issue and provides guidance on how to combine clusters.
The issue of not being able to identify the target parameter within each cluster is common in empirical settings. Difference-in-differences (DID) is a leading example of this identification issue within clusters. This is primarily because the treatment indicator equals zero in the control clusters. Some examples include bloometal2012qje for firm-level treatment, burdelinden2013aejae for village-level treatment, dinceccokatz2016ej for country-level treatment, goodman2016restat for state-level treatment, and alsanclaudia2019jpe for municipality-level treatment. In addition, if each cluster contains only one unit with multiple periods, then the identification issue arises in all clusters when cluster-by-cluster regressions are used. This is because there are more coefficients than observations per cluster when time fixed effects are included.
This paper makes two contributions to the literature. First, I analyze the local asymptotic power of CRS. canayetal2017ecta showed that their test controls size but did not study local asymptotic power. Analyzing the local asymptotic power of CRS is non-trivial because the number of clusters is fixed and does not go to infinity. Hence, standard large-sample analysis of local power for randomization tests, such as hoeffding1952ams, cannot be applied.
Second, I develop computationally efficient procedures to combine clusters by maximizing local asymptotic power. The procedures depend on the chosen significance level. I show that the procedure can be written as a binary linear program when the significance level is “small.” For other significance levels, I develop a heuristic to combine clusters. The procedures can be solved quickly and can help eliminate randomness in deciding how to combine the clusters. In the simulation study and empirical exercise, I find that the procedures perform well in various settings.
Similar to CRS, this paper is related to the literature on inference with a small number of clusters. bertrandetal2004qje pointed out the issue with standard inference when there is a small number of clusters. Different tests have been proposed in this literature apart from CRS. Some recent methods include besteretal2011joe, ibraimovmuller2010jbes, ibraimovmuller2016restat, and hagemann2022wp. The Wild cluster bootstrap popularized by cameronetal2008restat has been shown to be valid under strong homogeneity conditions when the number of clusters is fixed by canaysantosshaikh2021restat. The Wild cluster bootstrap has also been studied in other settings, such as mackinnonwebb2023jae, djogbenouetal2019joe, and webb2023cje. caietal2023jem provide a users-guide for CRS. cai2024wp studies finite sample properties of the sign test. Recently, caoetal2022wp proposed a data-driven procedure to partition observations into clusters and conduct valid inference. Their procedure focuses on spatially indexed data and assumes there exists a dissimilarity measure between observations. This paper is different in that I take the clusters as given and focus on combining the clusters when the target parameter cannot be identified within each cluster. This paper follows the above papers and uses the model-based approach in clustering. See abadieetal2022qje for a design-based analysis. For surveys on the literature of conducting inference with a fixed number of clusters, see, for instance, cameronmiller2015jhr, conleyetal2018jar, and mackinnonetal2023joe.
The remainder of this paper is organized as follows. Section (ref) reviews the approximate randomization test proposed by CRS and outlines two motivating examples for combining clusters. Section (ref) analyzes the local asymptotic power of CRS in order to develop the necessary tools for combining clusters. Section (ref) presents the data-driven procedures to combine clusters. Sections (ref) and (ref) present results from Monte Carlo simulations and an empirical application, respectively. Section (ref) concludes. All proofs can be found in the appendix.
Consider the following linear regression model of clustered data. Clusters are indexed by $j \in \cJ \equiv \{1, \ldots, q\}$ and units in the $j$th cluster are indexed by $i \in \cI_{n,j} \equiv \{1, \ldots, n_j\}$:
where $Y_{i,j}$ is an outcome variable, $X_{i,j} \in \bR^{d_x}$ is a vector of covariates, $U_{i,j} \in \bR$ is a residual, and $\beta \in \bR^{d_x}$ is an unknown parameter. The total number of observations is $n \equiv \sum_{j \in \cJ} n_j$.
The goal is to test the hypothesis:
for a given vector $c \in \bR^{d_x} \backslash \{0_{d_x}\}$ (where $0_{d_x}$ is a vector of $d_x$ zeros) and $\lambda \in \bR$ at level $\alpha \in (0, 1)$.
This section briefly reviews the approximate randomization test by CRS in the context of the model in Section (ref). The following five-step procedure is taken from caietal2023jem.
CRS requires the clusters to satisfy a joint convergence in distribution requirement and the corresponding limiting random variables to be invariant to sign changes. With these weak assumptions, the test is valid for a fixed number of clusters as $n \longrightarrow \infty$ under the null hypothesis. See Section 3 in canayetal2017ecta and Section 4 in caietal2023jem for more details. The assumptions are stated as follows.
Assumption (ref).(ref) requires that none of the clusters have a negligible size. Assumption (ref).(ref) requires that $\widehat\beta_{n,j}$ obtained from cluster-by-cluster regression is asymptotically normal. It is satisfied when there is weak dependence within clusters so that an appropriate central limit theorem applies. Note that Assumption 3.1 of canayetal2017ecta imposes a weaker version in that it requires the limiting random variable to be symmetric, not necessarily normal. Assumption (ref).(ref) is a standard assumption that requires the clusters to be independent of each other.
Assumption (ref) with $\cJ$ being the set of clusters is appropriate when the target parameter can be identified within each cluster. It is not immediately satisfied when there is an identification problem within clusters. However, suppose the clusters can be combined into a coarser level $\omega$ so that the target parameter is identified in each combined cluster. Then, Assumption (ref) holds when $\cJ$ is replaced by $\omega$. In addition, CRS controls size when $\cJ$ is replaced by $\omega$.
Recall that Step 1 of Algorithm (ref) requires the researcher to obtain $\widehat \beta_{n,j}$ using cluster-by-cluster regressions. This is not possible when the target parameter cannot be identified within each cluster. This section outlines two common empirical settings with this identification issue.
This section analyzes the local asymptotic power of CRS in order to develop the tools for combining the clusters in Section (ref).
Let $\delta \in \bR$ be a local parameter such that
The object of interest is the limiting rejection probability of CRS under the local alternative as described in (ref), i.e.,
where the critical value $\widehat{\text{cv}}_n(1-\alpha)$ is the $(1-\alpha)$-quantile of $\{T_n(g) : g \in \bG\}$ as defined in Step 4 of Algorithm (ref), and $\bP_\delta$ is used to emphasize the dependence of the distribution of the data on the local parameter $\delta$.
As mentioned in the introduction, the asymptotic framework here is $n \longrightarrow \infty$, while holding the number of clusters $q$ fixed. With $q$ being fixed, $|\bG|$ is also fixed. Hence, the “randomization distribution” may not settle down even when $n \longrightarrow \infty$. Therefore, results on the large-sample properties of the local power for randomization tests such as hoeffding1952ams cannot be applied here. See Remark 3.5 of canayetal2017ecta for more discussion.
In order to compute the local asymptotic power of CRS, note that there can be at most $2^{q-1}$ unique values in $\{T_n(g) : g \in \bG\}$. This follows by the definition of $T_n(g)$ in equation (ref) that
for any $g \in \bG$. Therefore, I focus on the set of sign changes $\bG_U \equiv \{g \in \bG : g_1 = 1\}$ that give the unique values of $T_n(g)$ to avoid duplications as in (ref). By construction, $|\bG_U| = 2^{q-1}$. The critical value $\widehat{\text{cv}}_n(1-\alpha)$ can also be evaluated as the $(1 - \alpha)$-quantile of $\{T_n(g) : g \in \bG_U\}$ instead of $\{T_n(g) : g \in \bG\}$.
For notational convenience, define $\bG_{U, -1} \equiv \bG_U \backslash \{1_q\}$, where $1_q \equiv (1, \ldots, 1)$ is the identity transformation. In addition, let $L \equiv |\bG_{U, -1}| = 2^{q-1} - 1$ and $\{\overline{g}_1, \ldots, \overline{g}_L\}$ be the set of all distinct sign changes in $\bG_{U, -1}$.
Following the discussion in the last section, the event $\{T_n > \widehat{\mathrm{cv}}_n(1-\alpha)\}$ can be interpreted in terms of the sign changes in $\bG_U$. In particular, this event collects all arrangements of $\{T_n(g) : g \in \bG_U\}$ such that $T_n$ is one of the largest $K \equiv \lfloor \alpha |\bG_U|\rfloor$ terms. But recall in (ref) that each $T_n(g)$ is defined as first summing the random variables $\widehat{S}_{n,j}$ over $j \in \cJ$ multiplied by the sign changes, and then taking the absolute value. Therefore, computing $\{T_n > \widehat{\mathrm{cv}}_n(1-\alpha)\}$ is a non-trivial task because all the terms in $\{T_n(g) : g \in \bG_U\} $ are dependent on each other.
Observe that computing the probability that $T_n$ is the $k$th largest term in $\{T_n(g) : g \in \bG_U\}$ amounts to comparing $T_n$ with $T_n(g)$ for each $g \in \bG_{U, -1}$. In particular, each comparison involves checking whether $T_n > T_n(g)$ or not for all $g \in \bG_{U, -1}$. The following lemma shows that these comparisons can be expressed explicitly in terms of the partial sums of $\widehat{S}_{n,j}$ depending on the sign change $g \in \bG_{U, -1}$.
The above lemma shows that each comparison of $T_n(h) > T_n(g)$ can be written as two disjoint events. The lemma helps to simplify the computation of local asymptotic power.
Next, recall from the beginning of this subsection that computing the local asymptotic power requires computing the probability that $T_n$ is the $k$th largest term in $\{T_n(g) : g \in \bG_U\}$ for each $k = 1, \ldots, \lfloor \alpha |\bG_U| \rfloor$. The next lemma considers the probability that $T_n$ is less than a subset of terms from $\{T_n(g) : g \in \bG_{U, -1}\}$ and is larger than the remaining terms. It uses Lemma (ref) to represent this as a comparison using partial sums.
In the above lemma, the partial sums $V_{\text{same}}(g, \delta)$ and $V_{\text{diff}}(g, \delta)$ are based on the partial sums in (ref). They are similar to those in Lemma (ref) but taking $n \longrightarrow \infty$, and are evaluated under the local alternative as in (ref). Since the comparisons are about $T_n$ and $T_n(g)$ for $g \in \bG_{U, -1}$, the lemma sets $h = 1_q$ directly from Lemma (ref). They describe the event that $T_n$ is larger than or smaller than another $T_n(\overline{g}_l)$ for $l = 1, \ldots, L$.
Next, the following is a mild assumption used to control the ties of $\{T_n(g) : g \in \bG_U\}$.
See Lemmas S.5.1 to S.5.3 in canayetal2017ecta-supp on how Assumption (ref) is related to various commonly-used test statistics. As they have stated in their lemmas, the assumption is satisfied when the distribution of $Z_j$ is absolutely continuous with respect to the Lebesgue measure for each $j \in \cJ$.
The following theorem states the main result that expresses the local asymptotic power using the above two lemmas. The main idea of the theorem is to consider different ordering of the terms in $\{T_n(g) : g \in \bG_U\}$ and keep those such that $T_n$ is one of the $\lfloor \alpha |\bG_U|\rfloor$ largest terms.
In words, equation (ref) considers the ordering of the terms in $\{T_n(g) : g \in \bG_U\}$ and keeps those such that $T_n$ is one of the $\lfloor \alpha |\bG_U|\rfloor$ largest terms. Hence, there are three summations in the local asymptotic power expression. The first summation (over $k$) sums over the events that $T_n$ is the $k$th largest term. For each $k$, the second summation (over $\cH$) sums over different combinations of the events where $T_n(g) > T_n$ for all $g \in \cH \subseteq \bG_{U, -1}$ with $|\cH| = k-1$. The third summation (over $m$) considers different combinations of the events $\{\cF^{(m_l)}(\overline{g}_l, \cH, \delta)\}_{l =1}^L$ as in Lemma (ref). The expression in Theorem (ref) involves computing the joint probability of multivariate normal distributions. Closed-form expressions may not be always available, but the expression can be computed numerically.
When $K = 1$, the expression in (ref) simplifies greatly. The reason is that with $K = 1$, $\bH_0 = \emptyset$, and the corresponding event $\{T_n > \widehat{\mathrm{cv}}_n(1 - \alpha)\}$ means that $T_n > T_n(g)$ for all $g \in \bG_{U, -1}$. The following corollary shows that when $K = 1$, the local asymptotic power in Theorem (ref) simplifies to a sum of two products of univariate normal cumulative distribution functions. Since $|\bG_{U}| = 2^{q-1}$, the range of significance levels $\alpha$ such that $\lfloor \alpha |\bG_U|\rfloor = 1$ is $\alpha \in [ \frac{1}{2^{q-1}} , \frac{1}{2^{q-2}})$.
A few useful properties of the local asymptotic power for $K=1$ above are summarized below.
Property (ref).(ref) shows that $\pi_{\mathrm{L}}(\delta)$ and $\pi_{\mathrm{R}}(\delta)$ move in opposite directions in $\delta$. Property (ref).(ref) shows that $\pi_{\mathrm{L}}(\delta)$ and $\pi_{\mathrm{R}}(\delta)$ intersect once at $\delta = 0$. Property (ref).(ref) is implied by the first two properties. It states that when $\delta \neq 0$, then one of $\pi_{\mathrm{L}}(\delta)$ or $\pi_{\mathrm{R}}(\delta)$ is larger than the other. The smaller term is bounded above by $\frac{1}{2^q}$, which is usually a small value. Finally, the functions $\pi_{\text{L}}(\delta)$ and $\pi_{\text{R}}(\delta)$ can be interpreted as the local asymptotic power of one-sided tests.
This section develops data-driven procedures to combine clusters based on the theory developed in Section (ref). The procedures aim to find the best combination of clusters to maximize the local asymptotic power.
Assume the clusters described by $\cJ$ can be classified into two groups such that the identification issue can be solved by combining one cluster from each group. This is formalized in the assumption below.
In the case where treatments are at the cluster level, as in Examples (ref) and (ref), Assumption (ref) is automatically satisfied with $\cJ_0$ being the set of control clusters and $\cJ_1$ being the set of treated clusters. For exposition purposes, I also assume that the numbers of clusters in $\cJ_0$ and $\cJ_1$ are equal. Hence, the goal of this section is to pair the clusters. This assumption is relaxed in Section (ref).
Under Assumption {(ref).(ref)}, the size of $\Omega$ is $\overline q!$. For each method of combining clusters $\omega \in \Omega$, write $\omega \equiv \{ \{1, \omega_1\}, \{2, \omega_2\}, \ldots, \{\overline q, \omega_{\overline q}\}\} $, where each pair $\{j, \omega_j\}$ represents the clusters $j \in \cJ_0$ and $\omega_j \in \cJ_1$ being grouped together. The following example illustrates the notations.
As mentioned in Section (ref), under the scenario where clusters have to be combined, the weak assumptions required on the clusters remain the same by replacing $\cJ$ in Assumption (ref) with $\omega \in \Omega$. To emphasize the dependence on the pair of clusters $\{j, \omega_j\}$ from $\omega \in \Omega$ and that the pair of clusters are being considered instead of the individual clusters, the notations $n_j$, $\xi_j$, $\widehat\beta_{n,j}$, and $\sigma_{j}$ that are indexed by $j$ in the last two sections are updated to $n_{j, \omega_j}$, $\xi_{j,\omega_j}$, $\widehat\beta_{n,j,\omega_j}$, and $\sigma_{j,\omega_j}$ respectively.
With the above notations, the goal is to solve the following problem in a computationally efficient way:
where $\pi(\omega, \delta, \alpha)$ is the local asymptotic power at significance level $\alpha \in (0, 1)$ with local alternative parameter $\delta \in \bR$ from (ref) for a given method of combining clusters $\omega \in \Omega$. A simple approach to solve the optimization problem (ref) would be to compute $\pi(\omega, \delta, \alpha)$ for each $\omega \in \Omega$. However, enumerating all possible groupings of clusters and computing $\pi(\omega, \delta, \alpha)$ for each $\omega \in \Omega$ can be computationally intensive when $\overline q$ is large because there are $\overline q!$ different ways in the case of pairing clusters. This motivates the development of computationally efficient procedures for solving (ref).
Let $\widehat\pi_n(\omega, \delta, \alpha)$ be the estimator of the objective function in (ref). In addition, let $\widehat\omega_n$ be the solution to the sample analog of problem (ref), i.e., using $\widehat\pi_n(\omega, \delta, \alpha)$ as the objective in the optimization problem. Note that $\Omega$ remains unchanged when solving the sample analog of (ref). To evaluate $\widehat\pi_n(\omega, \delta, \alpha)$, it requires estimating parameters and variances. This means one has to estimate $\xi_{j,r}$ and $\sigma_{j,r}$ for all $j \in \cJ_1$ and $r \in \cJ_0$. The term $\xi_{j,r}$ is the ratio of cluster size. For $\sigma_{j,r}$, it requires the researcher to specify a working model on the dependence structure and apply a variance estimator. For instance, variance estimators are available for time-dependent data (such as neweywest1987ecta) or spatial data (such as conley1999joe). Alternatively, one can also estimate variance using quasi-maximum likelihood estimation as in caoetal2022wp. Note that even though $\sigma_{j,r}$ has to be estimated, these variance estimators are not used in Algorithm (ref). Thus, specifying the working model and the variance estimator does not mean conducting inference using the variance estimator directly.
In establishing the validity of the data-driven procedure below, Theorem (ref) ahead does not require the correct specification of the working model. However, the choice of the working model can affect the power performance as it affects the objective function of (ref).
To justify the validity of the test that combines the clusters according to the optimal solution to the sample analog of (ref), consider the following assumption that requires how the clusters are combined has a negligible impact on the test. This is similar to the condition imposed in caoetal2022wp.
The above assumption can be satisfied in different settings. For example, this holds through sample-splitting where part of the data is used to obtain $\widehat\omega_n$, and another part is used to conduct inference. Note that Assumption (ref) is also satisfied under some weak assumptions as follows.
Assumption (ref).(ref) is a mild condition requiring the estimator to converge for each $\omega \in \Omega$. As before, this assumption does not require correctly specifying the local asymptotic power function. Assumptions (ref).(ref) and (ref).(ref) require that $\omega^\star \in \Omega$ is the unique solution to (ref). The following proposition shows that Assumption (ref) is satisfied under these conditions.
The following theorem establishes the validity of the test that chooses $\widehat\omega_n \in \Omega$ by solving problem (ref).
The above theorem shows that the data-driven procedure controls size under some high-level conditions when the sample analog of problem (ref) is used. The validity of the procedure does not require the local asymptotic power function to be correctly specified, although having an incorrect specification can lead to power loss.
This section develops a procedure to combine clusters for $K = 1$, i.e., when the significance level is chosen such that $\alpha \in [\frac{1}{2^{\overline q-1}}, \frac{1}{2^{\overline q - 2}} )$. The procedure utilizes the local asymptotic power derived in Corollary (ref). Using the notations from Section (ref), the local asymptotic power for $K = 1$ with a given $\omega \equiv \{ \{1, \omega_1\}, \{2, \omega_2\}, \ldots, \{\overline q, \omega_{\overline q}\}\} \in \Omega$ and $\delta \in \bR$ can be written as
for $\alpha \in [\frac{1}{2^{\overline q-1}}, \frac{1}{2^{\overline q - 2}} )$.
To develop a computationally feasible approach for solving optimization problem (ref) with objective (ref), recall from Assumption (ref) that each cluster in $j \in \cJ_0$ is matched with a cluster in $r \in \cJ_1$. Moreover, each cluster from $\cJ_0$ and $\cJ_1$ is matched only once. As a result, there are only $\overline q$ possible terms $\sigma_{j, r}$ and $\xi_{j,r}$ for each $j \in \cJ_0$. Hence, there are $\overline q^2$ possible terms to compute in total. This is less than $\overline q!$ when $\overline q \geq 4$. See the example below.
Let $\Psi_{j, r} \equiv \Phi \left( -\frac{\xi_{j, r} \delta}{\sigma_{j, r}} \right)$ for each $j \in \cJ_0$ and $r \in \cJ_1$. There are $\overline q^2$ possible combinations of this term by the above reasoning. Following the above discussion, these terms can also be precomputed in advance. Let $z_{j, r} \in \{0, 1\}$ be an indicator variable that equals 1 if clusters $j \in \cJ_0$ and $r \in \cJ_1$ are paired together, and equals 0 otherwise. Then, the optimization problem (ref) with objective (ref) can be formulated as the following integer program:
While the optimization problem (ref) has transformed the optimization problem (ref) with objective (ref) into an integer program with linear constraints, the objective function is nonlinear in $z_{j,r}$. To obtain a computationally feasible program, recall from Property (ref).(ref) that if $\delta < 0$, then
where the relationship is reversed for $\delta > 0$. Since $z_{j, r}$ is a binary variable, taking logarithm on inequality (ref) gives
Note that $\log \Psi_{j,r}$ and $\log (1 - \Psi_{j,r})$ have been pre-computed in advance. Hence, the terms on the left-hand side and the right-hand side of (ref) are linear in the unknowns $z_{j,r}$. Next, inequality (ref) implies that the RHS is in the small interval $[0, \frac{1}{2^{\overline q}}]$. The interval $[0, \frac{1}{2^{\overline q}}]$ can be partitioned into $A$ subintervals $[\epsilon_0, \epsilon_1], [\epsilon_1, \epsilon_2], \ldots, [\epsilon_{A-1}, \epsilon_A]$ where $0 < \epsilon_0 < \epsilon_1 < \cdots < \epsilon_A = \frac{1}{2^{\overline{q}}}$. Since the log of RHS of (ref) is considered in (ref), the number $\epsilon_0$ can be chosen to be $\min_{j, r \in \{1, \ldots, \overline q\}} (1 - \Psi_{j,r})^{\overline q}$ or a very small number above 0 in order to have a finite value for $\log \epsilon_0$. If $\delta < 0$, the optimization problem (ref) can be approximated as solving multiple integer linear programs where each program maximizes the left-hand side of (ref) and restricts the right-hand side of (ref) in a small interval $[\log \epsilon_{a-1}, \log \epsilon_a]$ for each $a = 1, \ldots, A$:
Since (ref) is merely a linear program with binary variables, solving $A$ of them is fast even for large values of $A$. After solving the optimization problem (ref) for each $a = 1, \ldots, A$, the solution can be chosen to be the one that gives the largest local asymptotic power. In practice, (ref) is replaced with its sample analog, i.e., by replacing $\Psi_{j,r}$ with $\widehat{\Psi}_{j,r}$ for each $j \in \cJ_1$ and $r \in \cJ_0$. This is done by estimating $\xi_{j,r}$ and $\sigma_{j,r}$ as discussed in Section (ref).
If $\delta > 0$, the procedure is the same except the 0-1 linear integer program in (ref) becomes the following program:
The procedure described in this section is summarized as follows.
Algorithm (ref) requires the researcher to choose a value of $A$ in Step 2. The following proposition establishes that when $A$ is large enough, then the approximation in Algorithm (ref) yields the same solution as in the original optimization problem (ref). Hence, by partitioning $[0, \frac{1}{2^{\overline q}}]$ into fine enough intervals, the solution from the approximation step coincides with the solution from (ref).
The above proposition assumes the researcher partitions the interval $[0, \frac{1}{2^{\overline q}}]$ into equally-spaced intervals. This assumption is for convenience purposes only as it seemed to be a natural way to run Algorithm (ref).
This section ends with three remarks on the above data-driven procedure for combining clusters.
The previous section focuses on developing an algorithm that uses the local asymptotic power of $K = 1$ as the objective function. In practice, researchers may be interested in conducting a test that involves $K > 1$. The reason for having a simple 0-1 integer linear program in the previous section is because the local asymptotic power for $K = 1$ has a simple expression as stated in Corollary (ref). The other cases require evaluating the expression in Theorem (ref) that involves more complicated expressions but can be evaluated using numerical methods.
Similar to the discussion in Section (ref) on the estimation of variances, evaluating the rejection probability in Theorem (ref) requires the researcher to specify a working model on the dependence structure and apply a variance estimator. Recall that Theorem (ref) does not impose any requirement on $K$. Hence, the choice of the working model and the variance estimator does not affect the validity of the procedure, but can affect the power performance.
Algorithm (ref) is found to be effective in computing the optimal grouping of clusters for $K > 1$ in simulations. Hence, it can be used to obtain an initial solution. When $\overline q$ is large and $K>1$, evaluating the local asymptotic power for each $\omega \in \Omega$ can be costly. The following algorithm provides a heuristic to proceed from the initial solution. The main idea is to start from the initial point provided by Algorithm (ref). Then, one proceeds to see if it is possible to improve upon the existing solution using local search. This algorithm below is motivated by the 2-opt algorithm for the traveling salesman problem (see, for instance, flood1956or, croes1958or, and linkernighan1973or).
Algorithm (ref) presents a procedure by conducting pairwise swaps. This can be easily extended for swapping more groupings.
This section briefly discusses the case with an unequal number of clusters in $\cJ_0$ and $\cJ_1$ for Algorithm (ref). Assume without loss of generality that $|\cJ_0| < |\cJ_1|$.
Following the notations in Section (ref), let $\Omega$ be the collection of all methods of combining clusters. In addition, let $\cJ_0 = \{1, 2, \ldots, \overline{q}\}$. Hence, the goal is to form $\overline{q}$ groups of clusters, where each group has one cluster from $\cJ_0$, and at least one cluster from $\cJ_1$. Moreover, all clusters from $\cJ_1$ have to be used. As before, each cluster can only be used once in the grouping. Therefore, each $\omega \in \Omega$ can be written as $\{\{1, \omega_1\}, \{2, \omega_2\}, \ldots, \{\overline{q}, \omega_{\overline{q}}\}\}$ as before, with the exception that $\omega_j$ can represent more than one cluster from $\cJ_1$ here for each $j = 1, \ldots, \overline{q}$. Let $\cM$ be the set of all possible ways of partitioning clusters in $\cJ_1$ into $\overline{q}$ groups, so that $\omega_j \in \cM$ for each $j = 1, \ldots, \overline{q}$, and that each $m \in \cM$ is nonempty and satisfies $m \subset \cJ_1$. For example, if $\cJ_0 = \{1, 2\}$ and $\cJ_1 = \{3, 4, 5\}$ and one wishes to form two groups of clusters, the set $\cM$ can be written as $\cM = \{\{3\}, \{4\}, \{5\}, \{3, 4\}, \{3, 5\}, \{4, 5\}\}$.
Similar to Section (ref), let $z_{j,m} \in \{0, 1\}$ be a binary variable that equals 1 when cluster $j \in \cJ_0$ is combined with $m \in \cM$. Define $\Psi_{j,m} \equiv \Psi \left( -\frac{\xi_{j,m}\delta}{\sigma_{j,m}}\right)$, where $\xi_{j,m}$ is the ratio of cluster size and $\sigma^2_{j,m}$ is the asymptotic variance in the group of cluster $\{j, m\}$, where $j \in \cJ_0$ and $m \in \cM$.
For each $a = 1, \ldots, A$ and the corresponding interval $[\log\epsilon_{a-1}, \log\epsilon_a]$, the optimization problem analogous to problem (ref) for $\delta < 0$ in the scenario of $|\cJ_0| < |\cJ_1|$ is as follows:
Similarly, if $\delta > 0$, the optimization problem becomes
The above optimization problems are binary integer programs as Section (ref), and are fast to solve.
This section explores the finite-sample performance of the data-driven procedures to combine clusters in Section (ref). Here, I consider the data generating process (DGP) based on the simulation design in hagemann2022wp, which is a version of the simulation design in conleytaber2011restat.
Let there be $q = 12$ clusters, which consists of 6 treated clusters represented by $\cJ_1 = \{1, 2, \ldots, 6\}$ and 6 control clusters represented by $\cJ_0 = \{7, 8, \ldots, 12\}$. In the following, $j \in \cJ \equiv \cJ_1 \cup \cJ_0$ indices all the clusters and $t \in \cT \equiv \{1, 2, \ldots, T\}$ indices time. Let $D_{t,j} \equiv \ind[j \in \cJ_1 \text{ and } t > t_0]$ for $j \in \cJ$ and $t \in \cT$ and $I_t \equiv \ind[t > t_0]$ for $t \in \cT$. Each simulated data is generated from the following model:
with $\theta_0 = \gamma_1 = \gamma_2 = \gamma_3 = \xi_j = 1$, $\rho = 0.5$, $\gamma_4 = 0.8$, $t_0 = 10$, and $T = 20$ as in hagemann2022wp. I consider three different DGPs below. DGP 1 generates the random variables in the same way as hagemann2022wp. DGPs 2 and 3 have different variance structures and assume there is the same number of heterogenous treated and control clusters. In all the specifications, the parameter $h \in \{1, 2, 3, 4\}$ governs the amount of heterogeneity in data, where a larger value of $h$ represents a higher degree of heterogeneity. The goal is to test the two-sided hypothesis $H_0 : \beta = 0$ vs $H_1 : \beta \neq 0$. I consider two significance levels $\alpha = 0.05$ and $\alpha = 0.1$ in order to study the finite-sample properties of Algorithms (ref) and (ref) separately. The specification for each $j \in \cJ$ are as follows.
For each simulation design, I estimate the rejection rate obtained from the procedure in Section (ref) and use CRS to conduct inference for each of the $6!=720$ different methods of combining clusters. Apart from CRS, I also apply other recent methods for conducting inference with a fixed number of clusters that were mentioned in Section (ref). The methods and their abbreviations are summarized in Table (ref). Some details of considering the two versions of CRS in Table (ref) are as follows. First, CRS-Data is the data-driven procedure that uses optimization problem (ref) for $\delta < 0$ and problem (ref) for $\delta > 0$ when $\alpha = 0.05$. For these two optimization problems, I partition the interval $[0, \frac{1}{2^5}]$ into $A = 200$ subintervals. When $\alpha = 0.1$, the heuristic in Algorithm (ref) is used. One of the two programs (ref) and (ref) is used to obtain an initial solution, depending on the sign of $\delta$. Next, CRS-Random chooses one of the 720 ways of combining clusters randomly in each draw. This aims to represent the scenario where the researcher chooses a random grouping of clusters that can solve the identification problem for each data set. For each specification, I use 20,000 Monte Carlo simulations. The data-driven procedure CRS-Data uses the local alternative parameter $\delta = 2\sqrt{q T}$ under $\beta \geq 0$ and $\delta = -2\sqrt{q T}$ under $\beta < 0$. I consider $\beta \in \{-3, -2, \ldots, 3\}$ in the simulation exercises.
Figure (ref) presents the rejection rate curves for various DGPs, heterogeneity parameters $h$, and methods for conducting inference at significance level $\alpha=0.05$. Each column corresponds to one of the three DGPs, with different amounts of heterogeneity. In each figure, the grey region is formed by overlapping 720 different rejection rate curves based on each of the different ways of combining clusters, and conducting inference using CRS. The other colored lines correspond to various methods as indicated in the legend with the abbreviations specified in Table (ref).
For DGP 1 (column 1 of Figure (ref)), the treated clusters are homogeneous for the different values of $h$ because the heterogeneity in variance only affects the control clusters. Hence, how the clusters are combined should not affect the population local asymptotic power as discussed in Remark (ref). Indeed, both CRS-Data and CRS-Random give similar performance in finite-samples. On the other hand, BCH and WCB over-rejects under the null when $h = 4$.
For the other DGPs (columns 2 and 3 of Figure (ref)), there is heterogeneity in both the treated and control clusters. Hence, there is more variation in the CRS rejection rates depending on how the clusters are combined. CRS controls size in all designs. In addition, CRS-Data (as indicated by the solid black lines) always belong to the top regions formed by the grey curves. This shows that Algorithm (ref) can lead to close to oracle performance in the simulations. Besides, CRS-Data also performs better than CRS-Random (as indicated by the solid blue lines). This shows that the data-driven procedures using (ref) or (ref) perform better than choosing a random grouping of clusters. On the other hand, WCB can over-reject in designs with more heterogeneity.
Next, I examine the performance of Algorithm (ref) by considering $\alpha = 0.1$ for the same set of simulations. The results are summarized in Figure (ref). The observations for $\alpha = 0.1$ are similar to those when $\alpha = 0.05$ as in Figure (ref). CRS controls size in all designs. CRS-Data continues to be able to lead to a close-to-oracle grouping of clusters as the solid black lines belong to the top parts of the grey regions. In addition, CRS-Data compares favorably to other methods in various DGPs and heterogeneity parameters.
To conclude, the simulation results in this section show that the data-driven procedures to combine clusters perform well under various settings.
dinceccokatz2016ej study the impact of fiscal and administrative powers of states on economic performance, using data from 11 European countries. Table 3 of dinceccokatz2016ej regresses real per capita GDP growth on two binary variables representing fiscal centralization and limited government, together with other variables and fixed effects. Their results are clustered at the country level.
The baseline regression they consider is:
where $j$ indices country, $t$ indices year, $Y_{t,j}$ is the (logarithm) annual growth rate of real per capita GDP in country $j$ between years $t-1$ and $t$, $C_{t,j} \in \{0, 1\}$ equals 1 for fiscal centralization in country $j$ and year $t$, $L_{t,j} \in \{0, 1\}$ equals 1 for limited government in country $j$ and year $t$, and $\mu_j$ is the country fixed-effects. The goal is to conduct inference on $\widehat \beta_1$ and $\widehat \beta_2$.
I focus on columns (1) to (3) of Table 3 in dinceccokatz2016ej in order to focus on the issue that there exist clusters in which the binary variables $C_{t,j}$ or $L_{t,j}$ do not have variation over time. The specifications for the three columns are as follows.
In this situation of having only 11 clusters, CRS can be used to conduct valid inference. However, there are clusters with no variation in the main variables of interest. There are three countries with $C_{t,j} = 1$ in all periods (Belgium, England, and Piedmont). There are two countries with $L_{t,j} = 1$ in all periods (Belgium and Denmark). As a result, researchers would have to combine clusters in order to ensure that all parameters can be identified in each combined cluster. In addition, the data is at the aggregate level in that there is one observation for each year in each cluster. Therefore, researchers would also have to combine clusters in the presence of time fixed effects. In this section, I form five groups of clusters across all specifications.
In order to study the performance of the data-driven procedures in Section (ref), this section performs a calibrated simulation exercise. I first estimate the model using the actual data. Then, I simulate data by using various distributions of the error terms. In each specification, let $\widehat\beta_{1}$ and $\widehat\beta_{2}$ be the estimates of $\beta_1$ and $\beta_2$ obtained from data respectively. I take them as the true values of $\beta_1$ and $\beta_2$ in the calibrated simulation procedure. In addition, I perform inference on the coefficients of the two binary variables fiscal centralization $C_{t,j}$ and limited government $L_{t,j}$ separately. Hence, in conducting calibrated simulation exercises for the coefficient on $C_{t,j}$, I set $\beta_1 = \widehat\beta_{1}$ under the null, $\beta_1 = \widehat\beta_{1} + \Delta$ under the alternative, and keep $\beta_2 = \widehat\beta_{2}$. Similarly, for the coefficient on $L_{t,j}$, I set $\beta_2 = \widehat\beta_{2}$ under the null, $\beta_2 = \widehat\beta_{2} + \Delta$ under the alternative, and keep $\beta_1 = \widehat\beta_{1}$. The significance level is $\alpha=0.1$ for all simulations. There are 9,660 different ways to combine the clusters in this calibrated simulation exercise. I use 2,000 Monte Carlo simulations for each column and specification. I consider $\Delta \in \{-3, -2, \ldots, 3\}$.
In the calibrated simulation exercise, I first use the residuals $\{U_{t,j}\}$ to fit an AR(1) model \[ U_{t,j} = \rho_j U_{t,j-1} + \epsilon_{t,j}, \] where $\epsilon_{t,j} \sim N(0, \nu_j^2)$ for each country $j \in \cJ$. Using the estimated parameters $\{\widehat\rho_j\}_{j \in \cJ}$ and $\{\widehat\nu_j\}_{j \in \cJ}$, I consider the following specifications of the error terms in conducting the calibrated simulation exercise.
More details of the calibrated simulation exercise can be found in Appendix (ref).
Figures (ref) and (ref) report the results for inference on the coefficients of the variables $C_{t,j}$ and $L_{t,j}$ respectively. Each figure reports the rejection rates against $\Delta$, where $\Delta$ was defined two paragraphs ago. $\Delta = 0$ corresponds to the results under the null. Based on the findings in Section (ref) that CRS has favorable performance across various designs, I focus on examining the performance of the data-driven approach here. See Table (ref) for the definitions of CRS-Data and CRS-Random. Similar to Section (ref), the grey regions are formed by the 9,660 rejection rate curves that are based on the different ways of forming five groups out of 11 clusters. As in the previous section, the data-driven procedure CRS-Data uses the local alternative parameter $\delta = 2\sqrt{q T}$ under $\Delta \geq 0$ and $\delta = -2\sqrt{q T}$ under $\Delta < 0$, where $T$ is the number of periods in the calibrated simulation exercise.
The simulation results show that CRS-Data controls size and chooses powerful grouping of clusters in various designs. It also continues to perform better than choosing a random grouping of clusters.
The approximate randomization test in canayetal2017ecta imposes weak assumptions on clusters and can be used to conduct valid inference when there is a small number of clusters. The test requires researchers to perform cluster-by-cluster regressions. When the target parameters cannot be identified within each cluster, researchers would have to combine clusters to perform the test.
In this paper, I first analyzed the local asymptotic power of CRS. Using the local asymptotic power as a criterion, I developed computationally efficient algorithms that can be used to guide how to combine clusters when there is an identification issue within cluster. Monte Carlo simulations and an empirical application show that the data-driven procedure performs well.