 # High-dimensional general linear hypothesis tests via non-linear spectral shrinkage

We are interested in testing general linear hypotheses in a high-dimensional multivariate linear regression model. The framework includes many well-studied problems such as two-sample tests for equality of population means, MANOVA and others as special cases. A family of rotation-invariant tests is proposed that involves a flexible spectral shrinkage scheme applied to the sample error covariance matrix. The asymptotic normality of the test statistic under the null hypothesis is derived in the setting where dimensionality is comparable to sample sizes, assuming the existence of certain moments for the observations. The asymptotic power of the proposed test is studied under various local alternatives. The power characteristics are then utilized to propose a data-driven selection of the spectral shrinkage function. As an illustration of the general theory, we construct a family of tests involving ridge-type regularization and suggest possible extensions to more complex regularizers. A simulation study is carried out to examine the numerical performance of the proposed tests.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In multivariate analysis, one of the fundamental inferential problems is to test a hypothesis involving a linear transformation of regression coefficients under a linear model. Suppose

is a matrix of observations modeled as

 Y =BX+Σ1/2pZ , (1.1)

where  (i)  is a matrix of regression coefficients;  (ii)  is a design matrix of rank ;  (iii)  is a

matrix with i.i.d. entries having zero mean and unit variance; and  (iv)

, a nonnegative definite matrix, is the population covariance matrix of the errors, with a “square-root” of so that . General linear hypotheses involving the linear model (1.1) are of the form

 H0:BC=0vs.Ha:BC≠0, (1.2)

for an arbitrary “constraints matrix” , subject to the requirement that

is estimable. Without loss of generality,

is taken to be of rank . Throughout, we assume that and are fixed, even as observation dimension and sample size increase to infinity. Henceforth,

is used to denote the effective sample size, which is also the degree of freedom associated with the sample error covariance matrix.

With various choices of and , the testing formulation incorporates many hypotheses of interest. For example, multivariate analysis of variance (MANOVA) is a special case. When the sample size is substantially larger than the dimension of the observations, this problem is well-studied. Anderson (1958) and Muirhead (2009) are among standard references. Various classical inferential procedures involve the matrices

 ˆΣp= 1nY(I−XT(XXT)−1X)YT, (1.3) ˆHp= 1nYXT(XXT)−1C[CT(XXT)−1C]−1CT(XXT)−1XYT, (1.4)

so that is the residual covariance of the full model, an estimator of , while is the hypothesis sums of squares and cross products matrix, scaled by . In a one-way MANOVA set-up, and are, respectively, the within-group and between-group sums of squares and products matrices, scaled by . In the rest of the paper, we shall refer to as the sample covariance matrix.

The testing problem (1.2) is well-studied in the classical multivariate analysis literature. Three standard test procedures are the likelihood ratio test (LR), Lawley–Hotelling trace test (LH) and Bartlett–Nanda–Pillai trace (BNP) test. They are called invariant tests, since under Gaussianity the null distributions of the test statistics are invariant with respect to . One common feature is that all test statistics are linear functionals of the spectrum of . Since this matrix is asymmetric, for convenience, a standard transformation is applied, giving the expressions of the invariant tests as follows. Define

 Qn =XT(XXT)−1C[CT(XXT)−1C]−1/2, (1.5) M0 =1nQTnYTˆΣ−1pYQn.

The matrix

is the “hat matrix” of the reduced model under the null hypothesis. Note that the non-zero eigenvalues of

are the same as those of . The test statistics for the LR, LH and BNP tests can be expressed as

 TLR0=∑qi=1log{1+λi(M0)}, TLH0=∑qi=1λi(M0), TBNP0=∑qi=1λi(M0)/{1+λi(M0)}.

The symbol denotes the -th largest eigenvalue of a symmetric matrix, further using the convention that and indicate the largest and smallest eigenvalue.

In contemporary statistical research and applications, high-dimensional data whose dimension is at least comparable to the sample size is ubiquitous. In this paper, focus is on the interesting boundary case when dimension and sample sizes are comparable. Primarily due to inconsistency of conventional estimators of model parameters — such as

—, classical test procedures for the hypothesis (1.2) — such as the LR, LH and BNP tests — perform poorly in such settings. When the dimension is larger than the degree of freedom , the invariant tests are not even well-defined because is singular. Even when is strictly less than , but the ratio is close to , these tests are known to have poor power behavior. Asymptotic results when were obtained in Fujikoshi et al. (2004) under Gaussianity of the populations, and more recently in Bai et al. (2017) under more general settings that only require the existence of certain moments.

Pioneering work on modifying the classical solutions in high dimension is in Bai et al. (2013), who corrected the scaling of the LR statistic when but , and are proportional to . The corrected LR statistic was shown to have significantly more power than its classical counterpart. In contrast, in this paper, we focus on the setting where and are fixed even as so that . In the multivariate regression problem, this corresponds to a situation where the response is high-dimensional, while the predictor is finite-dimensional. In the MANOVA problem, this framework corresponds to high-dimensional observations belonging to one of a finite number of populations.

To the best of our knowledge, when , the linear hypothesis testing problem has been studied in depth only for specific submodels of (1.1), primarily for the important case of two-sample tests for equality of population means. For the latter tests, a widely used idea is to construct modified statistics based on replacing with an appropriate substitute. This approach was pioneered in Bai and Saranadasa (1996) and further developed in Chen and Qin (2010). Various extensions to one-way MANOVA (Srivastava and Fujikoshi, 2006; Yamada and Himeno, 2015; Srivastava and Fujikoshi, 2006; Hu et al., 2017)

and a general multi-sample Behrens–Fisher problem under heteroscedasticity

(Zhou et al., 2017) exist. Other notable works for the two-sample problem include Biswas and Ghosh (2014); Chang et al. (2017); Chen et al. (2014); Guo and Chen (2016); Lopes et al. (2011); Srivastava et al. (2016); Wang et al. (2015). A second approach aims to regularize to address the issue of its near-singularity in high dimensions; see Chen et al. (2011) and Li et al. (2016) for ridge-type penalties in two-sample settings. Finally, another alternative line of attack consists of exploiting sparsity; see Cai et al. (2014); Cai and Xia (2014).

In this paper, we seek to regularize the spectrum of by flexible shrinkage functions. For a symmetric matrix and a function on , define

where

is the matrix of eigenvectors associated with the ordered eigenvalues of

. Now, consider any real-valued function on that is analytic over a specific domain associated with the limiting behavior of the eigenvalues of , as elaborated in Section 2. The proposed statistics are functionals of eigenvalues of the regularized quadratic forms

 M(f)=1nQTnYTf(ˆΣp)YQn.

Specifically, we propose regularized versions of LR, LH and BNP test criteria, namely

 TLR(f) =∑qi=1log{1+λi(M(f))}, TLH(f) =∑qi=1λi(M(f)), TBNP(f) =∑qi=1λi(M(f))/{1+λi(M(f))}.

These test statistics are designed to capture possible departures from the null hypothesis, when is replaced by , while suitable choices of the regularizer allow for getting around the problem of singularity or near-singularity when is comparable to .

Notice that has the same non-zero eigenvalues as . Thus, the proposed test family is a generalization of the classical statistics based on . Importantly, — and consequently the proposed statistics — is rotation-invariant

, which means if a linear transformation is applied to the observations with an arbitrary orthogonal matrix, the statistic remains unchanged. It is a desirable property when not much additional knowledge about

and is available. It should be noted that the two-sample mean tests by Bai and Saranadasa (1996) and Li et al. (2016), together with their generalization to MANOVA, are special cases of the proposed family with and , , respectively.

The present work builds on the work by Li et al. (2016). The theoretical analysis also involves an extension of the analytical framework adopted by Pan and Zhou (2011) in their study of the asymptotic behavior of Hotelling’s

statistic for non-Gaussian observations. However, the current work goes well beyond the existing literature in several aspects. We highlight these as the key contributions of this manuscript: (a) We propose new families of rotation-invariant tests for general linear hypotheses for multivariate regression problems involving high-dimensional response and fixed-dimensional predictor variables that incorporate a flexible regularization scheme to account for the dimensionality of the observations growing proportional to the sample size. (b) Unlike

Li et al. (2016), who assumed sub-Gaussianity, here only the existence of finite fourth moments of the observations is required. (c) Unlike Pan and Zhou (2011), who assumed , is allowed to be fairly arbitrary and subjected only to some standard conditions on the limiting behavior of its spectrum. (d) We carry out a detailed analysis of the power characteristics of the proposed tests. The proposal of a class of local alternatives enables a clear interpretation of the contributions of different parameters in the performance of the test. (e) We develop a data-driven test procedure based on the principle of maximizing asymptotic power under appropriate local alternatives. This principle leads to the definition of a composite test that combines the optimal tests associated with a set of different kinds of local alternatives. The latter formulation is an extension of the data-adaptive test procedure designed by Li et al. (2016) for the two-sample testing problem.

The rest of the paper is organized as follows. Section 2 introduces the asymptotics of the proposed test family both under the null hypothesis and under a class of local alternatives. Using these local alternatives, in Section 3 a data-driven shrinkage selection methodology based on maximizing asymptotic power is developed. In Section 4, an application of the asymptotic theory and the shrinkage selection method is given for the ridge-regularization family. An extension of ridge-regularization to higher orders is also discussed. The results of a simulation study are reported in Section 5. In the Appendix, a proof outline of the main theorem is presented, while technical details and proofs of other theorems are collected in the Supplementary Material, which is available at anson.ucdavis.edu/%7Elihaoran/.

## 2 Asymptotic theory

After giving necessary preliminaries on Random Matrix Theory (RMT), the asymptotic theory of the proposed tests under the null hypothesis and under various local alternative models is presented in this section. For any symmetric matrix , define the Empirical Spectral Distribution (ESD) of by

 FA(τ)=p−1∑pi=11{λi(A)≤τ}.

In the following, stands for the maximum absolute value of the entries of a matrix. The following assumptions are employed.

• (Moment conditions) has i.i.d. entries such that , , ;

• (High-dimensional setting) and are fixed, while such that and ;

• (Boundedness of spectral norm) is non-negative definite; ;

• (Asymptotic stability of ESD) There exists a distribution with compact support in , non-degenerate at zero, such that , as , where denotes the Wasserstein distance between distributions and , defined as

 DW(F1,F2)=supf{∣∣∫fdF1−∫fdF2∣∣:f is 1-Lipschitz}.
• (Asymptotically full rank) is of full rank and converges to a positive definite matrix. Moreover, ;

• (Asymptotically estimable) .

### 2.1 Preliminaries on random matrix theory

Recall that the Stieltjes transform of any function of bounded variation on is defined by

 mG(z)=∫∞−∞dG(x)x−z,z∈C+\coloneqq{u+iv:v>0}.

Minor modifications of a standard RMT result imply that, under Conditions C1C6, the ESD converges almost surely to a nonrandom distribution at all points of continuity of . This limit is determined in such a way that for any , the Stieltjes transform of is the unique solution in of the equation

 m(z)=∫dLΣ(τ)τ(1−γ−γzm(z))−z. (2.1)

Equation (2.1) is often referred to as the Marčenko–Pastur equation. Moreover, pointwise almost surely for , converges to . The convergence holds even when (negative reals) with a smooth extension of to . Readers may refer to Bai and Silverstein (2004) and Paul and Aue (2014) for more details. From now on, for notational simplicity, we shall write as and write as . Note that

 mn,p(z)=p−1tr(ˆΣp−zIp)−1

and define

 Θ(z,γ)={1−γ−γzm(z)}−1. (2.2)

It is known that , for any fixed , has a deterministic equivalent (Bai and Silverstein (2004); Liu et al. (2015); Li et al. (2016)), given by

 {Θ−1(z,γ)Σp−zI}−1,

in the sense that for symmetric matrices bounded in operator norm, as ,

 p−1tr[(ˆΣp−zIp)−1A]−p−1tr[{Θ−1(z,γ)Σp−zI}−1A]→0,~{}~{}~{}with probability 1.

Resolvent and deterministic equivalent will be used frequently in this paper. They will appear for example as Cauchy kernels in contour integrals in various places.

### 2.2 Asymptotics under the null hypothesis

To begin with, for , denote by the Gaussian Orthogonal Ensemble (GOE) defined by (1)  ; (2) , , ; (3) ’s are jointly independent for . Throughout this paper, is assumed to be analytic in an open interval containing

 X\coloneqq[0, limsupp→∞λmax(Σp)(1+√γ)2].

Let to be a closed contour enclosing such that has a complex extension to the interior of . Further use to denote .

###### Theorem 2.1

Suppose C1C6 hold. Under the null hypothesis ,

 √n{M(f)−Ω(f,γ)Iq}⟹Δ1/2(f,γ)W,

where denotes weak convergence and and are as follows.

 Ω(f,γ)=−12πi∮Cf(z)(Θ(z,γ)−1)dz.

See (2.2) for the definition of . For any two analytic functions and ,

 Δ(f1,f2,γ)=2(2πi)2\oiintC2f1(z1)f2(z2)δ(z1,z2,γ)dz1dz2,

and is written as for simplicity. The kernel is such that

 δ(z1,z2,γ) =Θ(z1,γ)Θ(z2,γ)[z1Θ(z1,γ)−z2Θ(z2,γ)z1−z2−1], δ(z,z,γ) =γ{1+zm(z)}Θ3(z,γ)+γz{m(z)+zm′(z)}Θ4(z,γ).

The contour integral is taken counter-clockwise.

Using knowledge of the eigenvalues of the GOE leads to the following statement.

###### Corollary 2.1

Under the conditions of Theorem 2.1, assume further that . Let

 ~λi=√nΔ1/2(f,γ){λi(M(f))−Ω(f,γ)},i=1,…,q.

Then, the limiting joint density function of at is given by

 (2q/2q∏i=1Γ(i/2))−1∏i

Although without closed forms, and do not depend on the choice of used to compute the contour integral. With the resolvent as kernel can be expressed as the integral of on any contour , up to a scaling factor. The quadratic form is then shown to concentrate around , which consequently serves as the integral kernel in . The kernel of is the limit of .

###### Remark 2.1

Two sufficient conditions for are

• ;

• , with for some , and

It would be convenient if and had closed forms in order to avoid computational inefficiencies. Closed forms are available for special cases as shown in the following lemma.

###### Lemma 2.1

When with , the contour integrals in Theorem 2.1 have closed forms, namely, for , , ,

 1(2πi)2\oiintC2∂j1f(z1,ℓ1)∂ℓj11∂j2f(z2,ℓ2)∂ℓj22δ(z1,z2,γ)dz1dz2=∂j1+j2δ(ℓ1,ℓ2,γ)∂ℓj11∂ℓj22.

The results continue to hold when .

Lemma 2.1 indicates that it is possible to have convenient and accurate estimators of the asymptotic mean and variance of under ridge-regularization. The result easily generalizes to the setting when is a linear combination of functions of the form , for any finite collection of ’s. We elaborate on this in Section 4.

To conduct the tests, consistent estimators of and are needed.

###### Lemma 2.2

Let and be the plug-in estimators of and , with estimated by . For general , , , we can estimate and by replacing and with and . Denote the resulting estimators by and . Then,

 √n|ˆΩ(f,γn)−Ω(f,γ)|P⟶0, √n|ˆΔ(f1,f2,γn)−Δ(f1,f2,γ)|P⟶0,

where indicates convergence in probability. Again, we write as .

For the special case of , and , using Lemma 2.1, natural estimators in closed forms are

 ˆΩ(f(j)(x,ℓ),γn)=∂j(ˆΘ(ℓ,γn)−1)∂ℓj, ˆΔ(f(j1)(x,ℓ1),f(j2)(x,ℓ2),γn)=∂j1+j22ˆδ(ℓ1,ℓ2,γn)∂ℓj11∂ℓj22.

In particular, for ,

 ˆΩ(f(x,ℓ),γn)=ˆΘ(ℓ,γn)−1, ˆΔ(f(x,ℓ1),f(x,ℓ2),γn)=2ˆδ(ℓ1,ℓ2,γn).

The estimators are consistent, for any fixed and . Given the eigenvalues of , the computational complexity of calculating the above estimators is .

Recall the definitions of and from Section 1.

###### Theorem 2.2

Suppose C1C6 hold and . Under ,

 ˆTLR(f) \coloneqq√n{1+ˆΩ(f,γn)}q1/2ˆΔ1/2(f,γn)[TLR(f)−qlog{1+ˆΩ(f,γn)}]⟹N(0,1), ˆTLH(f) ˆTBNP(f) \coloneqq√n{1+ˆΩ(f,γn)}2q1/2ˆΔ1/2(f,γn){TBNP(f)−qˆΩ(f,γn)1+ˆΩ(f,γn)}⟹N(0,1).

For any of the three tests, the null hypothesis is rejected at asymptotic level , if , where is the quantile of the standard normal distribution.

### 2.3 Asymptotic power under local alternatives

This subsection deals with the behavior of the proposed family of tests under a host of local alternatives. We start with deterministic alternatives, a framework commonly used in the literature to study the asymptotic power of inferential procedures. Next, we consider a Bayesian framework, using a class of priors that characterize the structure of the alternatives. Because the results to follow simultaneously hold for , and , the unifying notation will be used to refer to each of the test statistics.

#### 2.3.1 Deterministic local alternatives

Consider a sequence of such that, on an open subset of containing ,

 √nCTBT{Θ−1(z,γ)Σp−zI}−1BC⟶D(z,γ) pointwise, as n,p→∞. (2.3)

Observe that and define

 (2.4)
 T=limn→∞CT(n−1XXT)−1C. (2.5)

Note that exists and is non-singular under C5 and C6. If further for any , is non-negative definite.

###### Theorem 2.3

Suppose C1C6 and (2.3) hold, and . Then, as ,

 √nΔ1/2(f,γ){M(f)−Ω(f,γ)Iq}⟹W+H(D,f)Δ1/2(f,γ).

Denote the power functions of at asymptotic level , conditional on , by

 Υ(BC,f)=P(ˆT(f)>ξα∣BC).

The asymptotic behavior of the power functions is described in the following corollary.

###### Corollary 2.2

Under the assumptions of Theorem 2.3, as ,

 Υ(BC,f)⟶Φ(−ξα+tr(H(D,f))q1/2Δ1/2(f,γ)),

where is the standard normal CDF.

###### Remark 2.2

Corollary 2.2 indicates the three proposed statistics have identical asymptotic powers under the assumed local alternatives. This is because the first-order Taylor expansions of , and coincide at . However, the respective empirical powers may differ considerably for moderate sample sizes.

The following remark provides a sufficient condition under which (2.3) is satisfied. Denoting the columns of by , it follows that

 √nCTBT{Θ−1(z,γ)Σp−zI}−1BC=√n[μTi{Θ−1(z,γ)Σp−zIp}−1μj]qi,j=1.
###### Remark 2.3

(a) Let denote the eigen-projection associated with . Suppose that there exists a sequence (in ) of mappings from to , satisfying , , and a mapping continuous on such that, as and for ,

 ∫|Bij;p(x)−Bij;∞(x)|dFΣp(x)→0.

Then, under C4, it follows that (2.3) holds with and

 dij(z,γ)=∫Bij;∞(x)dLΣ(x)xΘ−1(z,γ)−z=∫Bij;∞(x)dLΣ(x)x{1−γ−γzm(z)}−z.

(b) If , then (2.3) is satisfied if , for some constants , . In this case, .

#### 2.3.2 Probabilistic local alternatives

While deterministic local alternatives provide useful information, they are somewhat restrictive for the purpose of a systematic investigation of the power characteristics. Therefore, probabilistic alternatives are considered in the form of a sequence of prior distributions for . This has the added advantage of providing flexibility for incorporating structural information about the regression parameters and the constraints matrices. The proposed formulation of probabilistic alternatives can be seen as an extension of the proposal adopted by Li et al. (2016) in the context of two-sample tests for equality of means. One challenge associated with formulating meaningful alternatives to the hypothesis (1.2), when compared to the two-sample testing problem, is that there are many more plausible ways in which the null hypothesis can be violated. Considering this, we propose a class of alternatives, that on one hand can incorporate a multitude of structures of the parameter , while on the other hand retains analytical tractability in terms of providing interpretable expressions for the local asymptotic power.

Assume the following prior model of with separable covariance

 BC=n−1/4p−1/2RVST, (2.6)

where is a stochastic matrix ( fixed) with independent elements such that , and for some ; is a deterministic matrix and is a fixed matrix. Moreover, let and suppose there is a nonrandom function such that, as , on an open subset of containing ,

 p−1tr{(Θ−1(z,γ)Σp−zI)−1RRT}→h(z,γ) pointwise. (2.7)

Recalling that is the deterministic equivalent of the resolvent , existence of the limit (2.7) also implies that converges pointwise in probability to . Notice also that is the Stieltjes transform of a measure supported on the eigenvalues of .

Model (2.6) leads to a fairly broad covariance design for multi-dimensional random elements, encompassing structures commonly encountered in many application domains, especially in spatio-temporal statistics. We give some representative examples by considering various functional forms of the matrix . Denote by the columns of and by the columns of .

###### Example 2.1

In all that follows takes values in .

• Independent: ;

• Longitudinal: ;

• Moving average: for constants .

Taking the MANOVA problem to illustrate, suppose that the columns of

represent group mean vectors, and suppose

is the matrix that determines successive contrasts among them. Then, is the difference between the means of group and group . Parts (a)–(c) of Example 2.1 correspond then to respectively following an independent, a longitudinal and a moving average process. The row-wise covariance structure is assumed to be such that each has a covariance matrix proportional to . The factor provides the scaling for the tests to have non-trivial local power.

A sufficient condition that leads to (2.7), similar to Remark 2.3, is to postulate the existence of functions satisfying , , and

 ∫|~Bp(x)−~B∞(x)|dFΣp(x)→0

for some function continuous on , where is the th eigenvalue of and is the eigen-projection associated with . Then

 h(z,γ)=∫~B∞(x)dLΣ(x)x{1−γ−γzm(z)}−z. (2.8)

Equations (2.7) and (2.8) indicate that effectively captures the distribution of the total spectral mass of across the spectral coordinates of , also taking into account the dimensionality effect through the aspect ratio . Later, we shall discuss specific classes of the matrices that lead to analytically tractable expressions for , with the structure of linking the parameter under the alternative through (2.6) to the structure of .

Another important feature of the probabilistic model is that it incorporates both dense and sparse alternatives through different specifications of the innovation variables . We consider two special cases.

1. Dense alternative: ;

2. Sparse alternative: , for some , where

is the discrete probability distribution assigning mass

to and mass to the points .

Note that the usual notion of sparsity corresponds to the setting where in addition, . More generally, the second specification above formulates a prior model for that is sparse in the coordinate system determined by . In particular, if is a polynomial in (see Section 3.2 for a discussion), can be seen as sparse in the spectral coordinates of .

###### Theorem 2.4

Suppose that C1C6 hold and . Also suppose that, under , has a prior distribution given by (2.6). Then, the power function of each of the three test statistics satisfies

 Υ(BC,f)L1⟶Φ(−ξα+tr(SSTT−1)q1/2Δ1/2(f,γ)∮C−12πif(z)h(z,γ)dz), (2.9)

as , where is as in (2.5) and indicates -convergence (with respect to the prior measure of ).

###### Remark 2.4

Even if the quantity