An Assumption-Free Exact Test For Fixed-Design Linear Models With Exchangeable Errors

TL;DR

The paper introduces the Cyclic Permutation Test (CPT), an exact, assumption-free method for testing linear hypotheses in fixed-design linear models with exchangeable errors, valid in finite samples and applicable to arbitrary designs.

Contribution

It presents a novel non-randomized test that guarantees exact Type I error control and extends existing methods by using cyclic permutations and solving a traveling salesman problem for enhanced power.

Findings

CPT achieves comparable power to existing tests in simulations.

The test provides exact confidence intervals for single coefficient contrasts.

It is valid for arbitrary fixed designs and exchangeable errors in finite samples.

Abstract

We propose the Cyclic Permutation Test (CPT) to test general linear hypotheses for linear models. This test is non-randomized and valid in finite samples with exact Type I error $\alpha$ for an arbitrary fixed design matrix and arbitrary exchangeable errors, whenever $1 / \alpha$ is an integer and $n / p \ge 1 / \alpha - 1$. The test involves applying the marginal rank test to $1 / \alpha$ linear statistics of the outcome vector, where the coefficient vectors are determined by solving a linear system such that the joint distribution of the linear statistics is invariant with respect to a non-standard cyclic permutation group under the null hypothesis.The power can be further enhanced by solving a secondary non-linear travelling salesman problem, for which the genetic algorithm can find a reasonably good solution. Extensive simulation studies show that the CPT has comparable power to existing tests. When testing for a single contrast of coefficients, an exact confidence interval can be obtained by inverting the test. Furthermore, we provide a selective yet extensive literature review of the century-long efforts on this problem, highlighting the novelty of our test.

Equations

$$y_i={\beta }_0+\sum _{j=1}^p{x}_{ij}{\beta }_j+{ϵ}_i,i=1,\dots ,n,$$

$$y={\beta }_0\text{1}+X\beta +ϵ,$$

$$H_0:R^T\beta =0,\text{mbox}whereR\in {\mathbb{R}}^{p\times r}\text{mbox}isafixedmatrixwithrankr.$$

$$H_0:{\beta }_1=\dots ={\beta }_r=0.$$

$$y={\beta }_0\text{1}+X{U}_R(U_R^T\beta )+X{V}_R(V_R^T\beta )+ϵ={\beta }_0\text{1}+\sum _{j=1}^r{\tilde{X}}_j{\tilde{\beta }}_j+\sum _{j=r+1}^p{\tilde{X}}_j{\tilde{\beta }}_j+ϵ.$$

$$S\stackrel{d}{=}{\pi }_L(S)\stackrel{d}{=}{\pi }_L^2(S)\stackrel{d}{=}\cdots \stackrel{d}{=}{\pi }_L^m(S),$$

$$\text{mbox}where{\pi }_L^k(S)=(S_k,S_{k+1},\dots ,S_m,S_0,S_1,\dots ,S_{k-1}),k=1,2,\dots ,m.$$

$$\mathcal{L}(R_0)⪰\mathrm{Unif}([m+1])⟹\text{ if }p\triangleq \frac{R_0}{m+1},\mathcal{L}(p)⪰\mathrm{Unif}([0,1])$$

$${\tilde{S}}_j=|S_j-\text{mbox}median\{(S_j{)}_{j=0}^m\}|.$$

$$\tilde{S}=\{g(S_0;S),g(S_1;S),\dots ,g(S_m;S)\}$$

$$g(z;w)=g(z;{\pi }_{L}w).$$

$$S_j=y^T{\eta }_j,j=0,1,\dots ,m,$$

$$y^T{\eta }_j=({\text{1}}^T{\eta }_j){\beta }_0+(X_{[r]}^T{\eta }_j{)}^T{\beta }_{[r]}+(X_{[-r]}^T{\eta }_j{)}^T{\beta }_{[-r]}+{ϵ}^T{\eta }_j.$$

$$({ϵ}_1,\dots ,{ϵ}_n)\stackrel{d}{=}({ϵ}_{\pi (1)},\dots ,{ϵ}_{\pi (n)}).$$

$$y^T{\eta }_j=\underset{\text{deterministic part}}{\underbrace{({\text{1}}^T{\eta }_j){\beta }_0+(X_{[-r]}^T{\eta }_j{)}^T{\beta }_{[-r]}}}+\underset{\text{stochastic part}}{\underbrace{{ϵ}^T{\eta }_j}}.$$

$$X_{[-r]}^T{\eta }_j={\gamma }_{[-r]},(j=0,1,\dots ,m).$$

$${\eta }_0=({\eta }_1^{\ast },{\eta }_2^{\ast },{\eta }_3^{\ast },{\eta }_4^{\ast },{\eta }_5^{\ast },{\eta }_6^{\ast }{)}^T,{\eta }_1=({\eta }_3^{\ast },{\eta }_4^{\ast },{\eta }_5^{\ast },{\eta }_6^{\ast },{\eta }_1^{\ast },{\eta }_2^{\ast }{)}^T,{\eta }_2=({\eta }_5^{\ast },{\eta }_6^{\ast },{\eta }_1^{\ast },{\eta }_2^{\ast },{\eta }_3^{\ast },{\eta }_4^{\ast }{)}^T.$$

$$\displaystyle(\epsilon^{T}\eta_{0},\epsilon^{T}\eta_{1},\epsilon^{T}\eta_{2})^{T}=\left(\begin{array}[]{cccccc}\epsilon_{1}&\epsilon_{2}&\epsilon_{3}&\epsilon_{4}&\epsilon_{5}&\epsilon_{6}\\ \epsilon_{5}&\epsilon_{6}&\epsilon_{1}&\epsilon_{2}&\epsilon_{3}&\epsilon_{4}\\ \epsilon_{3}&\epsilon_{4}&\epsilon_{5}&\epsilon_{6}&\epsilon_{1}&\epsilon_{2}\end{array}\right)\eta^{*},$$

$$\displaystyle(\epsilon^{T}\eta_{1},\epsilon^{T}\eta_{2},\epsilon^{T}\eta_{0})^{T}=\left(\begin{array}[]{cccccc}\epsilon_{5}&\epsilon_{6}&\epsilon_{1}&\epsilon_{2}&\epsilon_{3}&\epsilon_{4}\\ \epsilon_{3}&\epsilon_{4}&\epsilon_{5}&\epsilon_{6}&\epsilon_{1}&\epsilon_{2}\\ \epsilon_{1}&\epsilon_{2}&\epsilon_{3}&\epsilon_{4}&\epsilon_{5}&\epsilon_{6}\end{array}\right)\eta^{*}.$$

$$\left(\begin{array}[]{cccccc}\epsilon_{1}&\epsilon_{2}&\epsilon_{3}&\epsilon_{4}&\epsilon_{5}&\epsilon_{6}\\ \epsilon_{5}&\epsilon_{6}&\epsilon_{1}&\epsilon_{2}&\epsilon_{3}&\epsilon_{4}\\ \epsilon_{3}&\epsilon_{4}&\epsilon_{5}&\epsilon_{6}&\epsilon_{1}&\epsilon_{2}\end{array}\right)\stackrel{{\scriptstyle d}}{{=}}\left(\begin{array}[]{cccccc}\epsilon_{5}&\epsilon_{6}&\epsilon_{1}&\epsilon_{2}&\epsilon_{3}&\epsilon_{4}\\ \epsilon_{3}&\epsilon_{4}&\epsilon_{5}&\epsilon_{6}&\epsilon_{1}&\epsilon_{2}\\ \epsilon_{1}&\epsilon_{2}&\epsilon_{3}&\epsilon_{4}&\epsilon_{5}&\epsilon_{6}\end{array}\right)$$

$$⟹({ϵ}^T{\eta }_0,{ϵ}^T{\eta }_1,{ϵ}^T{\eta }_2)\stackrel{d}{=}({ϵ}^T{\eta }_1,{ϵ}^T{\eta }_2,{ϵ}^T{\eta }_0).$$

$${\eta }_j={\pi }_L^{tj}({\eta }^{\ast })$$

$${\eta }_j={\left[{\pi }_L^{tj}\{({\eta }_1^{\ast },\dots ,{\eta }_{(m+1)t}^{\ast })\},{\eta }_{(m+1)t+1}^{\ast },\dots ,{\eta }_n^{\ast }\right]}^T,$$

$$\left(\begin{array}[]{ll}-I_{p-r}&X_{[-r]}^{T}\\ -I_{p-r}&X_{[-r]}^{T}\Pi_{1}\\ \vdots&\vdots\\ -I_{p-r}&X_{[-r]}^{T}\Pi_{m}\end{array}\right)\left(\begin{array}[]{c}\gamma_{[-r]}\\ \eta_{*}\end{array}\right)=0.$$

$$n/(p-r)>m.$$

$$(y^T{\eta }^{\ast },y^T{\Pi }_1{\eta }^{\ast },\cdots y^T{\Pi }_m{\eta }^{\ast })$$

$$n>19(p-r).$$

$$y^T{\eta }_j=(X_1^T{\eta }_j){\beta }_1+W_j$$

$$X_1^T{\eta }_j={\gamma }_1(j=1,2,\dots ,m),X_1^T{\eta }_0={\gamma }_1+\delta .$$

$$\bigg{(}-e_{1,p(m+1)}\,\,\vdots\,\,A(X)^{T}\bigg{)}\left(\begin{array}[]{l}\delta\\ \gamma\\ \eta\end{array}\right)=0,$$

$$A(X)=\left(\begin{array}[]{cccc}-I_{p}&-I_{p}&\cdots&-I_{p}\\ X&\Pi_{1}^{T}X&\cdots&\Pi_{m}^{T}X\end{array}\right)\in\mathbb{R}^{(n+p)\times p(m+1)}.$$

Full text

Keywords: assumption-free test, exact test, fixed-design, linear model, linear hypothesis, marginal rank test, non-linear travelling salesman problem

1 Introduction

In this article, we consider the following fixed-design linear model

[TABLE]

where the $\epsilon_{i}\,(i=1,\ldots,n)$ are stochastic errors and the $x_{ij}\,(i=1,\ldots,n,j=1,\ldots,p)$ are treated as fixed quantities. Throughout we will use the following compact notation

[TABLE]

where $y=(y_{1},\ldots,y_{n})^{T}$ denotes the response vector, $X=(x_{ij})\in\mathbb{R}^{n\times p}$ denotes the design matrix, $\epsilon=(\epsilon_{1},\ldots,\epsilon_{n})^{T}$ denotes the error terms and $\textbf{1}\in\mathbb{R}^{n}$ denotes the vector with all entries equal to one. Our focus is on testing a general linear hypothesis:

[TABLE]

Testing linear hypotheses in linear models is ubiquitous and fundamental in numerous areas. One important example is to test whether a particular coefficient is zero, i.e. $H_{0}:\beta_{1}=0$, a special case where $R=(1,0,\ldots,0)^{T}\in\mathbb{R}^{p\times 1}$. Another important example is to test the global null, i.e. $H_{0}:\beta=0$, equivalent to the linear hypothesis with $R=I_{p\times p}$. We refer to Chapter 7 of Lehmann & Romano (2006) for an extensive discussion of other examples. By inverting a test with valid Type I error control, we can obtain a confidence interval/region for $R^{T}\beta$. This is of particular interest when $r=1$, which corresponds to a single linear contrast of the regression coefficient.

Testing linear hypotheses in linear models is one of the most fundamental and long-lasting problems in statistics, as well as a convenient powerful prototype to motivate methods for more complicated statistical problems. In the past century, several types of methods have been proposed: normal theory-based tests (Fisher, 1922, 1924), permutation tests (Pitman, 1937b, 1938), rank-based tests (Friedman, 1937), tests based on regression R-estimates (Hájek, 1962), M-estimates (Huber, 1973) and L-estimates (Bickel, 1973), resampling-based tests (Freedman, 1981), median-based tests (Theil, 1950a; Brown & Mood, 1951), symmetry-based tests (Hartigan, 1970) and non-standard tests (Meinshausen, 2015). Here we list only the earliest reference we could track down for each category to highlight the chronology of the methodological development; an extensive literature review is provided in Appendix B.

For a given confidence level $1-\alpha$, a test is exact if the Type I error is below or equal to $\alpha$, in finite samples without any asymptotics. Exact tests are intellectually and practically appealing because they provide strong error control without the requirement of a large sample or artificial asymptotic regimes. However, perhaps surprisingly, there is no test that is exact under reasonably general assumptions to the best of our knowledge. A brief summary of the conditions under which the existing tests are exact is as follows:

• •

Regression t- and F-tests are exact with normal errors;

• •

Permutation tests are exact for the global null or certain null hypotheses for certain analysis of variance (ANOVA) problems (e.g. Brown & Maritz, 1982);

• •

Rank-based tests are exact for ANOVA problems;

• •

Tests based on regression R-, M- or L-estimates can be made exact for the global null;

• •

Hartigan (1970)’s test is exact for certain forms of balanced ANOVA problems with symmetric errors and $r=1$;

• •

Meinshausen (2015)’s test is exact for rotationally invariant errors with known noise level, and if the $\epsilon_{i}$s are independent and identically distributed (i.i.d.), rotation invariance implies the normality of the $\epsilon_{i}$ (Maxwell, 1860);

• •

Other tests are exact either for the global null or under restrictive assumptions or require excessive computation.

In this article we develop an exact test, which we refer to as the Cyclic Permutation Test (CPT), that is valid in finite samples, and can accommodate an arbitrary fixed design matrix and arbitrary error distributions, provided that the error terms are exchangeable. Exchangeability is weaker than the frequently made assumption of i.i.d. random variables. Further, the test is non-randomized if $1/\alpha$ is an integer and $n/(p-r)>1/\alpha-1$. The former condition is true for all common choices of $\alpha$, e.g. $0.1,0.05,0.01,0.005$. The latter requirement is also reasonable in various applications. For instance, when $\alpha=0.05$, the condition reads $n/(p-r)>19$, which is true if $n/p>19$ or $p-r$ is small; both are typical in social science applications. Admittely, it may be stringent in areas like genetics where $p$ is often larger than $n$. However, valid inference, or even identification, in those problems would require extra assumptions on the sparsity of $\beta$, geometry of $X$, and distribution of $\epsilon$, which are not in accordance with the goal of this paper to develop assumption-free tests. We demonstrate the power of the CPT through extensive simulation studies and show it is comparable to the existing ones. Although exchangeability may not be valid in certain applications, the CPT is the first procedure that is provably exact with reasonable power under such weak assumptions. We want to emphasize that the goal of this paper is not to propose a procedure that is superior to existing tests, but rather to expand the toolbox of exact inference and, hopefully, motivate the development of novel methods for other problems.

2 Cyclic Permutation Test

2.1 Main idea

Throughout the article we denote the set $\{1,\ldots,n\}$ by $[n]$. First we show that it is sufficient to consider the sub-hypothesis:

[TABLE]

In fact, let $U_{R}\in\mathbb{R}^{p\times r}$ be an orthonormal basis of the column span of $R$ and $V_{R}\in\mathbb{R}^{p\times(p-r)}$ be an orthonormal basis of the orthogonal complement. Then $\beta=U_{R}U_{R}^{T}\beta+V_{R}V_{R}^{T}\beta$. Let $\tilde{X}=(XU_{R}\,\vdots\,XV_{R})$, where $\vdots$ marks the partition of columns, and $\tilde{\beta}^{T}=(\beta^{T}U_{R},\beta^{T}V_{R})$. Then the linear model (1) can be re-formulated as

[TABLE]

On the other hand, since $R$ has full column rank, the null hypothesis (2) is equivalent to $H_{0}:\tilde{\beta}_{1}=\ldots=\tilde{\beta}_{r}=0$, which is typically referred to as a sub-hypothesis (e.g. Adichie, 1978). For this reason, we will focus on (3) without loss of generality throughout the rest of the paper.

Our idea is to construct a pool of linear statistics $S=(S_{0},S_{1},\ldots,S_{m})$ such that $S$ is distributionally invariant under the left shifting operator $\pi_{L}$ under the null, in the sense that

[TABLE]

[TABLE]

Let $\mathrm{Id}$ denote the identity mapping. Then $\mathcal{G}=\{\mathrm{Id},\pi_{L},\ldots,\pi_{L}^{m}\}$ forms a group, which we refer to as the cyclic permutation group. We say a pool of statistics $S$ is invariant under the cyclic permutation group if $S$ satisfies (4). The following proposition describes the main property of statistics that are invariant under the cyclic permutation group.

Proposition 1**.**

Assume that $S=(S_{0},S_{1},\ldots,S_{m})$ is invariant under the cyclic permutation group. Let $R_{0}$ be the rank of $S_{0}$ in descending order, i.e. $R_{0}=|\{j\geq 0:S_{j}\geq S_{0}\}|$. Then

[TABLE]

where $\mathcal{L}$ denotes law, $\succeq$ denotes stochastic dominance, and $\mathrm{Unif}([0,1])$ denotes the uniform distribution on $[0,1]$. Furthermore, $R_{0}\sim\mathrm{Unif}([m+1])$ if $S$ has no ties with probability $1$.

Based on the p-value defined in (6), we can derive a test that rejects the null hypothesis if $p\leq\alpha$. We refer to this simple test as a marginal rank test. The following proposition shows that the marginal rank test is valid in finite samples and can be exact under mild conditions.

Proposition 2**.**

Suppose $S=(S_{0},S_{1},\ldots,S_{m})$ is invariant under the cyclic permutation group under $H_{0}$ and let the p-value be defined as in (6). Then $\mathbb{P}_{H_{0}}(p\leq\alpha)\leq\alpha$. If $1/\alpha$ is an integer, $m+1$ is divisible by $1/\alpha$, and $S$ has no ties almost surely, then $\mathbb{P}_{H_{0}}(p\leq\alpha)=\alpha$.

In practice, the reciprocals of commonly-used confidence levels (e.g. $0.1,0.05,0.01,0.005$) are integers. In these cases it is sufficient to set $m=1/\alpha-1$ to obtain an exact test.

The rank used in the marginal rank test only gives one-sided information and may not be suitable for two-sided tests. More concretely, $S_{0}$ may be significantly different from $S_{1},\ldots,S_{m}$ under the alternative but the sign of the difference may depend on the true parameters. An intuitive remedy is to apply the marginal rank test on the following modified statistics

[TABLE]

Whenever $S_{0}$ is significantly different from $S_{1},\ldots,S_{m}$, $\tilde{S}_{0}$ is significantly larger than $\tilde{S}_{1},\ldots,\tilde{S}_{m}$. The following proposition guarantees the validity of the transformation (7). In particular, the transformation in (7) satisfies the condition.

Proposition 3**.**

If $S=(S_{0},S_{1},\ldots,S_{m})$ is invariant under the cyclic permutation group,

[TABLE]

is invariant under the cyclic permutation group for every $g$ such that

[TABLE]

In this article, we consider linear statistics

[TABLE]

and apply the marginal rank test on $\tilde{S}_{0},\ldots,\tilde{S}_{m}$ defined in (7). Partition $X$ into $(X_{[r]}\,\,X_{[-r]})$ and $\beta$ into $(\beta_{[r]},\beta_{[-r]})$. The linear model (1) implies that

[TABLE]

In the next three subsections we will show how to construct the $\eta_{j}$s to guarantee the Type I error control and to enhance power. Surprisingly, the only distributional assumption on $\epsilon$ is exchangeability:

1. A1

The error vector $\epsilon$ has exchangeable components, i.e. for any permutation $\pi$ on $[n]$,

[TABLE]

2.2 Construction for Type I Error Control

Under $H_{0}$, (8) can be simplified as

[TABLE]

To ensure the distributional invariance of $\{y^{T}\eta_{0},\ldots,y^{T}\eta_{m}\}$ under the cyclic permutation group, it is sufficient to construct $\eta_{j}$s such that the deterministic parts are identical for all $j$ and the stochastic parts are invariant under the cyclic permutation group. To match the deterministic parts, we can simply set $X_{[-r]}^{T}\eta_{j}$ to be equal, as stated in the following condition.

1. C1

There exists $\gamma_{[-r]}\in\mathbb{R}^{p-r}$ such that

[TABLE]

To ensure the invariance of the stochastic parts, intuitively the $\eta_{j}$s should be left shifted transforms of each other. To be concrete, consider the case where $n=6$ and $m=2$. Then given any $\eta^{*}=(\eta^{*}_{1},\eta^{*}_{2},\eta^{*}_{3},\eta^{*}_{4},\eta^{*}_{5},\eta^{*}_{6})^{T}$, the following construction would imply the invariance of $\{\epsilon^{T}\eta_{j}:j=0,1,\ldots,m\}$ under the cyclic permutation group:

[TABLE]

To see this, note that

[TABLE]

By assumption A1,

[TABLE]

[TABLE]

Using the same argument we can show $(\epsilon^{T}\eta_{0},\epsilon^{T}\eta_{1},\epsilon^{T}\eta_{2})\stackrel{{\scriptstyle d}}{{=}}(\epsilon^{T}\eta_{2},\epsilon^{T}\eta_{0},\epsilon^{T}\eta_{1})$ and thus the invariance of $(\epsilon^{T}\eta_{0},\epsilon^{T}\eta_{1},\epsilon^{T}\eta_{2})$ under the cyclic permutation group.

In general, if $n$ is divisible by $m+1$ with $n=(m+1)t$, then we can construct $\eta_{j}$ as a left shifted transform of a vector $\eta^{*}$, i.e.

[TABLE]

where $\pi_{L}$ is the left shifting operator defined in (5). More generally, if $n=(m+1)t+s$ for some integers $t$ and $0\leq s\leq m$, we can leave the last $s$ components to be the same across the $\eta_{j}$s while shifting the first $(m+1)t$ entries as in (10), as stated in the following condition.

1. C2

There exists $\eta_{*}\in\mathbb{R}^{n}$ such that

[TABLE]

where $t=\lfloor n/(m+1)\rfloor$.

Proposition 4**.**

Under assumption A1, $(y^{T}\eta_{0},\ldots,y^{T}\eta_{m})$ is distributionally invariant under the cyclic permutation group if $(\eta_{0},\ldots,\eta_{m})$ satisfies C1 and C2.

Now we discuss the existence of $(\eta_{*},\gamma_{[-r]})$. Note that $\eta_{j}$ is a linear transformation of $\eta^{*}$. Let $I_{p-r}$ denote the identity matrix of size $p-r$ and $\Pi_{j}\in\mathbb{R}^{n\times n}$ be the matrix such that $\eta_{j}=\Pi_{j}\eta^{*}$. Then C1 and C2 imply that

[TABLE]

The above linear system has $(m+1)(p-r)$ equations and $n+p-r$ unknowns. Therefore, a non-zero solution always exists if $(m+1)(p-r)<n+p-r$.

Theorem 1**.**

Under assumption A1,

1. (a)

The linear system (11) always has a non-zero solution if

[TABLE] 2. (b)

for any solution $(\gamma_{[-r]},\eta^{*})$ of (11),

[TABLE]

is invariant under the cyclic permutation group under $H_{0}$, where $\Pi_{j}\in\mathbb{R}^{n\times n}$ is the coefficient matrix that maps $\eta^{*}$ to $\eta_{j}$ defined in C2.

Suppose $\alpha=0.05$ for illustration and set $m=1/\alpha-1=19$. Then the condition (12) reads

[TABLE]

Even when $r=1$, this is satisfied in many applications. On the other hand, when $r$ is large but $p-r$ is small, then (12) can still be satisfied even if $p>n$. This is in sharp contrast to regression F-tests and permutation F-tests that require fitting the full model and thus $p\leq n$. Furthermore, we emphasize that Theorem 1 allows arbitrary design matrices. This is fundamentally different from the asymptotically valid tests which always impose regularity conditions on $X$.

2.3 Construction for high power when $r=1$

To guarantee reasonable power, we need $y^{T}\eta_{0}$ to be significantly different from the other statistics under the alternative. In this subsection we focus on the case where $r=1$ to highlight the key idea. The general case with $r>1$ is discussed in Appendix C.

When $\beta_{1}\not=0$, (8) implies that

[TABLE]

where $W_{j}=\epsilon^{T}\eta_{j}+(\textbf{1}^{T}\eta_{*})\beta_{0}+(X_{[-1]}^{T}\eta_{*})^{T}\beta_{[-1]}$ and $(W_{1},\ldots,W_{m})$ is invariant under the cyclic permutation group by Theorem 1. To enhance power, it is desirable that $X_{1}^{T}\eta_{0}$ lies far from $\{X_{1}^{T}\eta_{1},\ldots,X_{1}^{T}\eta_{m}\}$. In particular, we impose the following condition on the $\eta_{j}$s:

1. C3

there exists $\gamma_{1},\delta\in R$, such that

[TABLE]

Putting C1, C2 and C3 together, we obtain the following linear system,

[TABLE]

where $e_{1,p(m+1)}$ is the first canonical basis in $\mathbb{R}^{p(m+1)}$ and

[TABLE]

This linear system has $(m+1)p$ equations and $n+p+1$ variables. Thus it always has a non-zero solution if

[TABLE]

When $\alpha=0.05$ and $m=19$, this condition is still reasonable in many problems.

The normalized gap $\delta/\|\eta\|$ can be regarded as a proxy of power. Write $\gamma$ for \left(\begin{array}[]{c}\gamma_{1}\\ \gamma_{[-1]}\end{array}\right). It is natural to consider the following optimization:

[TABLE]

This linear programming problem can be solved by fitting a linear regression and it permits a closed-form solution. Let $O^{*}(X)$ denote the optimal value of the objective function, i.e. the maximum achievable value of $\delta$ in this case. Here we use the symbol $O^{*}(X)$ instead of $\delta(X)$ to distinguish the role of the objective value and the variable $\delta$. In spite of these coinciding when $r=1$, they are distinct when $r>1$; see Appendix C for details.

Theorem 2**.**

Assume that $n\geq pm$. Let

[TABLE]

Partition $B(X)$ into $[B(X)_{1}\,\,B(X)_{[-1]}]$ where $B(X)_{1}$ is the first column of $B(X)$. Further let

[TABLE]

where $+$ denotes the Moore-Penrose generalized inverse. Then $O^{*}(X)=\|\tilde{\eta}\|_{2}$ and one global maximizer of (18) is given by

[TABLE]

Remark 1**.**

When $B(X)_{[-1]}$ has full column rank, $\tilde{\eta}$ is the residual vector obtained by regressing $B(X)_{1}$ on $B(X)_{[-1]}$ and $\|\tilde{\eta}\|_{2}^{2}$ is the residual sum of squares. Both quantities can be easily computed using standard software. If $B(X)_{[-1]}$ does not have full column rank, then $\tilde{\eta}$ is the residual from minimum-norm least squares solution obtained by regressing $B(X)_{1}$ on $B(X)_{[-1]}$, which is the limit of the ridge estimator with the penalty level tending to zero and is the limiting solution of standard gradient descent initialized at zero (e.g. Hastie et al., 2019).

2.4 Pre-ordering rows of design Matrix

Given any $X$, we can easily calculate the proxy of signal strength $O^{*}(X)$ by Theorem 2. However, the optimal value is not invariant to row permutations of $X$. That is, for any permutation matrix $\Pi\in\mathbb{R}^{n\times n}$, $O^{*}(X)\not=O^{*}(\Pi X)$ typically. Roughly speaking, this is because $\delta^{*}(X)$ involves the left shifting operator, which depends on the arrangement of the rows of $X$. Figure 1(a) illustrates the variability of $O^{*}(\Pi X)$ as a function of $\Pi$ for a fixed matrix with $8$ rows and $3$ columns, generated with i.i.d. Gaussian entries.

Notably, even in such regular cases variability is non-negligible. This motivates the following secondary combinatorial optimization problem:

[TABLE]

This is a non-linear travelling salesman problem. Note that we aim at finding a solution with a reasonably large objective value instead of finding the global maximum of (20), which is NP-hard. For this reason, we solve (20) by the genetic algorithm, which is generally efficient for moderate $n$ albeit without a worst-case convergence guarantee. In a nutshell, a genetic algorithm maintains a population of permutations, generates new permutations by two operations: crossover and mutation, and evolves the population via a mechanism called selection, based on the objective value. We refer the readers to Michalewicz (2013) for more details.

We compare the genetic algorithm, implemented in R package gaoptim, with a simple competing algorithm that randomly selects ordering and keeps the one yielding the largest objective value. We refer to this method as SS. Although this competitor is arguably too weak and more efficient algorithms may exist, our goal here is simply to illustrate the effectiveness of the genetic algorithm instead of to claim the superiority of the genetic algorithm. We compare the performance of the genetic algorithm and SS on three matrices with $n=1000$ and $p=20$ as realizations generated from random one-way ANOVA matrices with exactly one entry in each row at a unifromly random position, random matrices with i.i.d. standard normal entries and random matrices with i.i.d. standard Cauchy entries. The results are plotted in Figure 1(b) where the y-axis measures $O^{*}(\Pi X)$, scaled by the maximum achieved by the genetic algorithm and SS for visualization, and the x-axis measures the number of random samples each algorithm accesses. The population size is set to be $10$ for the genetic algorithm in all scenarios. It is clear that the genetic algorithm consistently improves the solution while SS gets trapped after a few iterations.

2.5 Implementation of the CPT

We summarize the implementation of the CPT below:

1. Step 1

Compute a desirable pre-ordering $\Pi_{0}$ for the combinatorial optimization problem

[TABLE]

where $O^{*}(\cdot)$ is defined in Theorem 2 when $r=1$ or in Theorem 3 when $r>1$; 2. Step 2

Replace $y$ and $X$ by $\Pi_{0}y$ and $\Pi_{0}X$; 3. Step 3

Compute $\eta^{*}$ via Theorem 2 when $r=1$ or via Theorem 3 when $r>1$; 4. Step 4

Compute $S_{j}=y^{T}\eta_{j}$ for $j=0,1,\ldots,m$ where

[TABLE] 5. Step 5

Compute \tilde{S}_{j}=\big{|}S_{j}-\mbox{median}\{(S_{j})_{j=0}^{m}\}\big{|}; 6. Step 6

Compute the p-value $p=R_{0}/(m+1)$ where $R_{0}$ is the rank of $\tilde{S}_{0}$ in the set $\{\tilde{S}_{0},\tilde{S}_{1},\ldots,\tilde{S}_{m}\}$ in descending order; 7. Step 7

Reject the null hypothesis if $p\leq\alpha$.

The computational cost of Step 3 is the same as solving a linear regression with the sample size $n$ and dimension $pm$, as indicated by Theorem 2. As a result, the computational cost of Step 2 and Step 4-7 are negligible. If the computing budget is tight, a random ordering  can be used for Step 1, for which the computational cost is negligible. Otherwise, a genetic algorithm can be used instead, of which the computational cost is the same as solving $M$ linear regressions of the same size as in Step 3, where $M$ is the total number of samples in the solution path. Admittedly, the latter option is computationally intensive compared to regression t- or F-tests and permutation tests – the former involves solving a single linear regression with a smaller dimension $p$ and the latter involves solving $M$ linear regressions of the same size, where $M$ is the number of permutations. However, for moderate-sized problems, the computational time of our method is acceptable. On the other hand, if the genetic algorithm is replaced by a more efficient search algortihm, the computational cost can be drastically reduced. We discuss one potential algorithm in Section 4.3.

3 Experiments

To assess the power of our procedure, we conduct extensive numerical experiments. In all the experiments we fix the sample size $n=1000$ and consider three values $25,33,40$ for dimension $p$ such that the sample per parameter $n/p\approx 40,30,25$. Given a value of $p$, we consider the three types of design matrices considered in Figure 1(b). For each type of design matrices, we generate $50$ independent copies. Given each $X$, we generate 3000 copies of $\epsilon$ with independent entries from the standard normal distribution and standard Cauchy distribution.

We consider two variants of the CPT, one with random ordering and one with pre-ordering by the genetic algorithm, as well as five competing tests: (i) the t- or F-test; (ii) the permutation t- or F-test which approximates the null distribution of the t- or F-statistic by the permutation distribution with $X_{[r]}$ reshuffled; (iii) the Freedman-Lane test (Freedman & Lane, 1983; Anderson & Robinson, 2001; Toulis, 2019) which approximates the null distribution of the t- or F-statistic by the permutation distribution with reduced-form regression residuals reshuffled; (iv) the asymptotic z-test for least absolute deviation (LAD) regression; (v) the GroupBound method (Meinshausen, 2015). For methods (ii) and (iii), we calculate the test based on 1000 random permutations. To further demonstrate the importance of the pre-ordering step in the CPT, we consider a weaker pre-ordering with 1000 random samples and a stronger pre-ordering with 10000 random samples for the genetic algorithm. All tests will be performed with level $\alpha=0.05$ and the number of statistics $m+1$ is set to be $20$ for the CPT. All programs to replicate the results in this article can be found in https://github.com/lihualei71/CPT.

Owing to the space constraint, here we only present the results for testing a single coordinate, i.e. $H_{0}:\beta_{1}=0$, while leaving other results to Appendix D. Since all tests considered here are invariant with respect to $\beta_{[-1]}$, we assume $\beta_{[-1]}=0$ without loss of generality. Given a design matrix $X$ and an error distribution $F$, we start by computing a benchmark signal-to-noise ratio $\beta_{1}$ such that the t- or F-test has approximately 20% power, using Monte-Carlo simulation, where $y$ is generated from

[TABLE]

Then all tests are performed on $X$ and the following $18000$ outcome vectors $y_{s}^{(b)}$, respectively:

[TABLE]

For each $s$, the proportion of rejections among 3000 $\epsilon$’s is computed. When $s=0$, this proportion serves as an approximation of the Type I error and should be closed to or below $\alpha$ for a valid test; when $s>0$, it serves as an approximation of power and should be large for a powerful test. For each of the three types of design matrices, the above experiments are repeated on 50 independent copies of $X$’s.

Figure 2 presents the Type I error of all tests for three types of design matrices. The boxplots present the variation among 50 independent copies of design matrices. In all cases, the three variants of the CPT are valid, as guaranteed by our theory, while GroupBound is overly conservative. The permutation test and Freedman-Lane test appear to be valid in our simulation settings even though there is no theoretical guarantee for heavy-tailed errors. When errors are Gaussian, the t-test is valid, as guaranteed by theory, but can be conservative or anti-conseravative with heavy-tailed errors depending on the design matrix. On the other hand, the test based on the LAD regression is anti-conservative when $X$ is a realization of Gaussian matrices and the error distribution is Gaussian or Cauchy, although validity can be proved asymptotically under regularity conditions that are satisfied by realizations of Gaussian matrices with high probability (e.g. Pollard, 1991). This makes a case for the fragility of some asymptotic guarantees.

For power comparison, we only show results for the case where the design matrices are realizations of Gaussian (resp. Cauchy) matrices and errors are Gaussian (resp. Cauchy) in Figure 3; the results for other cases will be presented in Appendix D. All figures plot the median power ratio, obtained from $50$ independent copies of $X$’s, between each variant of the CPT and each competing test. First we see that GroupBound has zero power in all scenarios, so the power ratios are infinite, and hence missing in the plots. Second, the pre-ordering step plays an important role in raising the power of the CPT. Third, the relative power of the CPT, with ordering via the genetic algorithm, improves as $n/p$ increases. Furthermore, in the Gaussian case, it is not surprising that the t-test is the most powerful one because it is provably the uniformly most powerful unbiased test for linear models with Gaussian errors. The efficiency loss of the CPT against the t-test, permutation t-test and the test based on LAD regression is moderate in general and is low when the sample size per parameter and the signal-to-noise ratio is large. In the Cauchy case, the CPT is more powerful than the t-test.

4 Discussion

4.1 Confidence interval/region by inverting the test

It is straightforward to deduce a confidence region for $\beta_{[r]}$ by inverting the CPT. Specifically, the inverted confidence region is given by $\mathcal{I}\triangleq\left\{\beta_{[r]}:p(y-X\beta;X)>\alpha\right\}$, where $p(y;X)$ is the p-value produced by the CPT with a design matrix $X$ and an outcome vector $y$. Under the construction C3,

[TABLE]

Thus,

[TABLE]

Then $\mathcal{I}$ can be simplified as

[TABLE]

where $x_{\min}$ and $x_{\max}$ are the infimum and the superimum of $x$ such that

[TABLE]

and $m(x)=\mbox{median}\left[\{y^{T}\eta_{j}-xI(j=0)\}_{j=0}^{m}\right]$. When $r=1$, the confidence interval (21) gives a useful confidence interval simply as

[TABLE]

where $x_{\min}$ and $x_{\max}$ are the smallest and the largest solutions of (22). When $r>1$, the confidence region (21) may not be useful because it is unbounded. More precisely, $\beta_{[r]}\in\mathcal{I}$ implies that $\beta_{[r]}+\xi\in\mathcal{I}$ for any $\xi$ orthogonal to $\delta$. We leave the construction of more efficient confidence regions to future research.

4.2 Connection to knockoff based inference

Our test is implicitly connected to the novel idea of knockoffs, proposed by Barber & Candès (2015) to control the false discovery rate for variable selection in linear models. Specifically, they assume a Gaussian linear model and aim at detecting a subset of variables that control the false discovery rate in finite samples. Unlike the single hypothesis testing considered in this paper, multiple inference requires dealing with the dependence between test statistics for each hypothesis carefully. They proposed an innovative idea of constructing a pseudo design matrix $\tilde{X}$ such that the joint distribution of $(X_{1}^{T}y,\ldots,X_{p}^{T}y,\tilde{X}_{1}^{T}y,\ldots,\tilde{X}_{p}^{T}y)$ is invariant to the pairwise swapping of $X_{j}^{T}y$ and $\tilde{X}_{j}^{T}y$ all for $j$ with $\beta_{j}=0$. Then the test statistic for testing $H_{0j}:\beta_{j}=0$ is constructed by comparing $X_{j}^{T}y$ and $\tilde{X}_{j}^{T}y$ in an appropriate way, thereby obtaining a valid binary p-value $p_{j}$ that is uniformly distributed on $\{1/2,1\}$ under $H_{0j}$. The knockoffs-induced p-values marginally resemble the construction of statistics in the CPT with $m=2,\eta_{0}=X_{j},\eta_{1}=\tilde{X}_{j}$. On the other hand, the validity of knockoffs essentially rests on the distributional invariance of $\epsilon$ under the rotation group while the validity of the CPT relies on the distributional invariance of $\epsilon$ under the cyclic permutation group. This coincidence illustrates the charm of group invariance in statistical inference.

4.3 More efficient algorithm for pre-ordering

Although a genetic algorithm is able to solve (20) efficiently for moderate-sized problems, it is not scalable enough to handle big data. Since the exact minimizer is not required, we can resort to other heuristic algorithms. One heuristic strategy is proposed by Fogel et al. (2013) by relaxing permutation matrice into doubly stochastic matrices, with $\Pi\textbf{1}=\Pi^{T}\textbf{1}=0$ and $\Pi_{ij}\geq 0$, and optimizing the objective using continuous optimization algorithms. This may suggest an efficient gradient based algorithm. We leave this as a future direction.

Acknowledgement

Peter J. Bickel was supported by the National Science Foundation (DMS 82978). The authors are grateful to Peng Ding, William Fithian, editors and reviewers for their constructive feedback.

Appendix A Technical Proofs

Proof of Proposition 1.

Let $R_{j}$ be the rank of $S_{j}$ in descending order as defined in (6). Then the invariance of $S$ implies the invariance of $(R_{0},R_{1},\ldots,R_{m})$. As a result,

[TABLE]

Then for any $k$,

[TABLE]

Let $S_{(1)}\geq S_{(2)}\geq\cdots\geq S_{(m+1)}$ be the ordered statistics of $(S_{0},\ldots,S_{m})$, which may involve ties. Then by definition, $R_{j}\geq k$ whenever $S_{j}\leq S_{(k)}$ and thus,

[TABLE]

implying that $\mathcal{L}(R_{0})\succeq\mathrm{Unif}([m+1])$. When there is no tie, the set $\{R_{0},R_{1},\ldots,R_{m}\}$ is always $\{1,2,\ldots,m+1\}$ and thus

[TABLE]

∎

Proof of Proposition 2.

By Proposition 1, $\mathcal{L}(p)\succeq\mathrm{Unif}([0,1])$. Thus $\mathbb{P}_{H_{0}}(p\leq\alpha)\leq\alpha$. If $S$ has no ties almost suresly and $1/\alpha=(m+1)/b$ for some integer $b$, Proposition 1 implies that $R_{0}\sim\mathrm{Unif}([m+1])$ and thus

[TABLE]

∎

Proof of Proposition 3.

By definition,

[TABLE]

where the last line uses the invariance of $g$. The proof is completed by noting that $S\stackrel{{\scriptstyle d}}{{=}}\pi_{L}(S)$. ∎

Proof of Proposition 4.

It is left to prove the invariance of $(\epsilon^{T}\eta_{0},\ldots,\epsilon^{T}\eta_{m})$ under the cyclic permutation group. Further, since the last $n-(m+1)t$ terms are the same for all $j$, it is left to prove the case where $n$ is divisible by $m+1$. Let $\tilde{\Pi}$ be the permutation matrix corresponding to $\pi_{L}^{t}$. Then C2 implies that

[TABLE]

Repeating (23) for $m-1$ times, we prove the invariance of $(\epsilon^{T}\eta_{1},\ldots,\epsilon^{T}\eta_{m})$ under the cyclic permutation group. ∎

Proof of Theorem 1.

When $n/(p-r)>m$, the number of variables $n+p-r$ of (11) is larger than the number of equations $(m+1)(p-r)$. Part (a) is then proved. For any solution of (11), by (9),

[TABLE]

The proof is completed by noting that the deterministic parts are identical for all $j$ and the stochastic parts are invariant under the cyclic permutation group by Proposition 4. ∎

Proof of Theorem 2.

First, (13) can be equivalently formulated as

[TABLE]

This can be further rewritten as

[TABLE]

For any $\eta$ satisfying the second constraint,

[TABLE]

and thus

[TABLE]

As a result,

[TABLE]

In other words, we have shown that $\delta^{*}(X)\leq\|\tilde{\eta}\|_{2}$. On the other hand, the vector $\tilde{\eta}/\|\tilde{\eta}\|_{2}$ satisfies the constraint (24) and

[TABLE]

This shows that $\delta^{*}(X)\geq\|\tilde{\eta}\|_{2}$. In this case, it is obvious that $O^{*}(X)=\delta^{*}(X)$. Therefore, $O^{*}(X)=\|\tilde{\eta}\|_{2}$ and one maximizer is $\eta^{*}(X)=\tilde{\eta}/\|\tilde{\eta}\|_{2}$. ∎

Appendix B 1908-2018: A Selective Review of The Century-Long Effort

The linear model is fundamental in the history of statistics and has been developed for over a century. Nowadays it is still a widely-used model for data analysts to demystify complex data, as well as a powerful tool for statisticians to understand complicated methods and expand the toolbox for advanced tasks. It is impossible to exhaust the literature for this long-standing problem. We thus provide a selective yet extensive review to highlight milestones in the past century. In particular, we focus on the linear hypothesis testing problem, as well as the estimation problem which can yield the former, for vanilla linear models with general covariates, and briefly discuss special cases such as location problems and ANOVA problems when necessary. However, we exclude the topics such as Bayesian linear models, high dimensional sparse linear models, selective inference for linear models, linear models with dependent errors, high breakdown regression methods, linear time series, and generalized linear models. We should emphasize that these topics are at least equally important as those discussed in this section; they are excluded simply to avoid digression.

B.1 Normal theory-based tests

Motivated by the seminal work by Student (1908b) and Student (1908a) which propose the one-sample and two-sample t-test, Ronald A. Fisher derived the well-known t-distribution (Fisher, 1915) and applied it to testing a single regression coefficient in homoskedastic Gaussian linear models (Fisher, 1922). In his 1922 paper, he also derived a test that is equivalent to the F-test for testing the global null under the same setting. Later he derived the F-distribution (Fisher, 1924) which he characterized through “z”, the half logarithm of F-statistics, and proposed the F-test for ANOVA problems. Both tests were elaborated in his impactful book (Fisher, 1925), and the term “F-test” was coined by George W. Snedecor (Snedecor, 1934).

This line of foundational work established the first generation of rigorous statistical tests for linear models. They are exact tests of linear hypotheses in linear models with i.i.d. normal errors and arbitrary fixed-design matrices. Although the exactness of the tests requires no assumption on the design matrices, the normality assumption can rarely be justified in practice. Early investigations of the test validity with non-normal errors can be dated back to Egon S. Pearson (Pearson, 1929; Pearson & Adyanthāya, 1929; Pearson, 1931). Unlike the large-sample theory that is standard nowadays, the early works took an approximation perspective to improve the validity in small samples. It was furthered in the next few decades (e.g. Eden & Yates, 1933; Bartlett, 1935; Geary, 1947; Gayen, 1949, 1950; David & Johnson, 1951b, a; Box, 1953; Box & Watson, 1962; Pearson & Please, 1975) and it was mostly agreed that the regression t-test is extremely robust to non-normal errors with a moderately large sample size (e.g. $>30$) while the regression F-test is more sensitive to the deviation from normality. It is worth emphasizing that these results were either based on mathematically unrigorous approximation or based on the Edgeworth expansion theory that could be justified rigorously (e.g. Esseen, 1945; Wallace, 1958; Bhattacharya & Ghosh, 1978) in the asymptotic regime that the sample size tends to infinity while the dimension of the parameters stays relatively low (e.g. a small constant).

Later on, due to the popularization of rigorous large-sample theories in 1950s (e.g. LeCam, 1953; Chernoff, 1956), pioneered by Doob (1935), Wilks (1938), Mann & Wald (1943), and Wald (1949), statisticians started to investigate the validity of regression t- and F-tests in certain asymptotic regimes. This can be dated back to Friedhelm Eicker (Eicker, 1963, 1967), to the best of our knowledge, and developed by Peter J. Huber in his well-known and influential paper (Huber, 1973), which showed that the least squares estimate is jointly asymptotically normal if and only if the maximum leverage score tends to zero. This clean and powerful result laid the foundation to asymptotic analyses for the t- and F-test (e.g. Arnold, 1980). Notably these early works did not assume that the dimension $p$ stays fixed, as opposed to the simplified arguments in standard textbooks. Before 1990s, the large-sample theory for least squares estimators were well established in the regime where the sample size per parameter $n/p$ grows to infinity, under regularity conditions on the design matrices and on the errors, typically with i.i.d. elements and finite moments. It shows that both the t- and F-test are asymptotically valid and can be approximated by the z- and $\chi^{2}$-test, respectively. For the t-test, the robustness to non-normality was proved even without typical regularity conditions (e.g. Zellner (1976); Jensen (1979) for spherically invariant errors, Efron (1969); Cressie (1980); Benjamini (1983); Pinelis (1994) for orthant symmetric errors) or beyond the aforementioned regime (e.g. Lei et al., 2018). In contrast, though similar results exist for the F-test (e.g. Zellner, 1976), more non-robustness results were established. For instance, a line of work (e.g. Boos & Brownie, 1995; Akritas & Arnold, 2000; Calhoun, 2011; Anatolyev, 2012) showed that the F-test is asymptotically invalid, unless the errors are normal, in the moderate dimensional regime where $n/p$ stays bounded as $n$ approaches infinity, although correction is available under much stronger assumptions on the design matrix or the coefficient vectors. Even with normal errors, Zhong & Chen (2011) showed that the power of the F-test diminishes as $n/p$ approaches $1$. In sum, there have been tremendous efforts over the past century put into the robustness of the regression t- and F-test and it was agreed that the t-test is insensitive to non-normality, high dimensions and irregularity of design matrices to certain extent while the F-test is less robust in general.

B.2 Permutation tests

Despite tremendous attentions on the regression t- and F-test, other methodologies were developed in parallel as well. The earliest alternative is the permutation test, which justifies the significance of the test through the so-called “permutation distribution”. However, the early attempts to justify permutation tests were based on the “randomization model” in contrast to the “population model” that we considered in (1). The “randomization model” was introduced by Jerzy S. Neyman in his master thesis (Neyman, 1923) and coined by Ronald A. Fisher in 1926 (Fisher, 1926). It is also known as the Neyman-Rubin model (Rubin, 1974), or design-based inference (Särndal et al. (1978), in contrast to the model-based inference), or “conditional-on-errors” model (Kennedy (1995), in contrast to the “conditional-on-treatment” model). The theoretical foundation of permutation tests was laid by Edwin J. G. Pitman in his three seminal papers (Pitman, 1937a, b, 1938), with the last two focusing on regression problems, albeit under the “randomization model”. The early works viewed the permutation test as a better machinery in terms of the logical coherence and robustness to non-normality (e.g. Geary, 1927; Eden & Yates, 1933; Fisher, 1935). They found that the permutation distribution under the “randomization model” mostly agree with the normality-based distribution under the “population model”, until 1937 when Li B. Welch disproved the agreement for Latin-squares designs (Welch, 1937). In the next few decades, most works on permutation tests were established under the “randomization model” without being justified under the “population model”. We will skip the discussion of this period and refer to Berry et al. (2013) for a thorough literature review on this line of work, because our paper focuses on the “population model” like (1).

The general theory of permutation tests under the “population model” can be dated back to Hoeffding (1952) and Box & Andersen (1955), and was further developed by e.g. Romano (1989), Romano (1990), Chung & Romano (2013). For regression problems, early studies investigated special cases in ANOVA problems (Mehra & Sen, 1969; Brown & Maritz, 1982; Welch, 1990). For testing a single regression coefficient, Oja (1987) and Collins (1987) proposed permutation tests on a linear statistic and the F-statistic by permuting the covariates of interest. Whereas the procedure can be easily validated for univariate regressions, the validity was only justified under the “randomization model” when $p>1$. Manly (1991) proposed permuting the response vector $y$, which is valid for testing the global null $\beta=0$ but not for general linear hypotheses. Freedman & Lane (1983), Ter Braak (1992) and Kennedy & Cade (1996) proposed three different permutation tests on regression residuals. The theoretical guarantees of the aforementioned tests were established in a later review paper by Anderson & Robinson (2001). The main take-away message is that the permutation test should be performed on asymptotically pivotal statistics. For instance, to test for a single coefficient, the permutation t-test is asymptotically valid. This was further confirmed and extended by DiCiccio & Romano (2017) to heteroscedastic linear models with random designs.

B.3 Rank-based tests

Rank-based methods for linear models can be dated back to 1936, when Hotelling & Pabst (1936) established the hypothesis testing theory for rank correlation, nowadays known as the Spearman’s correlation. This work can be regarded as an application of rank-based methods for univariate linear models. Appealed by the normality-free nature of rank-based tests, Milton Friedman extended the idea to one-way ANOVA problems (Friedman, 1937). It can be identified as the first application of rank-based methods for multivariate linear models and was further developed by Kendall & Smith (1939) and Friedman (1940). Friedman’s test transforms continuous or ordinal outcomes into ranks. It was widely studied for ANOVA problems, started by the famous Kruskal-Wallis test for one-way ANOVA (Kruskal & Wallis, 1952) and extended to two-way ANOVA problems and factorial designs (Hodges & Lehmann, 1962; Puri & Sen, 1966; Sen, 1968b; Conover & Iman, 1976, 1981; Akritas, 1990; Akritas & Arnold, 1994; Brunner & Denker, 1994; Akritas et al., 1997). As of 90s, motivated by the advances of high dimensional asymptotic theories, further progresses have been made to refine the procedures in presence of large number of factors or treatments (Brownie & Boos, 1994; Boos & Brownie, 1995; Wang & Akritas, 2004; Bathke & Lankowski, 2005; Bathke & Harrar, 2008).

However the aforementioned works are restricted to ANOVA problems, with a few exceptions (e.g. Sen, 1968a, 1969), and fundamentally different from the modern rank tests based on regression R-estimates, themselves based on ranks of regression residuals. The first R-estimate-based test can be dated back to Hájek (1962), which derived the asymptotically most powerful rank test for univariate regressions when the error distribution is known. Adichie (1967a) extended the idea to testing the intercept and the regression coefficient simultaneously. It was further extended to testing the global null for multivariate regressions (Koul, 1969). Tests for general sub-hypotheses were first proposed by Koul (1970) and Puri & Sen (1973) for bivariate regressions. The general theory of testing sub-hypotheses were independently developed by Srivastava (1972), McKean & Hettmansperger (1976) and Adichie (1978). The underlying theory is based on the seminal work by Jana Jureckova (Jureckova, 1969), as a significant generalization of Hodges & Lehmann (1963) for location problems and Adichie (1967b) for univariate regressions. Her work was further extended by Jureckova (1971) and van Eeden (1972). However, these approaches are computationally intensive due to the discreteness of ranks. A one-step estimator was proposed by Kraft & Van Eeden (1972), which is asymptotically equivalent to the maximum likelihood estimators if the error distribution is known. Another one-step rank-based estimator, motivated by Bickel (1975) for M-estimators, was proposed by McKean & Hettmansperger (1978). On the other hand, Jaeckel (1972) proposed a rank-based objective function, later known as the Jaeckel’s dispersion function, that is convex in $\beta$ whose minimizer is asymptotically equivalent to Jureckova’s score-based estimators. Hettmansperger & McKean (1978) found an equivalent but mathematically more tractable formulation of the Jaeckel’s dispersion function as the sum of pairwise differences of regression residuals. A weighted generalization of the dispersion function was introduced by Sievers (1983), which unifies the Jaeckel’s dispersion function and Kendall’s tau-based dispersion function (Sen, 1968a; Sievers, 1978). Three nice survey papers were written by Adichie (1984), Aubuchon & Hettmansperger (1984), and Draper (1988). In 90s, motivated by the development of quantile regressions (Koenker & Bassett, 1978), Gutenbrunner & Jureckova (1992) found an important coincidence between the dual problem of the quantile regression and the “rank-score process”, which generalizes the notion introduced by Hajek & Sidak (1967) to linear models. Gutenbrunner et al. (1993) then developed a rank-score test for linear hypotheses; see also Koenker (1997) for a review. In the past two decades, there were much fewer works on rank-based tests for linear models (e.g. Feng et al., 2013).

B.4 Tests based on regression M-estimates

Regression M-estimates were introduced by Peter J. Huber in 1964 for location problems (Huber, 1964). The idea was soon extended to linear models by Relles (1968), who proved the asymptotic theory for Huber’s loss with $p$ fixed and $n$ tending to infinity. The theory was further extended to general convex loss functions by Yohai (1972). Despite the appealing statistical properties, the computation remained challenging in 1970s. Bickel (1975) proposed one-step M-estimates that are computational tractable with the same asymptotic property as full M-estimates. In addition, he proved the uniform asymptotic linearity of M-estimates, which is a fundamental theoretical result that laid the foundation for later works. Based on Bickel (1975)’s technique, Jureckova (1977) established the relation between regression M- and R-estimates. The asymptotic normality of M-estimates directly yields an asymptotically valid Wald-type test for general linear hypotheses. Schrader & Hettmansperger (1980) developed an analogue of the likelihood-ratio test based on M-estimators for sub-hypotheses. It was further extended to general linear hypotheses by Silvapulle (1992). However, both Wald-type tests and likelihood-ratio-type tests involve unknown nuisance parameters. To get rid of them, Sen (1982) proposed the M-test as an analogue of the studentized score test, which is able to test general linear hypotheses with merely an estimate of regression coefficients under the null hypothesis. It is known that the Rao’s score test may not be efficient in presence of nuisance parameters. Singer & Sen (1985) discussed an efficient test, which is essentially the analogue of Neyman’s $C(\alpha)$ test based on projected scores (Neyman, 1959), although it brings back nuisance parameters. M-tests were later investigated and generalized in a general framework based on influence functions (e.g. Boos, 1992; Markatou & Ronchetti, 1997).

As with the regression t- and F-test, the robustness to high dimensionality was investigated extensively for M-estimators in general linear models. In Huber’s 1972 Wald Lectures (Huber, 1972), he conjectured that the asymptotic normality of M-estimates proved by Relles (1968) can be extended to the asymptotic regime where $p$ grows with $n$. The conjecture was proved one year later in the regime $\kappa p^{2}=o(1)$, where $\kappa$ is the maximum leverage score, which implies $p=o(n^{1/3})$ (Huber, 1973). This was improved to $\kappa p^{3/2}=o(1)$ by Yohai & Maronna (1979), which implies that $p=o(n^{2/5})$, to $p=o(n^{2/3}/\log n)$ by Portnoy (1985) under further regularity conditions on the design matrix, and to $\kappa n^{1/3}(\log n)^{2/3}=o(1)$, which implies that $p=o(n^{2/3}/(\log n)^{2/3})$. All aforementioned results are derived for smooth loss functions. For non-smooth loss functions, Welsh (1989) obtained the first asymptotic result in the regime $p=o(n^{1/3}/(\log n)^{2/3})$. It was improved to $p=o(n^{1/2})$ by Bai & Wu (1994). For a single coordinate, Bai & Wu (1994) showed the asymptotic normality in the regime $p=o(n^{2/3})$. These works prove that the classical asymptotic theory holds if $p<\!\!<n^{2/3}$. However, in moderate dimensions where $p$ grows linear with $n$, the M-estimates are no longer consistent in $L_{2}$ metric. For certain random designs, the estimation error $\|\hat{\beta}-\beta\|_{2}^{2}$ converges to a non-vanishing quantity determined by $p/n$, the loss function and the error distribution through a complicated system of non-linear equations (El Karoui et al., 2011; Bean et al., 2012; El Karoui, 2013; Donoho & Montanari, 2016; El Karoui, 2018). This surprising phenomenon marks the failure of the classical asymptotic theory for M-estimators. For least squares estimators, Lei et al. (2018) showed that the classical t-test with appropriate studentization is still asymptotically valid under regularity conditions on the design matrix. Cattaneo et al. (2018) proposed a refined test for heteroscedastic linear models. However it is unclear how to test general linear hypotheses with general M-estimators in this regime, even for a single coordinate. Lei et al. (2018) provided the only fixed-design result for the asymptotic property of a single coordinate of general M-estimates in this regime. For the purpose of hypothesis testing, the null variance needs to be estimated but no consistent variance estimator is known at this moment, except for special random designs (e.g. Bean et al., 2012).

B.5 Tests based on regression L-estimates

L-estimators constitute an important class of robust statistics based on linear combination of order statistics. Frederick Mosteller proposed the first L-estimator for Gaussian samples (Mosteller, 1946). This was further developed in the following two decades (e.g. Hastings et al., 1947; Lloyd, 1952; Evans & Evans, 1955; Jung, 1956; Tukey, 1960; Bickel, 1965; Gastwirth, 1966). In particular, John W. Tukey advocated the trimmed mean and Winsorized mean in his far-reaching paper (Tukey, 1962), which he attributed to Charles P. Winsor based on their personal communication in 1941. One year later, the well-known Hodges-Lehmann estimator was developed (Hodges & Lehmann, 1963), which established the first connection between R- and L-estimates. For location problems, Bickel & Lehmann (1975) found the superiority of L-estimates over M- and R-estimates.

Despite the simplicity and the nice theoretical property of L-statistics, they are not easy to be generalized to linear models. The first attempt was made by Bickel (1973), which proposed a one-step L-estimate for general linear models. However, this estimator is not equivariant to affine transformations of the design matrix. Motivated by this paper, Welsh (1987) proposed a class of one-step L-estimators that are equivariant to reparametrization of the design matrix. Welsh (1991) further extended the idea to construct an adaptive L-estimator. Another line of thoughts were motivated by the pinoneering work of Koenker & Bassett (1978), which introduced the notion of regression quantiles as a natural analogue of sample quantiles for linear models. Although the quantile regression yields an M-estimator, it had been the driving force for the development of regression L-estimators since 80s. In this paper, they proposed another class of L-estimators by a discrete weighted average of regression quantiles and derived its asymptotic distribution. This idea was furthered by Koenker & Portnoy (1987) to L-estimators with continuous weights, by Portnoy & Koenker (1989) to adaptive L-estimators, and by Koenker & Zhao (1994) to heteroscedastic linear models. Another notable strategy of contructing L-statistics is based on weighted least squares with “outliers” removed. Ruppert & Carroll (1980) developed two equivariant one-step estimators as analogues of the trimmed mean. Both estimators can be formulated in the form of weighted least squares with extreme residuals removed. As with Ruppert & Carroll (1980), Jureckova (1983) proposed an analogue of the winsorized mean. The Bahadur representation of the trimmed mean least squares estimator was derived by Jureckova (1984). A nice review article of regression L-estimators was written by Alimoradi & Saleh (1998). The asymptotic results of L-estimators induce asymptotically valid Wald-type tests with consistent estimates of the asymptotic variance. Unlike M-estimators, we are not aware of other types of tests based on L-estimates.

B.6 Resampling-based tests

Resampling, marked by the jackknife (Quenouille, 1949, 1956; Tukey, 1958) and bootstrap (Efron, 1979), is a generic technique to assess the uncertainty of an estimator. Although both involving resampling, resampling-based tests are fundamentally different from permutation tests. The former approximates the sampling distribution under the truth while the latter approximates the sampling distribution under the null hypothesis, though they are asymptotically equivalent in many cases (e.g. Romano, 1989). Miller (1974) proposed the first jackknife-based estimator for general linear models. He showed that the estimator is asymptotically normal, the jackknife variance estimator is consistent, and thus the Wald-type test is asymptotically valid. Hinkley (1977) pointed out that Miller’s estimator is less efficient than the least squares estimator and proposed a weighted jackknife estimator to achieve efficiency. Wu (1986) proposed a general class of delete-$d$ jackknife estimators for estimating the covariance matrix of the least squares estimator. This was extended by Shao & Wu (1987), Shao (1988), Shao (1989), Peddada & Patwardhan (1992), and Liu & Singh (1992).

On the other hand, David A. Freedman first studied the bootstrapping procedures for linear models (Freedman, 1981). He studied two types of bootstrap: the residual bootstrap, where the regression residuals are resampled and added back to the fitted values, and the pair bootstrap, where the outcome and the covariates are resampled together. In the fixed-$p$ regime, he showed the consistency of the residual bootstrap under homoscedastic linear models and that of the pair bootstrap under general “correlation models” including heteroscedastic linear models. Navidi (1989), Hall (1989) and Qumsiyeh (1994) established the higher order accuracy of the pair bootstrap for linear models and the results were then presented under a broader framework in the influential monograph by Peter Hall (Hall, 1992). Wu (1986) found that the residual bootstrap fails in heteroscedastic linear models because its sampling process is essentially homoscedastic. To overcome this, he introduced another type of bootstrapping method based on random rescalings of regression residuals that match the first and second moment. Liu (1988) introduced a further requirement to match the third moment and improved the rate of convergence. Later Mammen (1993) coined this procedure the “wild bootstrap” and proved the consistency for least squares estimators under random-design homoscedastic and heteroscedastic linear models. Hu & Zidek (1995) proposed an alternative bootstrap procedure for heteroscedastic linear models that resamples the score function instead of the residuals. A wild bootstrap analogue of the score-based bootstrap was proposed by Kline & Santos (2012). In particular, they developed the bootstrap Wald tests and the boostrap score tests for general linear hypotheses.

The bootstrap techniques were also widely studied for regression M-estimates. The residual bootstrap was extended to M-estimators with smooth loss functions by Shorack (1982). Unlike the least squares estimator, it requires a debiasing step to obtain distributional consistency. Lahiri (1992) proposed a weighted residual bootstrap that does not require debiasing. He additionally showed the higher order accuracy of the weighted bootstrap and Shorack’s bootstrap for studentized M-estimators. However, this weighted bootstrap is hard to be implemented in general. On the other hand, motivated by Bayesian bootstrap (Rubin, 1981), Rao & Zhao (1992) proposed a bootstrapping procedure by randomly reweighting the objective function. This idea was extended by Chatterjee (1999) in a broader framework called “generalized bootstrap”. It was later reinvented by Jin et al. (2001) and referred to as “perturbation bootstrap”. The higher order accuracy of perturbation bootstrap was established by Das & Lahiri (2019). It was pointed out by Das & Lahiri (2019) that the perturbation bootstrap coincides with the wild bootstrap for least squares estimators. Hu & Kalbfleisch (2000) proposed another estimating function based bootstrap, as essentially a resampling version of Sen (1982)’s M-tests. The wild bootstrap was introduced for quantile regressions by Feng et al. (2011).

The robustness of bootstrap methods against high dimensions was widely studied in the literature. Bickel & Freedman (1983) proved the distributional consistency of the residual bootstrap for least squares estimators in the regime $p=o(n)$ for linear contrasts of $\beta$ and in the regime $p=o(n^{1/2})$ for the vector $\beta$, under fixed-design linear models with vanishing maximum leverage scores. They further showed the failure of bootstrap in moderate dimensions where $p/n\rightarrow c\in(0,1)$ and the usual variance rescaling does not help because the bootstrap distribution is no longer asymptotically normal. For M-estimators, Shorack (1982) showed that the debiased residual bootstrap is distributionally consistent in the regime $p=o(n^{1/3})$ for linear contrasts of $\beta$. The results were extended by Mammen (1989) to the regime $p=o(n^{2/3}/(\log n)^{2/3})$ for linear contrasts of $\beta$, and to the regime $p=o(n^{1/2})$ for the vector $\beta$. For random designs with i.i.d. observations, Mammen (1993) proved the distributional consistency of both the pair bootstrap and wild bootstrap for linear contrasts of $\beta$ in the regime $p=o(n^{a})$ for arbitrary $a<1$. He also proved the consistency under heteroscedastic linear models in the regime $p=o(n^{3/4})$ for the pair bootstrap, and in the regime $p=o(n^{1/2})$ for the wild bootstrap. This was further extended by Chatterjee (1999) to the generalized bootstrap, including the perturbation bootstrap (Rao & Zhao, 1992), $m$-out-of-$n$ bootstrap (Bickel & Sakov, 2008) and delete-d jackknife (Wu, 1990). On the other hand, extending Bickel & Freedman (1983)’s negative result, El Karoui & Purdom (2018) showed the failure of various bootstrap procedures for M-estimators in moderate dimensions, including the pair bootstrap, residual bootstrap, wild bootstrap and jackknife.

B.7 Other tests

A generic strategy for hypothesis testing is through pivotal statistics. Specifically, if there exists a statistic $S$ whose distribution is fully known, then the rejection rule $S\in\mathcal{R}^{c}$ for any region $\mathcal{R}$ with $P(S\in\mathcal{R})\geq 1-\alpha$ yields an exact test. For linear models, it is extremely hard to find a pivotal statistic under general linear hypotheses, except for Gaussian linear models under which the t- and F-statistics are pivotal. However, if the goal is to test all coefficients plus the intercept, i.e. $H_{0}:\beta_{0}=\gamma_{0},\beta=\gamma$, then one can recover the stochastic errors as $\epsilon_{i}=y_{i}-\gamma_{0}-x_{i}^{T}\gamma$ under the null and construct pivotal statistics based on $\epsilon$. Taking one step further, given a pivotal statistic, one can invert the above test to obtain a finite-sample valid confidence region $\mathcal{I}$ for $(\beta_{0},\beta)$, by collecting all $(\gamma_{0},\gamma)$’s at which the corresponding null hypothesis fails to be rejected. This induces a confidence region for $R^{T}\beta$ as $\mathcal{I}^{\prime}=\{R^{T}\beta:(\beta_{0},\beta)\in\mathcal{I}\}$. Finally, using the duality between the confidence interval and hypothesis testing again, the test which rejects the null hypothesis is valid for the linear hypothesis $H_{0}:R^{T}\beta=0$ in finite samples. If $r<\!\!<p$, this seemingly “omnibus test” is in general conservative and inferior to the tests discussed in previous subsections. Nonetheless, it stimulates several non-standard but interesting tests that are worth discussions.

The most popular strategy to construct pivotal statistics is based on quantiles of $\epsilon_{i}$s, especially the median. Assuming $\epsilon_{i}$s have zero median, Fisher (1925) first introduced the sign test for location problems, which was investigated and formalized later by Cochran (1937). Thirteen years later, Henri Theil proposed an estimator for univariate linear models (Theil, 1950a, b, c), later known as the Theil-Sen estimator (Sen, 1968a). Brown & Mood (1951) proposed a median test for general linear models by reducing the problem into a contingency table and applying the $\chi^{2}$-tests. The theoretical property of the Brown-Mood test was studied by Kildea (1981) and Johnstone & Velleman (1985). Daniels (1954) proposed a geometry-based test for univariate linear models, which can be regarded as a generalization of the Brown-Mood test. It was later connected to the notion of regression depth (Rousseeuw & Hubert, 1999) and applied in deepest regression methods (Van Aelst et al., 2002). The idea of inverting the sign test was exploited in Quade (1979). An analogue incorporating Kendall’s tau between the residuals and covariates was proposed by Lancaster & Quade (1985). The idea also attracted some attention in signal processing (e.g. Campi & Weyer, 2005; Campi et al., 2009) and econometrics (e.g. Chernozhukov et al., 2009). It should be noted that the approach is computationally infeasible even in low dimensions. Assuming further the symmetry of $\epsilon_{i}$s, Hartigan (1970) proposed a non-standard test based on an interesting notion of typical values. It was designed for location problems but can be applied to certain ANOVA problems. Furthermore, Siegel (1982) proposed the repeated median estimator and Rousseeuw (1984) proposed the least median squares estimators to achieve a high breakdown point.

The pivotal statistics can also be constructed in other ways. Parzen et al. (1994) proposed a bootstrap procedure based on inverting a pivotal estimating function at a random point. This procedure mimics the Fisher’s fiducial inference but can be justified under the frequentist framework. Recently Meinshausen (2015) proposed the GroupBound test for sub-hypotheses, which even works for high-dimensional settings where $p>\!\!>n$. However, the validity is only guaranteed for rotationally invariant errors with a known noise level. This assumption is extremely strong as shown by Maxwell (1860): a rotationally invariant random vector with i.i.d. coordinates must be multivariate Gaussian.

Appendix C Construction of $\eta$’s When $r>1$

Similar to C3, we impose the following restriction on $\eta$.

1. C3’

there exists $\gamma_{[r]},\delta\in\mathbb{R}^{r}$, such that

[TABLE]

Combining with (11), we obtain an analogue of (13) as follows.

[TABLE]

where $A(X)$ is defined in (14) and \gamma=\left(\begin{array}[]{c}\gamma_{[r]}\\ \gamma_{[-r]}\end{array}\right). This linear system involves $p(m+1)$ equations and $n+p+r$ variables. Therefore it always has a non-zero solution if

[TABLE]

Unlike the univariate case, there are infinite ways to characterize the signal strength since $\delta$ is multivariate. A sensible class of criteria is to maximize a quadratic form

[TABLE]

The following theorem gives the optimal solution given any weighting matrix $M$. Let $O^{*}(X)$ denote the optimal value of the objective function.

Theorem 3**.**

Assume that $n\geq pm-r+1$. Let $B(X)$ be defined in (19) in the main text. Partition $B(X)$ into $(B(X)_{[r]}\,\,B(X)_{[-r]})$ where $B(X)_{[r]}$ is the matrix formed by the first $r$ columns of $B(X)$. Let

[TABLE]

where

[TABLE]

Further let $\lambda_{\max}(M_{r}(X))$ denote the maximum eigenvalue, $u$ denote any eigenvector corresponding to it and $\tilde{\eta}=(I-H_{[-r]})u$. Then $O^{*}(X)=\lambda_{\max}(M_{r}(X))$ and

[TABLE]

is an optimal solution of (26).

Proof of Theorem 3.

Similar to the proof of Theorem 2, we first rewrite (25) as

[TABLE]

As a result, $\eta$ lies in the row null space of $B(X)_{[-r]}$ and thus

[TABLE]

Then

[TABLE]

Since $\|\eta\|_{2}\leq 1$,

[TABLE]

On the other hand, for any eigenvector $u$ of $M_{r}(X)$ corresponding to its largest eigenvalue, let $\tilde{\eta}=(I-H_{[-r]})u$ and $\eta=\tilde{\eta}/\|\tilde{\eta}\|_{2}$, then

[TABLE]

Thus, $\eta^{*}(X)=\tilde{\eta}/\|\tilde{\eta}\|_{2}$ is an optimal solution. As a result, $\delta^{*}(X)=B(X)_{[r]}^{T}\eta^{*}(X)$ and $O^{*}(X)=\lambda_{\max}(M_{r}(X))$. ∎

Although Theorem 3 gives the solution of (26) for arbitrary weight matrix $M$, it is not clear which $M$ is the best choice. Note that

[TABLE]

where $\tilde{W}_{j}=\gamma^{T}\beta+\eta_{j}^{T}\epsilon$ is invariant under the cyclic permutation group. Thus, $\delta^{T}\beta_{[r]}$ characterizes the signal strength. In principle, the “optimal” weight matrix should depend on the prior knowledge of $\beta_{[r]}$. For instance, for a Bayesian hypothesis testing problem with a prior distribution $Q$ on $\beta_{[r]}$ under the alternative, the optimal weight matrix is $M=\mathbb{E}_{Q}\left[\beta_{[r]}\beta_{[r]}^{T}\right]$.

Appendix D Complementary Experimental Results

D.1 Testing for a single coordinate

In this appendix we present experimental results that complement Section 3. Figure 4 - 7 present the power comparison for testing a single coordinate under the same setting as considered in Section 3 for four extra scenarios with realizations of Gaussian matrices + Cauchy errors, realizations of Cauchy matrices + Gaussian errors and realizations of random one-way ANOVA matrices + Gaussian or Cauchy errors, respectively.

D.2 Testing for multiple coordinates

Next we consider testing the first five coordinates $H_{0}:\beta_{1}=\ldots=\beta_{5}=0$, with a Bayesian alternative hypothesis

[TABLE]

All other settings are exactly the same as Section 3, except that the t-test and permutation t-test are replaced by the F-test and permutation F-test. For the CPT, we choose the weight matrix $M=\mathbb{E}[\beta_{[5]}\beta_{[5]}^{T}]$. Figure 8 presents the Monte-Carlo Type I error of all tests. The results are qualitatively the same as those in Section 3, though the F-test and LAD-based test become more invalid. Figure 9 - 14 present the power results under the same setting as Section D.2 for six scenarios with realizations of Gaussian matrices + Gaussian or Cauchy errors, realizations of Cauchy matrices + Gaussian or Cauchy errors and realizations of random one-way ANOVA matrices + Gaussian or Cauchy errors, respectively.