Permutation inference with a finite number of heterogeneous clusters

Permutation inference with a finite number of heterogeneous clusters

\address{University of Michigan, Stephen M. Ross School of Business, 701 Tappan Ave, Ann Arbor, MI 48109, USA. Tel.: +1 (734) 615-6663} \email{[email removed]} \urladdr{\href{https://umich.edu/ hagem}{umich.edu/ hagem}}

Introduction

It has become widespread practice in economics to conduct inference that is robust to within-cluster dependence. Typical examples of clusters are states, counties, cities, schools, firms, or stretches of time. Units within the same cluster are likely to influence one another or are influenced by the same external shocks. Several analytical and computationally intensive procedures such as the bootstrap are available to account for the presence of data clusters. Most of these procedures achieve consistency by requiring the number of clusters to go to infinity. Numerical evidence by bertrandetal2004, mackinnonwebb2014, and others suggests that this type of asymptotics often translates into heavily distorted inference in empirically relevant situations when the number of clusters is small or the clusters are heterogenous. In both situations, the overall finding is that true null hypotheses are rejected far too often. In this paper, I introduce an adjusted permutation procedure that is able to asymptotically control the size of tests about the effect of a binary treatment in the presence of finitely many large and heterogeneous clusters. The procedure applies to difference-in-differences estimation and other situations where treatment occurs in some but not all clusters and the treatment effect of interest is identified by between-cluster comparisons.

The main theoretical insight of this paper is that classical permutation inference can be adjusted to test the null hypothesis of equality of means of two finite samples of mutually independent but arbitrarily heterogeneous normal variables. This runs counter to classical permutation testing hoeffding1952, where the data under the null are presumed to be exchangeable. The adjustment corrects the significance level of the test downwards to account for heterogeneity. I prove that this is possible for empirically relevant levels of significance if both samples consist of more than three observations. The corrections needed for all standard levels of significance are tabulated in the paper. I also show that if a random vector of interest converges weakly to multivariate normal with diagonal covariance matrix, then permutation inference remains approximately valid for that vector. To exploit this result in a cluster context, I construct asymptotically normal statistics from each cluster and then apply adjusted permutation inference to the collection of these statistics. The resulting permutation test is consistent against all fixed alternatives to the null, powerful against local alternatives, and is free of user-chosen parameters.

The strategy of using cluster-level estimates as the basis for a test goes back at least to famamacbeth1973, who without formal justification run $t$ tests on regression coefficients obtained from year-by-year cross-sectional regressions. Their approach is generalized and formalized by ibragimovmueller2010, ibragimovmueller2016, who construct $t$ statistics from cluster-level estimates and show that for certain combinations of numbers of clusters and significance levels these statistics can be compared to Student $t$ critical values. The \citetalias{ibragimovmueller2016} test and the adjusted permutation test complement one another because they both rely on finite-sample inference with heterogeneous normal variables but apply to non-nested combinations of numbers of clusters and significance levels. The empirical example in this paper features a practically relevant situation where the \citetalias{ibragimovmueller2016} test does not apply but the adjusted permutation test does. If both tests apply, the Monte Carlo results in this paper indicate that neither test dominates the other in terms of power but the adjusted permutation test has clear advantages if the underlying data are heavy tailed.

Several other papers show that inference with a fixed number of clusters is possible under a variety of conditions: canayetal2014 permute the signs of cluster-level statistics under symmetry assumptions. This approach requires the parameter of interest to be identified within each cluster and clusters therefore have to be paired in an ad-hoc manner for difference-in-differences estimation. This pairing has a substantial impact on the test decision and requires a large number of choices on the part of the researcher. besteretal2014 use standard cluster-robust covariance matrix estimators but adjust critical values under homogeneity assumptions on the clusters. canayetal2018 show that certain cluster-robust versions of the wild bootstrap can be valid under strong homogeneity assumptions with a fixed number of clusters. In sharp contrast, the test developed here does not require pairing clusters or any other decisions on the part of the researcher and applies even if the clusters are arbitrarily heterogeneous.

I will use the following notation: $1\{\cdot\}$ is the indicator function, $\min\{a,b \} = a\wedge b$, and cardinality of a set $A$ is $|A|$. The smallest integer larger than $a$ is $\lceil a \rceil$ and the largest integer smaller than $a$ is $\lfloor a \rfloor$. Limits are as $n\to\infty$ unless noted otherwise.

All proofs can be found in the online appendix.

Permutation inference with heterogenous symmetric variables

In this section I show that classical permutation inference can be adjusted to test for the equality of location of two finite samples of independent symmetric variables with heterogeneous scales. The discussion focuses on heterogeneous normal variables but several of the results apply more generally.

Suppose the random vector $X = (X_1,\dots, X_q)\in\mathbb{R}^q$ has entries $X_k = \mu_1 + \sigma_k Z_k$ for $1\leqslant k\leqslant q_1$ and $X_k = \mu_0 + \sigma_k Z_k$ for $q_1 +1 \leqslant k\leqslant q_1 + q_0 = q$, where the $Z_1,\dots,Z_q$ are iid symmetric variables. The $\sigma_k$ are not known and no estimates are assumed to be available. The number of variables $q$ is taken as fixed throughout this paper. The goal is to construct an $\alpha$-level permutation test of the hypothesis $H_0\colon \mu_1=\mu_0$. This is a two-sample problem with “treatment” sample $X_1,\dots, X_{q_1}$ and “control” sample $X_{q_1+1},\dots, X_q$. The test statistic $T$ considered here is the comparison of means

No standardization is needed.

Let $\mathfrak{S}_q$ be the group of permutations of the set $\{1,\dots, q\}$. For $g\in\mathfrak{S}_q$, denote by $g(k)$ the value the permutation $g$ assigns to $k$ for $1\leqslant k\leqslant q$. The “group action” on $X$ in $\mathfrak{S}_q$ is the relabeling of the indices $gX = (X_{g(1)},\dots, X_{g(q)})$. A permutation test derives its critical values from the permutation statistics $T(gX)$. Because $x\mapsto T(x)$ is invariant to the ordering of the first $q_1$ and last $q_0$ entries of $x$, it suffices to compute the $T(gX)$ for the set of group actions with unique combinations of $g(1),\dots, g(q_1)$ and $g(q_1 + 1),\dots, g(q)$. One way of representing this set is

Denote by $T^{(1)}(X,\mathfrak{G}) \leqslant T^{(2)}(X,\mathfrak{G})\leqslant \cdots \leqslant T^{(|\mathfrak{G}|)}(X,\mathfrak{G})$ the ordered values of $T(gX)$ as $g$ varies over $\mathfrak{G}$ and define critical values

Classical permutation inference operates under the null hypothesis that $X$ has the same distribution as $gX$ for all $g\in\mathfrak{S}_q$. In the present context this would be equivalent to assuming that $\mu_1 = \mu_0$ and that all $\sigma_k$ are identical under the null. An argument due to hoeffding1952 would then show that $T^{\alpha}(X,\mathfrak{G})$ could be used as the critical value for an $\alpha$-level test against the alternative $H_1\colon \mu_1> \mu_0$. If the null hypothesis is weakened to $H_0\colon \mu_1=\mu_0$ without restrictions on $\sigma_k$, a natural question to ask if there exists any order statistic $j\mapsto T^{(j)}(X,\mathfrak{G})$, $\lceil(1-\alpha)|\mathfrak{G}|\rceil\leqslant j< |\mathfrak{G}|$, that can be used as a critical value for an $\alpha$-level test even if the classical permutation hypothesis $X\sim gX$ for all $g\in\mathfrak{S}_q$ fails. As I will discuss now, the answer to this question is affirmative for empirically relevant choices of $\alpha$ if $q_1$ and $q_0$ are larger than $3$.

Because $T(X)\in \{ T(gX) : g\in\mathfrak{G} \}$, it is always true that $T(X) \leqslant T^{(|\mathfrak{G}|)}(X,\mathfrak{G})$. The largest non-trivial critical value from $\{ T(gX) : g\in\mathfrak{G} \}$ is therefore the second largest order statistic $T^{(|\mathfrak{G}|-1)}(X,\mathfrak{G})$. The following theorem shows that the probability that $T(X)$ exceeds $T^{(|\mathfrak{G}|-1)}(X,\mathfrak{G})$ is necessarily small under $H_0\colon \mu_1=\mu_0$. In fact, this probability is so small that $T(X) > T^{(|\mathfrak{G}|-1)}(X,\mathfrak{G})$ is well below any standard choice of $\alpha$ for most values of $q_1$ and $q_0$. By monotonicity, the existence of a $j$ such that ${\mathord P} (T(X) > T^{(j)}(X,\mathfrak{G}))\leqslant \alpha$ is then guaranteed.

A byproduct of the theorem is a bound for the case where the scales $\sigma_1,\dots, \sigma_q$ are replaced by positive random variables independent of $Z_1,\dots,Z_q$. The $X_k$ are then called “scale mixtures” of a symmetric variable $Z$. The following corollary is immediately obtained from Theorem (ref) by conditioning on a given set of random scales.\footnote{A referee points out that szekely2006 studies one-sample Student $t$-tests for similar classes of distributions. szekely2006 does not deal with permutation inference and uses a fundamentally different proof technique but the results are also powers of two.}

Theorem (ref) shows that a test with critical value $T^{(|\mathfrak{G}|-1)}(X,\mathfrak{G})$ has size $1/2^{q_1\wedge q_0} = 0.0625$, $0.0313$, $0.0156$, $0.0078$, $0.0039$ as $q_1\wedge q_0$ increases from $4$ to $8$. Consequently, a 10%-level permutation test that relies only on symmetry is available with $q_1$ and $q_0$ as small as $4$. One can perform a 5%-level test with $q_1\wedge q_0\geqslant 5$, a 5%-level two-sided test (see the discussion below (ref) ahead) with $q_1\wedge q_0\geqslant 6$, a 1%-level test with $q_1\wedge q_0 \geqslant 7$, and a 1%-level two-sided test with $q_1\wedge q_0 \geqslant 8$.

More generally, Theorem (ref) implies that for many combinations of $q_1$, $q_0$, and $\alpha$ there exist $p\in (0,1)$ such that $\lceil (1-\alpha)|\mathfrak{G}|\rceil \leqslant \lceil (1-p)|\mathfrak{G}|\rceil < |\mathfrak{G}|$ and ${\mathord P} (T(X) > T^{p}(X,\mathfrak{G})) \leqslant \alpha$. The largest such value of $p$ maximizes power while still controlling the size of the test. Finding this $p$ is theoretically and computationally challenging. However, computation can be simplified if $Z$ is restricted to a single distribution. For normal distributions, the best possible $p$ is

where I suppress the dependence on $q_1$ and $q_0$ to prevent notational clutter. By construction, $\bar{\alpha}$ controls the size of the permutation test not only for arbitrarily heterogeneous normal variables but also for the entire class of scale mixtures of normals. This class includes all Student $t$ and Laplace distributions, as well as many other standard distributions gneiting1997. Moreover, because the critical value is from a permutation distribution, the test also controls size for all exchangeable distributions. The remainder of the paper therefore focuses on this $\bar{\alpha}$ and heterogeneous normal $X$ but other choices of distributions are possible.

A convenient feature of $\bar{\alpha}$ is that it does not depend on the data and can therefore be tabulated. To this end, I use a location-scale invariance argument to reduce the inner supremum in (ref) to a supremum over $(0,1]^q$, simulate ${\mathord P}$ over large random grids on $(0,1]^q$, and compute $\bar{\alpha}$ by iteratively searching over these grids. (See Online Appendix (ref) for details.) The search is not exhaustive and does not guarantee that the target quantity in (ref) is found. However, in experiments this method consistently replicated the theoretical result in Theorem (ref) up to a small approximation error, which indicates---but does not unequivocally establish---that this approximation of $\bar{\alpha}$ is reliable.

Table (ref) lists $\bar{\alpha}$ for common choices of $\alpha$ as a function of $q_1$ and $q_0$. As can be seen, the adjustment needed to make inference robust to variance heterogeneity is substantial if $q_1\wedge q_0$ is very small but disappears quickly as $q_1\wedge q_0$ increases. For example, for $q_1 = 4 = q_0$ a robust 10%-level test requires using the 95.62% quantile of the unadjusted test but for $q_1 = 9 = q_0$ the 91% quantile is already sufficient for a robust 10%-level test. For larger numbers of variables the need for adjustment nearly disappears at conventional levels of significance. This is also confirmed by results in hagemann2019, who shows that unadjusted permutation inference in this context with the statistic $T(X)$ is consistent if the number of treated and control units grows in a balanced manner.

The test decision is now simple. For $q_1 \wedge q_0 > 3$, choose $\bar{\alpha}$ for a feasible $\alpha$ from Table (ref) to ensure ${\mathord P} (T(X) > T^{\bar{\alpha}}(X,\mathfrak{G}))\leqslant \alpha$ under $H_0\colon \mu_1 = \mu_0$. The existence of such an $\bar{\alpha}$ for the comparison-of-means test statistic $T$ is guaranteed by Theorem (ref). For an $\alpha$-level test of the null hypothesis $H_0\colon \mu_1 = \mu_0$, reject in favor of the alternative $H_1\colon \mu_1 > \mu_0$ if

For a one-sided test of level $\alpha$ against $\mu_1 < \mu_0$, reject if $T(-X) > T^{\bar{\alpha}}(-X,\mathfrak{G})$ or, equivalently, $T(X) < T^{(\lfloor |\mathfrak{G}|\bar{\alpha} \rfloor)}(X,\mathfrak{G})$. For a two-sided test of level $2\alpha$ against $\mu_1 \neq \mu_0$, reject if $T(X) > T^{\bar{\alpha}}(X,\mathfrak{G})$ or $T(-X) > T^{\bar{\alpha}}(-X,\mathfrak{G})$. Test decisions can also be equivalently made with the p-value of the unadjusted test

because $T(X) > T^{p}(X,\mathfrak{G})$ if and only if $\hat{p}(X, \mathfrak{G}) \leqslant p$ for every $p \in (0,1)$. A $p$-value for a two-sided test can be defined as $2( \hat{p}(X, \mathfrak{G}) \wedge \hat{p}(-X, \mathfrak{G}) ).$ Reject the null hypothesis if the $p$-value does not exceed $\bar{\alpha}$ from Table (ref) to perform an $\alpha$-level test.

Online Appendix (ref) contains additional results on power, stochastic approximation of $\mathfrak{G}$, and large sample approximation of $X$. The next section applies Theorem (ref) to situations where $X$ is the distributional limit of cluster-level statistics.

Permutation inference with heterogenous clusters

In this section, I establish large sample results for an adjusted permutation test with finitely many clusters under a single high-level condition. I then outline how these results can be applied in empirical practice.

\addline

Suppose data from $q$ large clusters (e.g., counties, regions, schools, firms, or stretches of time) are available. Observations are independent across clusters but dependent within clusters. An intervention took place during which clusters $1\leqslant k\leqslant q_1$ received treatment and clusters $q_1 + 1\leqslant k\leqslant q$ did not. The quantity of interest is a treatment effect or an object related to a treatment effect that can be represented by a scalar parameter $\delta$. Because entire clusters receive treatment, this parameter is only identified up to a location shift $\theta_0$ within a treated cluster. Hence, only the left-hand side of \[ \theta_1 = \theta_0 + \delta \] can be identified from such a cluster. If the clusters have similar characteristics, then $\theta_0$ can be identified from an untreated cluster. Comparing the two clusters identifies $\delta$.

The identification strategy outlined in the preceding paragraph is the basis for differences-in-differences estimation---arguably the most popular identification strategy in economics today---and a variety of other models. The purpose of this section is to use the results from Section (ref) to develop a permutation test of the conventional (non-sharp) hypothesis \[H_0\colon \delta = 0, \] or, equivalently, $H_0\colon \theta_1 = \theta_0$. The idea is to obtain independent estimates \smash{$\hat{\theta}_{n,1},\dots, \hat{\theta}_{n,q_1}$} of $\theta_1$ and independent estimates \smash{$\hat{\theta}_{n,q_1 + 1},\dots, \hat{\theta}_{n,q}$} of $\theta_0$ so that $\hat{\theta}_n = (\hat{\theta}_{n,1},\dots, \hat{\theta}_{n,q})$ is approximately multivariate normal with diagonal covariance matrix. The following example outlines a simple situation where this is possible.

The cluster-level statistics $\hat{\theta}_n$ can be combined with the results in the previous section to perform a consistent permutation test as the sample size $n$ grows large. The test is not limited to the $\hat{\theta}_n$ constructed in the preceding example. Instead, the key high-level condition is that a centered and scaled version of some estimate $\hat{\theta}_n$ converges to a $q$-dimensional standard normal distribution,

The $\sigma_1,\dots, \sigma_q$ may depend on $\theta_1$ or $\theta_0$ but are not presumed to be known or estimable by the researcher. This is an important feature of the test because consistent covariance matrix estimation would require knowledge of an explicit ordering of the dependence structure within each cluster. While ordering the data is straightforward for time-dependent data, it may be difficult or impossible to infer or credibly assume an ordering of the data within villages or schools. In contrast, (ref) can be established under weak dependence assumptions where it is only presumed that there exists a possibly unknown ordering for which the dependence decays at a certain rate. machkouriaetal2013 present easy-to-use moment bounds and limit theorems for this situation; see also besteretal2014 for further results.

I now show that under the joint convergence (ref) a permutation test based on comparison of means of $\smash{\hat{\theta}_{n,1},\dots, \hat{\theta}_{n,q_1}}$ and \smash{$\hat{\theta}_{n,q_1 + 1},\dots, \hat{\theta}_{n,q}$} can be adjusted to be asymptotically of level $\alpha$ with a fixed number of clusters. This is possible for $q_1 \wedge q_0 > 3$ if $\bar{\alpha}$ in Table (ref) is available at the desired significance level $\alpha$. In that case, the test has power against fixed alternatives $\theta_1 = \theta_0 + \delta$ with $\delta > 0$ and local alternatives $\theta_1 = \theta_0 + \delta/\sqrt{n}$ converging to the null. In the latter situation, $\theta_0$ is fixed and $\theta_1$ implicitly depends on $n$. The convergence in (ref) is then no longer pointwise in $\theta = (\theta_1,\theta_0)$ but a statement about the sequence $\theta_n = (\theta_0 + \delta/\sqrt{n}, \theta_0)$. As before, the test can be made two-sided to have power against fixed and local alternatives from either direction. Let $x\mapsto\tilde{\Phi}_{\theta_0}(x) = \prod_{1\leqslant k\leqslant q_0}\Phi(x/\sigma_{k+q_1}(\theta_0))$.

Online Appendix (ref) provides more practical guidance for the implementation of the adjusted permutation test and applies the test to several standard econometric models.

Numerical results

This section studies the behavior of the adjusted permutation test and related methods in a Monte Carlo experiment and in data from a randomized trial. The discussion focuses on one-sided tests to the right but the results apply more generally. Online Appendix (ref) contains additional numerical examples and empirical applications.