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.
72,363 characters · 15 sections · 64 citation commands
A Robust Permutation Test for Subvector Inference in Linear Regressions
Inference in linear regressions is one of the oldest problems in statistics. The first tests, developed by Student and Fisher, are still in use today, with their nice features of being exact with independent, normally distributed unobserved terms but also asymptotically valid under homoskedasticity only. Heteroskedasticity is a common phenomenon, however. As a result, many applied researchers nowadays rely instead on robust t- and Wald tests based on robust variance estimators White(1980). These tests are not the panacea, yet: they are only asymptotically valid, and may suffer from important distortions in finite samples MacKinnon(2013), even if the unobserved terms are independent of the regressors.
Recently, DR2017 show that it is actually possible to construct a test sharing the desirable properties of both approaches. Specifically, they develop a permutation test that is exact in finite samples under independence, but also asymptotically valid under weak exogeneity conditions, allowing in particular for heteroskedasticity. However, under general conditions, the exact validity of the test only holds for inference on the whole vector of parameters. Most often, researchers are interested in subvector (e.g., scalar) inference. When applied to such subvectors, their test is exact only if the regressors corresponding to the subvector that is tested and other regressors are independent. This condition is seldom satisfied in practice.\footnote{An important exception is randomized experiments, where the treatment is often drawn independently of all observed variables.} DR2017 also develop a partial correlation test for specific components, analogous to the residual permutation test of Freedman-Lane(1983) (see Toulis(2022) for an extension). Their test is asymptotically valid but is not exact in finite samples, even if the unobserved terms are independent of the regressors.
The objective of this paper is to extend the results of DR2017 by developing a test on subvectors that is exact under a conditional independence assumption but also consistent under a weaker exogeneity condition. To this end, we consider a new, “stratified randomization” test (SR test hereafter) for a parameter vector $\beta$ based on a heteroskedasticity-robust Wald statistic in the partitioned regression
where $X$ is the regressors matrix of interest, $Z$ is the matrix of the nuisance regressors and $u$ is the vector of error terms. The key idea is to stratify the data according to the different values of $Z$ and permute the set of observations within each stratum. It is exact if $X$ is independent of the error term $u$ conditional on $Z$, without any restriction on the dependence between $X$ and $Z$.
The test is also asymptotically valid and has power against local alternatives under a weak exogeneity condition. Specifically, we assume that $X$ and $u$ are uncorrelated conditonal on $Z$. This condition is stronger than the usual condition of no correlation between $u$ and $(X,Z)$, but weaker than the mean independence condition $\operatorname{E}[u|X,Z]=0$. To obtain this result, we assume in particular that the number of strata $S$, equivalently the cardinality of the empirical support of $Z$, is small compared to $n$. This condition fails if $Z$ has a continuous distribution. But it holds if the distribution of $Z$ is discrete, with a finite support or even an infinite one (as with, e.g., a Poisson distribution or a multivariate extension of it) provided that some moment of $Z$ is finite.
The main technical challenge for proving our asymptotic results is that as the sample size tends to infinity, there may be a growing number of strata, some being large and others small. We then consider separately “large” and “small” strata. For large strata, whose number may still tend to infinity, we establish a permutation central limit theorem using Stein's method CGS(2011). To this end, we derive a permutation version of the Marcinkiewicz-Zygmund inequality which, to our knowledge, is also new. For small strata, the combinatorial central limit theorem does not apply because the strata sizes may not tend to infinity. Instead, we use independence of units belonging to different strata, and the central limit theorem for triangular arrays.
We also study the performance of the SR test and other tests, through simulations. We show that the exactness of the SR test may endow it with a power advantage in some DGPs where other tests are underpowered, at least for small to moderate $n$. Otherwise, the SR test seems to have comparable power to other tests when $S/n$ is small. In the heteroskedastic case we explore, the SR test has a level closer to the nominal level than most of the other tests we consider. We also study in simulations an approximate version of the SR test, which may be useful when $Z$ is not discrete or $S/n$ is large. The approximate SR test is the SR test based on the strata obtained by discretizing the index $Z\widehat{\gamma}$, where $\widehat{\gamma}$ is the OLS estimator of $\gamma$. The approximate SR test displays some level of distortions, but still outperforms the partial correlation permutation test or the standard heteroskedasticity-robust test in some cases.
Finally, we consider two applications. The first studies the effect of some policies on traffic fatalities in the US, whereas the second revisits the effect of class size on students' achievement, using data from the project STAR (Student-Teacher Achievement Ratio). In both cases, confidence intervals on the effect of the evaluated policy obtained by inverting the SR test are informative. In the second application, the SR confidence intervals are always smaller than the usual, so-called HC3, confidence intervals. We also present evidence that inference based on the SR test may be preferable than that based on robust standard errors with these data.
{\bf Related literature.} As mentioned above, the paper is most closely related to DR2017. We extend their work by developing subvector inference that is exact under conditional independence between the covariates of interest and the unobserved term. Another related work is Lei-Bickel(2019), who develop a cyclic permutation test on subvectors. Their test is exact under a stronger independence condition than ours, but without any restriction on the regressors. On the other hand, they do not study the asymptotic behavior of their test under weaker conditions, and our simulations strongly suggest that it is not asymptotically valid under heteroskedasticity.
The idea of a stratified (or a restricted) randomization appears in the context of experimental designs Edgington(1983), Good(2013) and evaluation of treatment with randomly assigned instruments Imbens-Rosenbaum(2005) and treatments Bugni-Canay-Shaikh(2018). While the theoretical results developed in this paper are for observed units randomly sampled from some population, they are also directly applicable in experimental settings where units are randomly assigned to treatments as considered by the aforementioned studies.\footnote{Lehmann(1975) calls the former the population model and the latter the randomization model. Regarding the terms “permutation test" and “randomization test" which are often used interchangeably, Edgington-Onghena(2007) indicate that the former typically refers to permutation methods in population model while the latter is used for permutation methods in randomization model. Related to this distinction, see AAIW(2020) for obtaining standard errors in regressions in the presence of design-based and/or sampling-based uncertainties.} In particular, when applied to Imbens-Rosenbaum(2005)'s setup, our results allow one to fully characterize the asymptotic distribution of their test statistic without restricting the strata sizes and, at the same time, render it heteroskedasticity-robust.
Yet another approach to constructing exact permutation tests in the presence of nuisance parameters is to condition, whenever available, on a sufficient statistic for the nuisance parameters. This approach has been pursued by Rosenbaum(1984) for testing sharp null of no treatment effect under the logit assumption for the propensity score. However, in the absence of parametric assumptions as in our context, a low-dimensional sufficient statistic cannot be obtained in general. The approximate SR test we consider in the simulations is nonetheless related to Rosenbaum(1984), in the sense that it depends only on the index $Z\widehat{\gamma}$, rather than on the full set of regressors $Z$.
{\bf Organization of the paper.} The paper is organized as follows. Section (ref) introduces the set-up and develops the test. Section (ref) studies the finite-sample and asymptotic properties of the test. In Section (ref), we compare the performance of our test with alternative procedures through simulations. The two applications are considered in Section (ref), while Section (ref) concludes. The appendix gathers all proofs and additional results on the project STAR.
Consider (ref), where $y=[y_{1},\dots, y_{n}]^{\prime}$ is a $n\times 1$ vector of dependent variables, $X=[X_1,\dots, X_n]^{\prime}$ and $Z=[Z_1,\dots, Z_n]^{\prime}$ are $n\times k$ and $n\times p$ regressors, where $Z$ is assumed to include the intercept, and ${u}=[u_1,\dots, u_n]^{\prime}$ is a $n\times 1$ vector of exogenous error terms (see Assumptions (ref) and (ref)(ref) below). $\beta$ and $\gamma$ are $k\times 1$ and $p\times 1$ vector of unknown regression coefficients, respectively. We consider tests of the restriction\footnote{We thus do not consider tests on the intercept. These would require a different approach from that considered below.}
To formally define our test, we introduce the following notation. Let $W_i=(X_i^{\prime}, Z_i^{\prime})^{\prime}$, $W= [X, Z]$ and for any $n\times m$ matrix $A$, $M_{A}=I_{n}-A(A^{\prime}A)^{-1}A^{\prime}$. The set of all permutations of $\{1,\dots, n\}$ is denoted by $\mathbb{G}_n$, with $\text{Id}\in\mathbb{G}_n$ corresponding to the identity permutation. For any $\pi\in\mathbb{G}_n$ and vector $c=[c_1,\dots,c_n]'$, we let $c_\pi=[c_{\pi(1)},\dots,c_{\pi(n)}]'$. Similarly, for a matrix $A$ with $n$ rows, $A_\pi$ is the matrix obtained by permuting the rows according to $\pi$.
Our test differs from that of DR2017 in that instead of considering $\mathcal{W}^\pi$ for all $\pi\in\mathbb{G}_n$, we focus on a subset of $\mathbb{G}_n$. Let $\{z_1,\dots,z_S\}$ denote the set of (distinct) observed values in the sample. Let $n_s=\vert \{i:Z_i=z_s\}\vert$ for $s\in\{1,\dots,S\}$. Let also $y^s$ denote the subvector of $y$ including the rows $i$ satisfying $Z_i=z_s$, and define $X^s$ and $Z^s$ similarly. Without loss of generality, assume that the vector $y$ and matrices $X$ and $Z$ are arranged such that
Let the corresponding partitioning of the error terms be $u=[u^1,\dots, u^S]'$, and $\tilde{X}=[\tilde{X}^{1\prime},\dots, \tilde{X}^{S\prime}]'$, where $\tilde{X}^s=M_{\bm{1}_s}X^s$ and $\bm{1}_s$ denotes the $n_s\times 1$ vector of ones.\footnote{On the other hand, we do not modify the labels $1,\dots,n$ of the units. This means that $y\ne [y_1,\dots,y_n]'$, for instance ($y$ is a permuted version of $[y_1,\dots,y_n]'$).} Also, for any vector $v=[v_1,\dots,v_n]'$, let $\Sigma(v)$ denote the diagonal matrix with $(i,i)$ element equal to $[Dv]_i^2$, with $D\equiv \mathrm{diag}(M_{{\bm{1}_1}},\dots, M_{{\bm{1}_S}})$ and define
We consider the following heteroskedasticity-robust Wald statistic for $H_0$:
We now construct a permutation test based on $\mathcal{W}$. The idea behind permutations tests is that the distribution of some test statistic ($\mathcal{W}$ here) remains invariant under $H_0$ if we permute the data in an appropriate way. Then, we reject $H_0$ at the level $\alpha\in(0,1)$ if the test statistic on the initial data is larger than $100\times(1-\alpha)\%$ of the test statistics obtained on all possible permuted data.
First, we define the permuted version of $\mathcal{W}$ as
Second, we focus on permutations $\pi$ such that for any $k\in \{1,\dots,n\}$, we have, for all $t\in \{1,\dots,S\}$, $$\sum_{s=1}^{t-1} n_s < k\leq \sum_{s=1}^{t} n_s \Rightarrow \sum_{s=1}^{t-1} n_s < \pi(k)\leq \sum_{s=1}^{t} n_s \quad \left(\text{with } \sum_{s=1}^0 n_s=0\right).$$ That is, $\pi$ shuffles rows only within each of the $S$ strata, ensuring that $Z_\pi=Z$. We denote by $\mathbb{S}_n$ the set of such “stratified” permutations. Clearly, $\text{Id}\in\mathbb{S}_n$.
With the test statistic $\mathcal{W}$, its permuted version $\mathcal{W}^\pi$ and the set of admissible permutations $\mathbb{S}_n$ in hand, we define a level-$\alpha$ stratified randomization test by following the general construction of permutation tests. Let $N\in\{1,\dots,|\mathbb{S}_n|\}$, possibly random but independent of $y$ conditional on $W$ and let $\mathbb{S}'_n \subset \mathbb{S}_n$ be such that (i) $\text{Id}\in\mathbb{S}'_n$; (ii) $\mathbb{S}'_n\backslash \{\text{Id}\}$ is obtained by simple random sampling without replacement of size $N-1$ from $\mathbb{S}_n\backslash \{\text{Id}\}$. Note that if $N=|\mathbb{S}_n|$, we simply have $\mathbb{S}'_n = \mathbb{S}_n$.\footnote{ A practical way to obtain $\mathbb{S}'_n$ is (i) to pick $N'-1$ permutations at random from $\mathbb{S}_n$, with replacement and with equal probability, (ii) to add $\text{Id}$ to this initial set, (iii) to delete the duplicates (if any) from this set.} Then, let
be the order statistics of $(\mathcal{W}^\pi)_{\pi\in \mathbb{S}'_n}$. Let $q=N-\left \lfloor{N\alpha}\right \rfloor$ (with $\left \lfloor{x}\right \rfloor$ the integer part of $x$), $N^{+}=\vert \{i\in \{1, \dots, N\}: \mathcal{W}^{(i)}> \mathcal{W}^{(q)}\}\vert $ and $N^{0} =\vert \{i\in \{1, \dots, N\}: \mathcal{W}^{(i)}= \mathcal{W}^{(q)}\}\vert$. We define the level-$\alpha$ test function $\phi_\alpha$ by
Increasing $N$ reduces the role of randomness but increases the computational cost of the test. The exact result (Theorem (ref) below) holds for any $N$ but for the asymptotic results, we will assume that $N\stackrel{p}{\longrightarrow}\infty$ as $n\to\infty$.
The power of the test is directly related to $|\mathbb{S}_n|=\prod_{s=1}^S n_s!$. If $|\mathbb{S}_n|=1$, which occurs if $Z$ includes a continuous component,\footnote{This is at least the case if the data are i.i.d.. If not, we may have $S<n$ with positive probability even if $Z$ includes a continuous component.} the test becomes trivial: $\phi_\alpha=\alpha$. On the other hand, if $|\mathbb{S}_n|>1$, the test is non-trivial, and we may have $\phi_\alpha=1$ as soon as $|\mathbb{S}_n|\ge 1/\alpha$: this occurs if $\mathcal{W}>\max_{\pi \in\mathbb{S}_n:\pi\ne\text{Id}} \mathcal{W}^\pi$. For instance, $S=n-4$, $n_{s_1}=3$, $n_{s_2}=n_{s_3}=2$ for some $(s_1,s_2,s_3)$ and $n_s=1$ otherwise is sufficient to induce $|\mathbb{S}_n|\ge 1/\alpha$ for $\alpha\ge 0.05$. More generally, even if $Z$ has an infinite support, as with count data (e.g., Poisson distributed variables), $|\mathbb{S}_n|$ may be large, and $S/n$ small. We refer to Lemma (ref) below for a formal result along these lines.
As explained above, our test is trivial when $Z$ is continuous. More generally, it may also have low power with many strata. To circumvent these issues, we consider here an approximate version of our initial proposal. Specifically, instead of constructing strata based on $\{Z_i\}_{i=1}^n$, we rely on a discretization of $I_i\equiv Z'_i\widehat{\gamma}$, with $\widehat{\gamma}$ the OLS estimator of $\gamma$. Letting $u^1=\min_{i=1,\dots,n} I_i$, $u^{S+1} =\max_{i=1,\dots,n} I_i$ and $u^s = u^1 + (u^S-u^1) (s-1)/S$ for $s=2,\dots,S$, one such discretization is simply $$\sum_{s=1}^S s 1(I_i\in [u^{s}, u^{s+1})),$$ where $1(\cdot)$ denotes the indicator function. Intuitively, the number of strata $S$ that one chooses trades off size distortion and power. With a low $S$, the test is distorted because there are still variations in $\{I_i\}$ within each stratum, but we can expect larger power since the effective sample size $n-S$ is larger.
Because of the discretization, the test is not exact in general, even if $X$ and $u$ are conditionally independent. While we leave the study of its asymptotic validity for future research, we evaluate in Subsection (ref) its performances and the impact of the choice of $S$ through Monte Carlo simulations.
The confidence region for the parameters can be obtained by test inversion. Given the duality between tests and confidence sets, the finite and large sample validity of the confidence regions follow from the corresponding test results.
We recommend using the same set of permutations $\mathbb{S}'_n$ for different tested values $\beta_0$ and the same random variable in case randomization is required for the test (namely, when $\mathcal{W}=\mathcal{W}^{(q)}$ in (ref) above), to avoid random fluctuations that could create “holes” in the confidence region.\footnote{ If the test is trivial ($S=n$), confidence intervals obtained by test inversion may be empty. To avoid this issue, one can define the confidence interval as $(-\infty,\infty)$ in such cases, instead of randomizing.} Even if we do so, inverting the test may not lead to a convex confidence region. In this case, taking the convex hull of the corresponding region is a simple and conservative fix. This issue is not specific to our test but potentially arises with all permutation tests for which the critical value depends on the parameter value that is tested ($\beta_0$ in our context). For discussions about permutation confidence intervals and finding their endpoints, we refer to Garthwaite(1996), Good(2013) and Wang-Rosenberger(2020).
First, we prove that the SR test is exact under conditional independence between $X$ and $u$. We rely on the following conditions.
Because all observations $i$ in stratum $s$ are such that $Z_i=z_s$, the first condition allows for any dependence between $Z$ and $u$. The first two conditions are thus weaker than unconditional exchangeability of $u$ and independence between $u$ and $W$, the conditions imposed by Lei-Bickel(2019) and DR2017 (DR2017, test statistic $U_n(X,Y)$) to establish the exactness of their tests. A particular case where Condition (ref) holds is stratified randomized experiment with homogeneous treatment effects. Let $X_i\in\{0,1\}$ denote the treatment variable of individual $i$, $Z_i$ be the vector of strata dummies in and $Y_i(x)$ the potential outcome corresponding to treatment value $x\in\{0,1\}$. With homogeneous treatment effects, we have $Y_i(1)= \beta + Y_i(0)$ for all $i$, with $Y_i(0)\perp\!\!\!\perp X_i | Z_i$. Then, letting $u_i=Y_i(0)-\operatorname{E}[Y_i(0)|Z_i]$, we obtain (ref) and Condition (ref). Condition (ref) is maintained for convenience, as the Wald statistic ((ref)) and its permutation versions would still be exchangeable when defined using a generalized inverse of the covariance matrix estimator.
Hence, the SR test is exact in finite samples. In particular, in stratified randomized experiments, the test is exact for testing sharp null hypotheses of the kind $Y(1)- Y(0)= \beta$. By focusing on permutations in $\mathbb{S}_n$, we therefore extend the result of DR2017 to designs where $X$ and $Z$ are not independent. While the details of the proofs are given in the appendix, the intuition of the result is as follows. First, one can show that for the test to be exact, it suffices to prove that the $(\mathcal{W}^\pi)_{\pi \in\mathbb{S}'_n}$ are exchangeable, conditional on $W$. Now, if $H_0$ holds, $D(y-X\beta_0)=D u$. Then, because $v$ in $g(W,v)$ is always premultiplied by $D$, we have $\mathcal{W}= g(W, u)$. Next, for any $\pi\in\mathbb{S}_n$, we have, still under $H_0$, $$D(y-X\beta_0)_{\pi}=D(Z_\pi\gamma+u_\pi)=Du_\pi,$$ since $Z_\pi=Z$. Therefore, $\mathcal{W}^\pi = g(W, u_\pi)$. Conditional independence between $X$ and $u$ and exchangeability of $u^s$ for all $s$ then imply that the $(\mathcal{W}^\pi)_{\pi \in\mathbb{S}'_n}$ are exchangeable, conditional on $W$.
Next, we study the asymptotic validity of the SR test under weaker conditions than conditional independence between $X$ and $u$. To this end, we partly build on DR2017, who show that a randomization test based on the Wald statistic ($U_n(X,Y)$ in their paper) is of asymptotically correct level and heteroskedasticity-robust if, in addition to usual moment and nonsingularity conditions, the data are i.i.d. and $\operatorname{E}[W_iu_i]=0$.\footnote{{A close inspection of the proof of Theorem 4.1 in DR2017 reveals that their partial correlation test, to which we compare our test below, also works if $\operatorname{E}[W_i u_i]=0$.}} We present below analogous large sample results for the SR test under similar conditions, displayed in Assumption (ref) below.
Before we present these conditions, let us make some preliminary remarks and additional definitions. First, the distribution of $\{(W_i^{\prime}, u_i)^{\prime}\}_{i=1}^n$ and $\beta$ are allowed to change with $n$. For notational convenience, we usually do not make this dependence on $n$ explicit throughout the main text, though we do it in the appendix. Second, in order to derive the asymptotic distribution of the stratified randomization statistic, we make a distinction between “large” and “small” strata as follows. Let $c_n$ be a sequence satisfying $c_n\geq n^{1/2}$ and $c_n/n\to 0$ as $n\to \infty$ e.g. $c_n=n^{1/2}$. Define
We explain below on why we separate strata this way. Now, let $X^s=[X_{1}^s,\dots, X_{n_s}^s]'$ with $X_{i}^{s}$ being $k\times 1$, and $u^s=[u_{1}^s,\dots, u_{n_s}^s]$ $(n_s\times 1)$ denote the regressor matrix and the vector of error terms in the $s$th stratum respectively. The covariance matrix for the Wald statistic is defined as $$ \Omega_n \equiv n^{-1}\sum_{i=1}^{n} \operatorname{E}\left[(X_{i}-\operatorname{E}[X_{i}|Z_{i}])(X_{i}-\operatorname{E}[X_{i}|Z_{i}])'u_{i}^{2}|Z_{i}\right].$$ Define also the covariance matrices for the permutation statistic as ${V}_{n\mathcal{I}} \equiv \sum_{s\in \mathcal{I}_n}(n_s/n)\sigma_{u}^{s2}Q^s$ and $V_{n\mathcal{J}} \equiv \sum_{s\in \mathcal{J}_n}(n_s/n)\sigma_{u}^{s2}Q^s$, where $\sigma_{u}^{s2}\equiv n_s^{-1}\sum_{i=1}^{n_s}\operatorname{E}[u_{i}^{2}|Z_{i}=z_{s}]$ and $$Q^s \equiv n_s^{-1}\sum_{i=1}^{n_s}\operatorname{E}[X_{i}X_{i}^{\prime}|Z_{i}=z_{s}] -n_s^{-1}\sum_{i=1}^{n_s}\operatorname{E}[X_{i}|Z_{i}=z_{s}]\,n_s^{-1}\sum_{i=1}^{n_s}\operatorname{E}[X_{i}^{\prime}|Z_{i}=z_{s}].$$ Let $\lambda_{\min}(\cdot)$ denote the smallest eigenvalue of a symmetric matrix. We impose the following assumption on the strata.
Conditions (ref)-(ref) are quite standard. Condition (ref) is a usual moment condition required for the central limit theorem with independent observations White(2001),Hansen(2021). Condition (ref) is slightly stronger than the usual unconditional absence of correlation between $W_{i}$ and $u_{i}$ but weaker that the usual mean independence condition $\operatorname{E}[u_{i}|W_{i}]=0$, and, as such, allows for any form of heteroskedasticity.\footnote{ {An example where Condition (ref) holds but $\operatorname{E}[u_{i}|W_{i}]\ne 0$, suppose that $X_{i}$ is continuous, independent of $Z_{i}$ (as in a standard randomized experiment) and has a nonlinear effect on $Y_{i}$, so that $\operatorname{E}[Y_{i}|X_{i}]=g(X_{i})+Z_{i}'\tilde{\gamma}$ for some vector $\tilde{\gamma}$ and nonlinear function $g$. If we consider the best linear prediction of $Y_{i}$ by $(X_{i}', Z_{i}')'$, it will take the form $X_{i}'\beta+Z_{i}'\gamma$, and $u_{i}\equiv Y_{i}-X_{i}'\beta-Z_{i}'\gamma$ will satisfy Condition (ref), though $\operatorname{E}[u_{i}|W_{i}]\ne 0$.}} Consider again the case of stratified randomized experiments, but now assume that treatment effects can be heterogeneous, still with $\operatorname{E}[Y_i(1)-Y_i(0)|Z_i]=\operatorname{E}[Y_i(1)-Y_i(0)]$. Then, (ref) holds, with $\beta=\operatorname{E}[Y_i(1)-Y_i(0)]$. Besides, $\operatorname{Var}[u_{i}|W_{i}]$ depends on $X_{i}$ because of treatment effect heterogeneity. However, $\operatorname{E}[u_{i}|W_{i}]=0$, so Condition (ref) holds.
Conditions (ref)-(ref) are more specific to our setup. The first imposes that the covariance matrices corresponding to each large stratum are not degenerate. Condition (ref) also imposes the invertibility of the limit covariance matrix for $\mathcal{J}_n$ when the number of small strata is large. The latter can be replaced by the condition that $\lim_{n\to\infty}V_{n\mathcal{J}}$ exists; then the limit could be degenerate.
To understand Condition (ref), note that because the SR test uses deviations from strata means, the “effective” sample size $n-S$ should tend to infinity for the test to be consistent. Condition (ref) reinforces this requirement, since under the condition, $n-S=n(1-n^{-1}S)\stackrel{a.s.}{\longrightarrow} \infty$. Condition (ref) obviously holds if the $\{Z_i\}_{i=1}^n$ are identically distributed; with distribution independent of $n$, and $Z_1$ has finite support. The following lemma shows that it also holds if the support of $Z_{i}$ is in $\mathbb{N}^p$, as with, e.g., Poisson or geometric distributions (or multivariate versions of them), provided that some moments of $Z_{i}$ are finite.
For any finite set $B$, let $\mathcal{U}(B)$ denote the uniform distribution over $B$. To establish the asymptotic properties of the SR test, we first study the asymptotic behavior of $\mathcal{W}^\pi$, with $\pi \sim \mathcal{U}(\mathbb{S}_n)$, conditional on the data.
The main technical difficulty, and the reason why we cannot apply the same proof as in DR2017, is that the number of strata may tend to infinity, and there may be only small strata. To deal with these issues, we consider separately large and small strata. For large strata, we prove the following combinatorial central limit theorem with possibly many strata. Hereafter, we let $P^\pi$ denote the probability measure of $\pi$.
The proof of this lemma relies on Stein's method with exchangeable pairs, and on the following permutation version of the Marcinkiewicz-Zygmund inequality, which, to our knowledge, is also new.
The usual combinatorial central limit theorem, used in DR2017, would apply to a finite number of large strata; but here, the number of large strata, $|\mathcal{I}_n|$, may tend to infinity. We can accomodate that using Lemma (ref). Still, to check Condition (c) therein, we need to restrict the growth of $|\mathcal{I}_n|$, by imposing that $|\mathcal{I}_n|\le C n^{1/2}$ for some $C>0$ (see (ref) in the proof of Theorem (ref)). This is why we imposed above $c_n\ge n^{1/2}$.
For small strata, the combinatorial central limit theorem does not apply because the strata sizes may not tend to infinity. Instead, we use the fact that all strata are independent. Then, we check that the assumptions underlying a conditional version of the Lindeberg CLT for triangular arrays hold. This is technical, however, and we have to rely again on the Marcinkiewicz-Zygmund inequality, among other tools. Also, we rely therein on the condition that strata are small enough in the sense that $c_n/n\to 0$. A similar condition is used by Hansen-Lee(2019) to establish asymptotic normality of a sample mean with potentially many clusters.
The asymptotic properties of the SR test are also based on the the asymptotic behavior of $\mathcal{W}$, given in the following theorem. Hereafter, we let $\chi^2_k(c)$ be the noncentral chi-squared distribution with degrees of freedom $k$ and noncentrality parameter $c$, for any $c\ge 0$.
As Theorem (ref), Theorem (ref) would be standard with a finite number of strata; but here there may be many small strata. In particular, the result does not immediately follow from the result of Wooldridge(2001), which is derived under the assumption that each stratum frequency has a nondegenerate limit. To prove the result, we show that the assumptions underlying a conditional Lindeberg CLT hold, see Lemma (ref) in Appendix (ref).
The two previous results imply the following asymptotic properties of the SR test. Below, $q_{1-\alpha}(\chi^2_k)$ denotes the $1-\alpha$ quantile of the $\chi^2_k$ distribution.
The first result shows that the SR test is of asymptotically correct level. Combined with Theorem (ref), this result implies that under i.i.d. sampling and technical restrictions, the SR test is exact under conditional independence between $X_i$ and $u_i$, and asymptotically valid under the weaker exogeneity restriction $\operatorname{E}[(X_{i}^{\prime}, 1)' u_{i}|Z_{i}=z_{s}]=0$. The second result in Corollary (ref) states that the test has nontrivial power to detect local alternatives. Finally, the third result implies that if $\operatorname{E}[n^{-1} \tilde{X}\tilde{X}']$ converges to symmetric positive definite matrix, the test is consistent for fixed alternatives, namely when $\beta_n=\beta\ne \beta_0$.
We first present some simulation evidence on the performance of the proposed test. We consider the following model:
where $\beta=\gamma_1=0$, $\gamma_j=1$ for $j\geq 2$. We consider $p\in\{2,4\}$ and sample sizes $n\in\{50,100,500\}$. Then, we consider three data generating processes (DGPs) with various distributions for $X_i$ and $u_i$. In the three cases, the $\{Z_{ij}\}_{j=2}^{p}$ are i.i.d. and follow a Poisson distribution with parameter 1. Then:
Remark that in the first two DGPs, $u_i$ is independent of $X_i$ given $Z_i$ and so the SR test has exact size. In DGP3, $u_i$ is heteroskedastic and the SR test may exhibit a finite-sample size distortion. As seen below, DGP3 leads to over-rejection of usual heteroskedasticity-robust test in finite samples. We aim at investigating the properties of the SR test (and others) in these situations.
Table (ref) shows some characteristics of $W$ related to the SR test. With one Poisson regressor, the largest stratum roughly corresponds to 40% of the whole sample, which implies that the number of distinct permutations in $\mathbb{S}_n$ ($\operatorname{E}[|\mathbb{S}_n|]$) is very large, even with $n=50$. With three Poisson regressors, on the other hand, many strata have just size one, and $S$ is quite large compared to $n$. It is then interesting to investigate the power of the test in this more difficult case.\footnote{ When increasing $p$ further, e.g. to 10, strata have only size 1, and the SR test becomes trivial in this case.}
We consider the power curve on the interval $[-0.5, 0.5]$. We construct $\mathbb{S}'_n$ as explained in Footnote (ref) above, with $N'=499$. However, if $|\mathbb{S}_n|<499$, we simply consider all possible permutations, so that $\mathbb{S}'_n=\mathbb{S}_n$.
We then compare our test with an asymptotic version of it (denoted by SRa in the figures below), where the test statistic is the same but we use the $\chi^2_1$ asymptotic critical value instead of the distribution of the permutation statistic. Apart from the SR test, we consider the cyclic permutation (CP) test of Lei-Bickel(2019), which is exact under independence between $W$ and $u$ (but not conditional independence between $X$ and $u$ given $Z$). Following Lei-Bickel(2019) (p.406), we use $19$ cyclic permutation samples and a stochastic algorithm to find an ordering of the data that improves power. We also consider the partial correlation permutation test, denoted as PC, proposed by DR2017, which is asymptotically valid under simply no correlation between $W$ and $u$. Finally, we consider the sign test of Toulis(2022), which is asymptotically valid under symmetry of $u$, which holds in the three DGPs. For these two permutation tests, we use 499 permutations drawn at random with replacement from $\mathbb{G}_n$. We also compare our test with the usual F-test, which is not heteroskedasticity-robust. We use the $F_{1, n-p-1}$ critical value for this test. Finally, we compare our test with the heteroskedasticity-robust Wald tests. For the latter, we use the so-called HC1 and HC3 versions of the test. The former is the default in Stata and is widely used, whereas HC3 is often the recommended version of the test long2000using. In both tests, we use the $\chi^2_1$ critical value.
In the first DGP and with $p=2$, the performances of the permutation tests are overall similar (see Figure (ref)). We just note that the CP test is slightly less powerful than the others. Compared to standard tests, the SR test is slightly less powerful but no difference can be detected when $n=500$. With $p=4$ and $n<500$, the SR test is less powerful than the PC and sign tests. This could be expected because basically, the test relies only on $n-S$ observations, and for $p=4$ and $n\in\{50,100\}$, $n-S$ is much smaller than $n$ (see Table (ref) above).\footnote{ The tails of $X$ also seem to affect the power of the SR test. Considering $X_i=\exp(X_i^*)$ instead of $X_i=X_i^*$, the power of the SR deteriorates compared to that of the CP and sign test (but improves compared to the PC test, at least for $p=2$).} However, the difference in power becomes very small when $n=500$. Interestingly, the CP test does not seem to have power for $n=50$. We also notice that in this yet homoskedastic model, the HC1 and HC3 Wald tests overreject when $n=50$.
The second DGP is an instance where the SR test is the only exact one, among all the tests we consider. Even without size correction, it also has generally larger power than the PC test, except for $n=50$ and $p=4$ (see Figure (ref)). The sign test has better power (especially when $n\le 100$) but is also distorted, with levels around 10% for $n\le 100$. The HC1 and HC3 tests exhibit large distortions. With $p=2$, they respectively reject the null hypothesis in 33.0% and 16.1% of the samples with $n=100$, and in 25.5% and 15.7% of the samples with $n=500$.
\thispagestyle{empty}
The last DGP corresponds to a case where no test is exact because of heteroskedasticity, and all the tests we consider overreject in finite samples (see Figure (ref)). Overall, the SR test exhibits reasonable level of distortion, with a level that never exceeds 7.4% over the six combinations of $n$ and $p$, and equal to 6.6% on average. Other tests have average rejection rates of 52.4% (CP), 5.7% (PC), 17.1% (sign), 30.3% (HC1), 11.2% (HC3) and 62.8% (non-robust). The CP test exhibits high level of distortions, which also increase with $n$, confirming our conjecture that it is not heteroskedasticity-robust. Compared to the PC test, which has a similar level, the SR test has a larger power both when $p=2$ and $p=4$.
Of course, the ranking between the different tests in terms of level and power may vary for other DGPs. Nevertheless, the simulations above suggest that the SR test can be a good competitor to existing tests, especially in models with few additional covariates.
We now investigate the approximate version of our SR test which handles continuous $Z$s, see Subsection (ref) above. To this end, we consider the same DGPs as above but now assume that the $\{Z_{ij}\}_{j=2}^{p}$ are i.i.d. and follow a standard normal distribution, instead of a Poisson distribution. We refer to the corresponding DGPs as DGP1', DGP2' and DGP3'. We also have to choose the number of strata $S$. As mentioned above, this involves a trade-off between size distortion and power. Also, one should consider a larger number of strata when the correlation between $X_i$ and $Z_i'\gamma$ is high, since then the correlation between $y_i-X_i\beta_0$ and $X_i$ within strata becomes larger. Guided by this, we use the following data-driven $S$:
where $\lceil x\rceil$ is the smallest integer greater than or equal to $x$ and $\widehat{\text{corr}}$ denotes the empirical correlation coefficient. The rule in (ref) generally leads to small strata sizes, around four in DGP1' and five in DGP2' and DGP3', as Table (ref) shows.
The power curves are displayed in Figures (ref)-(ref); here we compare our test with the PC and HC3 tests. In DGP1', the test exhibits almost no distortion. It is slightly less powerful than the two other tests, but the difference gets attenuated as $n$ increases. In DGP2', the test is hardly distorted with $p=2$, with a level of at most 5.2%, and has, in general, better power than the PC test. It exhibits some distortion when $p=4$ and $n=50$, with a level of 7.4%, but this distortion quickly vanishes as $n$ increases. Its power is higher than that of the PC test for $\beta<0$ and $n=100$ but slightly lower otherwise. In DGP3', the test appears to have good power compared to the PC and HC3 tests, even with $p=4$. It is also less distorted than the HC3 test but more than the PC test, with an average over the six cases of 7.1% vs respectively 9.9% and 4.7%. Finally, additional simulations with higher $p$ (e.g., $p=10$), not presented here, suggest that the approximate SR test is not really sensitive to the dimension of $Z$.
Overall, the results suggest that with the data-driven choice of $S$ above, the test is asymptotically valid and consistent, though we leave this question for future research.
We first apply the proposed randomization inference method to analyse the effect of driving regulations on traffic fatalities in the US, using the same data and model as Wooldridge(2015) (Wooldridge(2015), Chapter 13). Specifically, we consider the linear model
where for any US state $i$, year $t$ and random variable $A_{i,t}$, we let $\Delta A_i=A_{i,1990}-A_{i,1985}$. In (ref), $dthrte$ denotes traffic fatality rate, $open$ is a dummy variable for having an open container law, which illegalizes for passengers to have open containers of alcoholic beverages, and $admn$ is a dummy variable for having administrative per se laws, allowing courts to suspend licenses after a driver is arrested for drunk driving but before the driver is convicted. The OLS estimate is $(\widehat{\beta},\widehat{\gamma})= (-0.42,-0.15)$, pointing towards a deterrent effect of the two types of laws on alcohol consumption by drivers.
The coefficient $\gamma$ is never significant for any usual level, so we focus below on $\beta$. We compute the confidence intervals based on the SR test inversion (SR confidence interval hereafter) and those based on the CP and PC test inversion. We also consider the inversion of the test of the full vector of parameters, considered by DR2017 in their Section 3 (PR confidence interval hereafter). By projecting the corresponding confidence region over $\beta$, we obtain a confidence interval that is conservative under independence, or asymptotically conservative under weaker conditions. Finally, we consider the standard, non-robust and robust confidence intervals.
For all confidence intervals based on permutation tests, we invert the tests of $\beta=\beta_0$ for $\beta_0\in\{-1.7,-1.69,\dots,0.3\}$. As recommended above, we use the same set of permutations for all values of $\beta_0$ that we test. This way, we obtain proper intervals for the four permutation methods. We draw $\mathbb{S}'_n$ as explained in Footnote (ref). For the other tests, we draw uniformly and with replacement $N'$ permutations from $\mathbb{G}_n$. To limit the effect of randomness, we use $N'=99,999$ instead of 499 as in the simulations.\footnote{In this application, $n=51$, $S=3$ with $\max_{s} n_s=41$, so the set $\mathbb{S}_n$ is large ($|\mathbb{S}_n|\simeq 1.2\times 10^{55}$).} Even though we invert 201 tests and use a large number of permutations, the SR and PC confidence intervals take on our computer just a few seconds to compute (see the last column of Table (ref) for computational times). The CP and PR confidence intervals, on the other hand, are more computationally intensive. The reason for the CP method is that improving its power requires finding an optimal ordering, a difficult optimization problem Lei-Bickel(2019). The PR confidence interval is very costly to compute because it requires testing for values of $(\beta,\gamma)$, and thus considering a grid in $\mathbb{R}^2$ instead of $\mathbb{R}$.
The results are reported in Table (ref). We consider confidence intervals with nominal levels of 90% and 95%. The CP confidence intervals are by far the largest. Not surprisingly, the PR confidence interval is also large for the 95% nominal level confidence interval, though it remains much shorter than the CP confidence interval. And actually, it is close to the SR, PC and robust confidence intervals when considering a nominal level of 90%. For both nominal levels, the SR and PC confidence intervals are very close. They are also close to the robust confidence interval for the nominal level of 90%, but around 11% larger than this confidence interval for the nominal level of 95%. Finally, the non-robust confidence intervals are the shortest of all intervals, being roughly 13% smaller than the robust confidence intervals.
The SR confidence interval may be more reliable than the robust confidence interval. To see why, note first that there is no evidence of heteroskedasticity: the White and Breusch-Pagan tests have p-values of respectively 0.56 and 0.34, respectively. If independence holds, the SR confidence interval has exact coverage, whereas the robust confidence interval may exhibit some distortion. To evaluate this, we ran simulations, assuming that $u_i$ is independent of $W_i$ and is distributed according to the empirical distribution of the residuals $\widehat{u}_i$. For nominal coverage of 95% and 90%, the robust confidence interval includes the true parameter in only 89.7% and 85.7% of the samples, respectively. Finally, there is strong evidence of non-normal errors: the Shapiro-Wilk and Jarque-Bera tests of normality on residuals have p-values of 0.0015 and 0.0012, respectively. As a result, the non-robust confidence intervals based on normality may also exhibit distortion. If one drops the normality assumption but maintains independence, the SR confidence interval is the only one with exact coverage. Then, one cannot exclude even at the 10% level that the open container law has no effect on traffic fatalities.
Finally, we apply the SR test to the well-known dataset of the Project STAR experiment. Imbens-Rubin(2015) (Chapter 9) provide a detailed analysis of the data using several stratified randomization-based inference methods. We follow their regression analysis of stratified randomized experiments (Chapter 9.6) and focus on schools with at least two regular classes and two small classes ignoring classes with teacher's aides. In the specifications considered below, there are at most 25 such schools that define the strata. The majority of schools have exactly 4 classes, and only a few have 5 or more classes.
The regression model considered by Imbens-Rubin(2015) is as follows:
where the treatment variable $X_{i}$ is the indicator for small classes, and the nuisance regressors $Z_{i}^s=1(i\in s), s=1,\dots, S,$ are the school (strata) indicators, and the outcome variable is the class-level (teacher-level) average math test scores for kindergarten children. We consider the class-level average reading test scores in addition to the math scores, and Grade 1 and Grade 2 as well. As pointed out by Imbens-Rubin(2015), restricting the analysis to class-level data avoids a possible violation of the no-interference part of the Stable Unit Treatment Assumption (SUTVA).
We deviate in a minor way from Imbens-Rubin(2015)'s analysis by not standardizing the outcome variables to have mean 0 and standard deviation 1, because doing so introduces a slight dependence in the observations although the exchangeability of $u^s, s=1,\dots, S,$ would still be preserved. The results for standardized test scores are nevertheless similar, see Table (ref) in Appendix (ref).\footnote{Also, we were unable to obtain exactly the same sample, and thus the same results, as Imbens-Rubin(2015). There are 66 classes (34 small and 32 regular) and 15 schools (strata) in our sample, as opposed to 68 classes (36 small and 32 regular) and 16 schools in theirs. However, our estimate of $\beta$ (0.22) and standard error (0.09), calculated following the variance formula in Theorem 9.1 of Imbens-Rubin(2015), are close to theirs Imbens-Rubin(2015).} The tests implemented are the same as those in Table (ref), except that we include the HC0 confidence interval, denoted as Robust HC0 (IR),\footnote{The HC0 variance estimate is numerically identical to a sample analog of the variance formula in Theorem 9.1 of Imbens-Rubin(2015).} but do not include the projection-based test PR. The latter is computationally prohibitive in the current application, as there are at least 15 regression coefficients not under the test. We invert the tests of $\beta=\beta_0$ for $\beta_0\in\{-40,-39.99,\dots, 39.99, 40\}$ using $N'=99,999$ for each point.
The 95% confidence intervals are reported in Table (ref). The number of possible permutations is at least $\vert \mathbb{S}_n\vert=9.47\times 10^{24}$ and $S/n$ is at most $25/109=0.23$ in the specifications. As a general pattern, we can notice that the CP confidence intervals include all the tested points in all of the specifications and the robust HC3 confidence intervals are the second widest, while the HC0 and PC confidence intervals are the shortest. The SR confidence intervals, though wider than the HC0 and PC confidence intervals, are comparable with the non-robust confidence intervals and always shorter than the HC3 confidence intervals.
The baseline results for the math test scores of kindergarten children analyzed by Imbens-Rubin(2015) are noteworthy. In this case, there is no evidence of either heteroskedasticity or non-normality. The Breusch-Pagan test for heteroskedasticity has p-value 0.32, so homoskedasticity assumption is supported at the conventional significance levels. The Jarque-Bera and Shapiro-Wilk tests for the residuals have p-values 0.89 and 0.53, respectively, pointing towards Gaussian errors. As such, the tests with finite sample validity i.e. the non-robust and SR tests, should be more reliable. And in fact, they turn out to be very similar: if anything, the non-robust confidence interval is slightly longer. It also includes 0, contrary to the SR confidence interval. The HC0 and PC confidence intervals also show a significant treatment effect.
For reading scores in kindergarten, all test results suggest that the class-size reduction has no signifcant effect, in contrast with the results for math test scores. But class size reduction does seem to have an effect on both test scores in Grade 1 and Grade 2: all tests suggest significant treatment effects. As in the first application, the Breusch-Pagan and normality tests point to homoskedasticity and non-normality of the error terms in Grade 2, in which case the SR test could be the most reliable.
Inasmuch as the permutation tests should be more reliable than the non-permutation tests in experimental datasets such as the current one, and the analysis using teacher-level samples guards effectively against a possible spillover in the students' performance, the results suggest that the class-size reduction has a small but significant effect on the math test scores for kindergarten children, and a bigger effect on the math and reading test scores for Grade 1 and Grade 2 students but no effect on the reading test scores for kindergarten children at teacher-level.
We develop a new permutation test for subvector inference in linear regressions. The test has exact size in finite samples if the error terms $u_i$ are independent of the regressors of interest $X_i$, conditional on other regressors $Z_i$. If independence fails but $\operatorname{Cov}(X_i,u_i|Z_i)=0$, the test remains asymptotically valid with power against local alternatives under some conditions. The main one is that the number $S$ of distinct rows of $(Z_1,\dots,Z_n)'$ is negligible compared to the sample size $n$. Monte Carlo simulations suggest that the test has good power compared to other tests when, indeed, $n^{-1}S$ is small, and that it exhibits limited distortion without conditional independence. The two applications confirm that in some realistic designs, the test is informative and can thus be an appealing alternative to existing methods.
A few questions are left for future research. First, some simulations we conducted (not reported above) suggest that the condition $n^{-1}S\to 0$ could be replaced by the weaker condition that the effective sample size $n-S$ tends to infinity. Second, while we show that the test is asymptotically valid if $\operatorname{Cov}(X_i,u_i|Z_i)=0$, we do not establish finite sample guarantees in this set-up.\footnote{ See Theorem 2 in Toulis(2022) for an example of such guarantees. Note however that his result is obtained under conditions for which our test is actually exact.} Finally, constructing a permutation test for subvectors that is both exact under independence and asymptotically heteroskedasticity-robust for any design remains an important challenge.