Statistical significance in high-dimensional linear models

Introduction

Many data problems nowadays carry the structure that the number $p$ of covariables may greatly exceed sample size $n$, i.e., $p\gg n$. In such a setting, a huge amount of work has been pursued addressing prediction of a new response variable, estimation of an underlying parameter vector and variable selection, see for example the books by Hastie, Tibshirani and Friedman 2009, Bühlmann and van de Geer 2011 or the more specific review article by Fan and Lv 2010. With a few exceptions, see Section 1.3.1, the proposed methods and presented mathematical theory do not address the problem of assigning uncertainties, statistical significance or confidence: thus, the area of statistical hypothesis testing and construction of confidence intervals is largely unexplored and underdeveloped. Yet, such significance or confidence measures are crucial in applications where interpretation of parameters and variables is very important. The focus of this paper is the construction of $p$-values and corresponding multiple testing adjustment for a high-dimensional linear model which is often very useful in $p\gg n$ settings:

We are interested in testing one or many null-hypotheses of the form:

where $G\subseteq\{1,\ldots,p\}$ is a subset of all the indices of the covariables. Of substantial interest is the case where $G=\{j\}$ corresponding to a hypothesis for the individual $j$th regression parameter ($j=1,\ldots,p$). At the other end of the spectrum is the global null-hypothesis where $G=\{1,\ldots,p\}$, and we allow for any $G$ between an individual and the global hypothesis.

We review in this section an important stream of research for high-dimensional linear models. The more familiar reader may skip Section 1.1.

has become tremendously popular for estimation in high-dimensional linear models. The three main themes which have been considered in the past are prediction of the regression surface (and for a new response variable) with corresponding measure of accuracy

estimation of the parameter vector whose quality is assessed by

and variable selection or estimating the support of $\beta^{0}$, denoted by the active set $S_{0}=\{j;\ \beta^{0}_{j}\neq 0,\ j=1,\ldots,p\}$ such that

is large for a selection (estimation) procedure $\hat{S}$.

Greenshtein and Ritov 2004 proved the first result closely related to prediction as measured in (3). Without any conditions on the deterministic design matrix $\mathbf{X}$, except that the columns are normalized such that $(n^{-1}\mathbf{X}^{T}\mathbf{X})_{jj}\equiv 1$, one has with high probability at least $1-2\exp(-t^{2}/2)$:

Such a slow rate of convergence can be improved under additional assumptions on the design matrix $\mathbf{X}$. The ill-posedness of the design matrix can be quantified using the

concept of “modified” eigenvalues. Consider the matrix $\hat{\Sigma}=n^{-1}\mathbf{X}^{T}\mathbf{X}$. The smallest eigenvalue of $\hat{\Sigma}$ is

where $\phi_{0}$ is the compatibility constant (smallest “modified” eigenvalue) of the fixed design matrix $\mathbf{X}$ (Bühlmann and van de Geer 2011, Bühlmann and van de Geer 2011, Cor. 6.2). Again, this holds by assuming Gaussian errors but the result can be extended to non-Gaussian distributions. From (7), we have two immediate implications: from an asymptotic point of view, using $\lambda\asymp\sqrt{\log(p)/n}$ and assuming that $\phi_{0}$ is bounded away from 0,

Furthermore, when making a restrictive assumption for the design, called neighborhood stability, or assuming the equivalent irrepresentable condition, and choosing a suitable $\lambda\gg\sqrt{\log(p)/n}$:

see Meinshausen and Bühlmann 2006, Zhao and Yu 2006, and Wainwright 2009 establishes exact scaling results. The “beta-min” assumption in (10) as well as the irrepresentable condition on the design are restrictive and non-checkable. Furthermore, these conditions are essentially necessary (Meinshausen and Bühlmann 2006 Meinshausen and Bühlmann 2006; Zhao and Yu 2006 Zhao and Yu 2006). Thus, under weaker assumptions, we can only derive a weaker yet useful result about variable screening. Assuming a restricted eigenvalue condition on the fixed design $\mathbf{X}$ and the “beta-min” condition in (10) we still have asymptotically that for $\lambda\asymp\sqrt{\log(p)/n}$:

The cardinality of the estimated active set (typically) satisfies $|\hat{S}(\lambda)|\leq\min(n,p)$: thus if $p\gg n$, we achieve a massive and often useful dimensionality reduction in the original covariates.

We summarize that a slow convergence rate for prediction “always” holds. Assuming some “constrained minimal eigenvalue” condition on the fixed design $\mathbf{X}$, we obtain the fast convergence rate in (8), and an estimation error bound as in (9); with the additional “beta-min” assumption, we obtain the practically useful variable screening property in (11). For consistent variable selection, we necessarily need a (much) stronger condition on the fixed design, and such a strong condition is questionable to be true in a practical problem. Hence variable selection might be a too ambitious goal with the Lasso. That is why the original translation of Lasso (Least Absolute Shrinkage and Selection Operator) may be better re-translated as Least Absolute Shrinkage and Screening Operator. We refer to Bühlmann and van de Geer 2011 for an extensive treatment of the properties of the Lasso.

1.2 Other methods

Of course, the three main inference tasks in a high-dimensional linear model, as described by (3), (4) and (5), can be pursued with other methods than the Lasso.

Quite different from estimation of the high-dimensional parameter vector are variable screening procedures which aim for an analogous property as in (11). Prominent examples include the “Sure Independence Screening” (SIS) method (Fan and Lv 2008), and high-dimensional variable screening or selection properties have been established for forward variable selection (Wang 2009) and for the PC-algorithm (Bühlmann, Kalisch and Maathuis 2010) (“PC” stands for the first names of its inventors, Peter Spirtes and Clark Glymour).

2 Assigning uncertainties and pp-values for high-dimensional regression

At the core of statistical inference is the specification of statistical uncertainties, significance and confidence. For example, instead of having a variable selection result where the probability in (5) is large, we would like to have measures controlling a type I error (false positive selections), including $p$-values which are adjusted for large-scale multiple testing, or construction of confidence intervals or regions. In the high-dimensional setting, answers to these core goals are challenging.

Wasserman and Roeder 2009 propose a procedure for variable selection based on sample splitting. Using their idea and extending it to multiple sample splitting, Meinshausen, Meier and Bühlmann 2009 develop a much more stable method for construction of $p$-values for hypotheses $H_{0,j}:\ \beta^{0}_{j}=0\ (j=1,\ldots,p)$ and for adjusting them in a non-naive way for multiple testing over $p$ (dependent) tests. The main drawback of this procedure is its required “beta-min” assumption in (10). And this is very undesirable since for statistical hypothesis testing, the test should control type I error regardless of the size of the coefficients, while the power of the test should be large if the absolute value of the coefficient would be large: thus, we should avoid assuming (10).

Up to now, for the high-dimensional linear model case with $p\gg n$, it seems that only Zhang and Zhang 2011 managed to construct a procedure which leads to statistical tests for $H_{0,j}$ without assuming a “beta-min” condition.

3 A loose description of our new results

for some known positive definite matrix $\Omega$ and some known constants $\Delta_{j}$. This is the key to derive $p$-values based on this stochastic upper bound. It can be used for construction of $p$-values for individual hypotheses $H_{0,j}$ as well as for more global hypotheses $H_{0,G}$ for any subset $G\subseteq\{1,\ldots,p\}$, including cases where $G$ is (very) large. Furthermore, Theorem 2 justifies a simple approach for controlling the familywise error rate when considering multiple testing of regression hypotheses. Our multiple testing adjustment method itself is closely related to the Westfall–Young permutation procedure (Westfall and Young 1993) and hence, it offers high power, especially in presence of dependence among the many test-statistics (Meinshausen, Maathuis, and Bühlmann 2011).

Our new method as well as the approach in Zhang and Zhang 2011 provide $p$-values (and the latter also confidence intervals) without assuming a “beta-min” condition. Both of them build on using linear estimators and a correction using a non-linear initial estimator such as the Lasso. Using e.g., the Lasso directly leads to the problem of characterizing the distribution of the estimator (in a tractable form): this seems very difficult in high-dimensional settings while it has been worked out for low-dimensional problems (Knight and Fu 2000). The work by Zhang and Zhang 2011 is the only one which studies (sufficiently closely) related questions and goals as in this paper.

The approach by Zhang and Zhang 2011 is based on the idea of projecting the high-dimensional parameter vector to low-dimensional components, as occurring naturally in the hypotheses $H_{0,j}$ about single components, and then proceeding with a linear estimator. This idea is pursued with the “efficient score function” approach from semiparametric statistics (Bickel et al. 1998). The difficulty in the high-dimensional setting is the construction of the score vector $z_{j}$ from which one can derive a confidence interval for $\beta^{0}_{j}$: Zhang and Zhang 2011 propose it as the residual vector from the Lasso when regressing $\mathbf{X}^{(j)}$ against all other variables $\mathbf{X}^{(\setminus j)}$ (where $\mathbf{X}^{(J)}$ denotes the design sub-matrix whose columns correspond to the index set $J\subseteq\{1,\ldots,p\}$). They then prove the asymptotic validity of confidence intervals for finite, sparse linear combinations of $\beta^{0}$. The difference to our work is primarily a rather different construction of the projection where we make use of Ridge estimation with a very simple choice of regularization. A drawback of our method is that, typically, it is not theoretically rate-optimal in terms of power.

Model, estimation and pp-values

Consider one or many null-hypotheses as in (2). We are interested in constructing $p$-values for hypotheses $H_{0,G}$ without imposing a “beta-min” condition as in (10): the statistical test itself will distinguish whether a regression coefficient is small or not.

We consider model (1) with fixed design. Without making additional assumptions on the design matrix $\mathbf{X}$, there is a problem of identifiability. Clearly, if $p>n$ and hence $\operatorname{rank}(\mathbf{X})\leq n<p$, there are different parameter vectors $\theta$ such that $\mathbf{X}\beta^{0}=\mathbf{X}\theta$. Thus, we cannot identify $\beta^{0}$ from the distribution of $Y_{1},\ldots,Y_{n}$ (and fixed design $\mathbf{X}$).

Shao and Deng 2012 give a characterization of identifiability in a high-dimensional linear model (1) with fixed design. Following their approach, it is useful to consider the singular value decomposition

where $A^{-}$ denotes the pseudo-inverse of a squared matrix $A$.

A natural choice of a parameter $\theta^{0}$ such that $\mathbf{X}\beta^{0}=\mathbf{X}\theta^{0}$ is the projection of $\beta^{0}$ onto $\mathcal{R}(\mathbf{X})$. Thus,

Then, of course, $\beta^{0}\in\mathcal{R}(\mathbf{X})$ if and only if $\beta^{0}=\theta^{0}$.

2 Ridge regression

where $\lambda=\lambda_{n}$ is a regularization parameter. By construction of the estimator, $\hat{\beta}\in\mathcal{R}(\mathbf{X})$; and indeed, as discussed below, $\hat{\beta}$ is a reasonable estimator for $\theta^{0}=P_{\mathbf{X}}\beta^{0}$. We denote by

The covariance matrix of the Ridge estimator, multiplied by $n$, is then

a quantity which will appear at many places again. We assume that

Consider the Ridge regression estimator $\hat{\beta}$ in (14) with regularization parameter $\lambda>0$. Assume condition (16), see also (17). Then,

for any $0<C<\infty$, and where $0<L_{C}<M_{C}<\infty$ are constants which depend on $C$ and on the design matrix $\mathbf{X}$ (and hence on $n$ and $p$).

The proof is straightforward using the expression (2.2). The statement 3. says that for a given data-set, the variances of the $\hat{\beta}_{j}$’s remain in a reasonable range even if we choose $\lambda>0$ arbitrarily small; the statement doesn’t imply anything for the behavior as $n$ and $p$ are getting large (as the data and design matrix change). From Proposition 1, we immediately obtain the following result.

Consider the Ridge regression estimator $\hat{\beta}$ in (14) with regularization parameter $\lambda>0$ satisfying

In addition, assume condition (16), see also (17). Then

3 The projection bias and corrected Ridge regression

As discussed in Section 2.1, Ridge regression is estimating the parameter $\theta^{0}=P_{\mathbf{X}}\beta^{0}$ given in (13). Thus, in general, besides the estimation bias governed by the choice of $\lambda$, there is an additional projection bias $B_{j}=\theta^{0}_{j}-\beta^{0}_{j}\ (j=1,\ldots,p)$. Clearly,

In terms of constructing $p$-values, controlling type I error for testing $H_{0,j}$ or $H_{0,G}$ with $j\in G$, the projection bias has only a disturbing effect if $\beta^{0}_{j}=0$ and $\theta^{0}_{j}\neq 0$, and we only have to consider the bias under the null-hypothesis:

The bias $B_{H_{0};j}$ is also the relevant quantity for the case under the non null-hypothesis, see the brief comment after Proposition 2. We can estimate $B_{H_{0};j}$ by

We then have the following representation.

A proof is given in Section .1. We infer from Proposition 2 a representation which could be used not only for testing but also for constructing confidence intervals:

The normalizing factors for the variables $Z_{j}$ bringing them to the $\mathcal{N}(0,1)$-scale are

4 Stochastic bound for the distribution of the corrected Ridge estimator: Asymptotics

We consider a triangular array of observations from a linear model as in (1):

where all the quantities and also the dimension $p=p_{n}$ are allowed to change with $n$. We make the following assumption.

There are constants $\Delta_{j}=\Delta_{j,n}>0$ such that

We will discuss in Section 2.4.1 constructions for such bounds $\Delta_{j}$ (which are typically not negligible). Our next result is the key to obtain a $p$-value for testing the null-hypothesis $H_{0,j}$ or $H_{0,G}$, saying that asymptotically,

where $Z_{1},\ldots,Z_{P}$ are as in Proposition 2.

A proof is given in Section .1. As written above already, due to the third statement in Lemma 1, the condition for $\lambda_{n}$ is reasonable. We note that the distribution of $\max_{j\in G_{n}}(a_{n,p;j}(\sigma)|Z_{j}|+\Delta_{j})$ does not depend on $\sigma$ and can be easily computed via simulation.

We discuss an approach for constructing the bounds $\Delta_{j}$. As mentioned above, they should not involve any unknown quantities so that we can use them for constructing $p$-values from the distribution of $|W|+\Delta_{j}$ or $\max_{j\in G_{n}}(a_{n,p;j}(\sigma)|Z_{j}|+\Delta_{j})$, respectively.

To proceed further, we consider the Lasso as initial estimator. Due to (7) we obtain

A proof follows from (24). We summarize the results as follows.

Assume the conditions of Theorem 1 without condition (A) and the conditions of Lemma 2. Then, when using the Lasso as initial estimator, the statements in Theorem 1 hold.

The construction of the bound in (25) requires the compatibility condition on the design and an upper bound for the sparsity $s_{0}$. While the former is an identifiability condition, and some form of identifiability assumption is certainly necessary, the latter condition about knowing the magnitude of the sparsity is not very elegant. When assuming bounded sparsity $s_{0,n}\leq M<\infty$ for all $n$, we can choose $\xi=0$ with an additional constant $M$ on the right-hand side of (25). In our practical examples in Section 5, we use $\xi=0.05$.

5 PP-values

Our construction of $p$-values is based on the asymptotic distributions in Theorem 1. For an individual hypothesis $H_{0,j}$, we define the $p$-value for the two-sided alternative as

Of course, we could also consider one-sided alternatives with the obvious modification for $P_{j}$. For a more general hypothesis $H_{0,G}$ with $|G|>1$, we use the maximum as test statistics (but other statistics such as weighted sums could be chosen as well) and denote by

where the latter is independent of $\sigma$ and can be easily computed via simulation ($Z_{1},\ldots,Z_{p}$ are as in Proposition 2). Then, the $p$-value for $H_{0,G}$, against the alternative being the complement $H_{0,G}^{c}$, is defined as

Error control follows immediately by the construction of the $p$-values.

Assume the conditions in Theorem 1. Then, for any $0<\alpha<1$,

Furthermore, for any sequence $\alpha_{n}\to 0\ (n\to\infty)$ which converges sufficiently slowly, the statements also hold when replacing $\alpha$ by $\alpha_{n}$.

A discussion about detection power of the method is given in Section 4. Further remarks about these $p$-values are given in Section .4.

We propose to use the estimator $\hat{\sigma}$ from the Scaled Lasso method (Sun and Zhang 2012). Assuming $s_{0}\log(p)/n=o(1)\ (n\to\infty)$ and the compatibility condition for the design, Sun and Zhang 2012 prove that $|\hat{\sigma}/\sigma-1|=o_{P}(1)\ (n\to\infty)$.

Multiple testing

recall that $S_{0}=\{j;\ \beta^{0}_{j}\neq 0\}$ is the set of true active variables. The number of false positives using the nominal significance level $\alpha$ is the denoted by

Consider the variables $Z_{1},\ldots,Z_{p}\sim\mathcal{N}_{p}(0,\sigma^{2}n^{-1}\Omega)$ appearing in Proposition 2 or Theorem 1. Consider the following distribution function:

We first derive familywise error control in an asymptotic sense. For a finite sample result, see Section 6. We consider the framework as in (22).

Assume the conditions in Theorem 1. For the $p$-value in (26) and using the correction in (28) with $\zeta>0$ we have: for $0<\alpha<1$,

[(Multiple testing correction in (28) with $\bolds{\zeta=0}$)] We could modify the correction in (28) using $\zeta=0$: the statement in Theorem 2 can then be derived when making the additional assumption that

where $F_{n,Z}(\cdot)=F_{Z}(\cdot)$ is the distribution function appearing in (28) which depends in the asymptotic framework on $n$ and (mainly on) $p=p_{n}$. Verifying (29) may not be easy for general matrices $\Omega=\Omega_{n,p_{n}}$. However, for the special case where $Z_{1},\ldots,Z_{p}$ are independent,

which is nicely bounded as a function of $u$, over all values of $p$.

2 Multiple testing of general hypotheses

The methodology for testing many general hypotheses $H_{0,G_{j}}$ with $|G_{j}|\geq 1$, $j=1,\ldots,m$ is the same as before. Denote by $S_{0,G}=\{j;\ H_{0,G_{j}}\ \mbox{does not hold}\}$ and by $S_{0,G}^{c}=\{j;\ H_{0,G_{j}}\ \mbox{holds}\}$; note that these sets are determined by the true parameter vector $\beta^{0}$. Since the $p$-value in (27) is of the form $P_{G_{j}}=1-J_{G_{j}}(\hat{\gamma}_{G_{j}})$, we consider

which can be easily computed via simulation (and it is independent of $\sigma$). We then define the corrected $p$-value as

Sufficient conditions for detection

We consider detection of alternatives $H_{0,j}^{c}$ or $H_{0,G}^{c}$ with $|G|>1$. We use again the notation $S_{0}$ as in Section 3 and denote by $a_{n}\gg b_{n}$ that $a_{n}/b_{n}\to\infty\ (n\to\infty)$.

Consider the setting and assumptions as in Theorem 1.

When considering individual hypotheses $H_{0,j}$: if $j\in S_{0}$ with

there exists an $\alpha_{n}\to 0\ (n\to\infty)$ such that

When considering individual hypotheses $H_{0,G}$ with $G=G_{n}$ and $|G_{n}|>1$: if $H_{0,G}^{c}$ holds, with

there exists an $\alpha_{n}\to 0\ (n\to\infty)$ such that

When considering multiple hypotheses $H_{0,j}$: if for all $j\in S_{0}$,

there exists an $\alpha_{n}\to 0\ (n\to\infty)$ such that

If in addition, $a_{n,p;j}(\sigma)\to\infty$ for all $j$ appearing in the conditions on $\beta_{j}^{0}$, we can replace in all the statements 1–3 the “$\gg$” relation by “$\geq\!C$”, where $0<C<\infty$ is a sufficiently large constant.

A proof is given in Section .1. Under the additional assumption of Lemma 2, where the Lasso is used as initial estimator and using the bounds in (25), we obtain the bound (for statement 1 in Theorem 3):

where $0<\xi<1/2$. This can be sharpened using the oracle bound, assuming known order of sparsity:

for some $D>0$ sufficiently large (for example, assuming $s_{0,n}$ is bounded, and replacing $s_{0,n}$ by $1$ and choosing $D>0$ sufficiently large). It then suffices to require

and analogously for the second statement in Theorem 3.

The order of $a_{n,p;j}(\sigma)$ is typically much larger than $\sqrt{n}$ since in high dimensions, $\Omega_{jj}$ is very small. This means that the Ridge estimator $\hat{\beta}_{j}$ has a much faster convergence rate than $1/\sqrt{n}$ for estimating the projected parameter $\theta^{0}_{j}$. This looks counter-intuitive at first sight: the reason for the phenomenon is that $\|\theta^{0}\|_{2}$ can be much smaller than $\|\beta^{0}\|_{2}$ and hence, Ridge regression (which estimates the parameter $\theta^{0}$) is operating on a much smaller scale. This fact is essentially an implication of the first statement in Lemma 1 (without the “$\min_{j}$” part). We can write

where the columns of $U=[U_{jr}]_{j,r=1,\ldots,p}$ contain the $p$ eigenvectors of $\mathbf{X}^{T}\mathbf{X}$, satisfying $\sum_{j=1}^{p}U_{jr}^{2}=1$. For $n\ll p$, only very few, namely $n$ terms, are left in the summation while the normalization for $U_{jr}^{2}$ is over all $p$ terms. For further discussion about the fast convergence rate $a_{n,p;j}(\sigma)^{-1}$, see Section .4.

While $a_{n,p;j}(\sigma)^{-1}$ is usually small, there is compensation with $(P_{\mathbf{X}})_{jj}^{-1}$ which can be rather large. In the detection bound in e.g., the first part of (4), both terms appearing in the maximum are often of the same order of magnitude; see also Figure 3 in Section .4. Assuming such a balance of terms, we obtain in e.g., the first part of (4):

The value of $\kappa_{j}=\max_{k\neq j}|(P_{\mathbf{X}})_{jk}|/|(P_{\mathbf{X}})_{jj}|$ is often a rather small number between 0.05 and 4, see Table 1 in Section 5. For comparison, Zhang and Zhang 2011 establish under some conditions detection for single hypotheses $H_{0,j}$ with $\beta_{j}^{0}$ in the $1/\sqrt{n}$ range. For the extreme case with $G_{n}=\{1,\ldots,p_{n}\}$, we are in the setting of detection of the global hypotheses, see for example Ingster, Tsybakov and Verzelen 2010 for characterizing the detection boundary in case of independent covariables. Here, our analysis of detection is only providing sufficient conditions, for rather general (fixed) design matrices.

Numerical results

For single testing, we construct $p$-values as in (26) or (27) with $\Delta_{j}$ from (25) with $\xi=0.05$. For multiple testing with familywise error control, we consider $p$-values as in (28) with $\zeta=0$ (and $\Delta_{j}$ as above).

We simulate from the linear model as in (1) with $\varepsilon\sim\mathcal{N}_{n}(0,I)$, $n=100$ and the following configurations:

For both $p\in\{500,2500\}$, the fixed design matrix is generated from a realization of $n$ i.i.d. rows from $\mathcal{N}_{p}(0,I)$. Regarding the regression coefficients, we consider active sets $S_{0}=\{1,2,\ldots,s_{0}\}$ with $s_{0}\in\{3,15\}$ and three different strengths of regression coefficients where $\beta^{0}_{j}\equiv b\ (j\in S_{0})$ with $b\in\{0.25,0.5,1\}$.

The same as in (M1) but for both $p\in\{500,2500\}$, the fixed design matrix is generated from a realization of $n$ i.i.d. rows from $\mathcal{N}_{p}(0,\Sigma)$ with $\Sigma_{jk}\equiv 0.8\ (j\neq k)$ and $\Sigma_{jj}=1$.

Here, a pair such as $(3,0.25)$ denotes the values of $s_{0}=3,\ b=0.25$ (where $b$ is the value of the active regression coefficients).

We consider the decision-rule at significance level $\alpha=0.05$

for testing single hypotheses where $P_{j}$ is as in (26) with plugged-in estimate $\hat{\sigma}$. The considered type I error is the average over non-active variables:

2 Values of P𝐗P_{\mathbf{X}}

The detection results in (30) and (4) depend on the ratio $\kappa_{j}=\max_{k\neq j}|(P_{\mathbf{X}})_{jk}|/|(P_{\mathbf{X}})_{jj}|$. We report in Table 1 summary statistics of $\{\kappa_{j}\}_{j}$ for various datasets. We clearly see that the values of $\kappa_{j}$ are typically rather small which implies good detection properties as discussed in Section 4. Furthermore, the values $\max_{k\neq j}|(P_{\mathbf{X}})_{jk}|$ occurring in the construction of $\Delta_{j}$ in Section 2.4.1 are typically very small (not shown here).

3 Real data application

We consider a problem about motif regression for finding the binding sites in DNA sequences of the HIF1$\alpha$ transcription factor. The binding sites are also called motifs, and they are typically 6–15 base pairs (with categorical values $\in\{A,C,G,T\}$) long.

When compared to the Bonferroni–Holm procedure for controlling FWER based on the raw $p$-values as shown in Figure 2(a), we have for the variables with smallest $p$-values:

Thus, for this example, the multiple testing correction as in Section 3 does not provide large improvements in power over the Bonferroni–Holm procedure; but our method is closely related to the Westfall–Young procedure which has been shown to be asymptotically optimal for a broad class of high-dimensional problems (Meinshausen, Maathuis, and Bühlmann, Meinshausen, Maathuis, and Bühlmann 2011).

Finite sample results

We present here finite sample analogues of Theorem 1 and 2. Instead of assumption (A), we assume the following:

There are constants $\Delta_{j}>0$ such that

Similarly, with probability at least $1-\kappa$, for any subset $G\subseteq\{1,\ldots,p\}$ and if $H_{0,G}$ holds:

Theorem 2 is a consequence of the following finite sample result.

A proof is given in Section .1. We immediately get the following bound for $\zeta\geq 0$:

Conclusions

We have proposed a novel construction of $p$-values for individual and more general hypotheses in a high-dimensional linear model with fixed design and Gaussian errors. We have restricted ourselves to max-type statistics for general hypotheses but modifications to e.g., weighted sums are straightforward using the representation in Proposition 2. A key idea is to use a linear, namely the Ridge estimator, combined with a correction for the potentially substantial bias due to the fact that the Ridge estimator is estimating the projected regression parameter vector onto the row-space of the design matrix. The finding that we can “succeed” with a corrected Ridge estimator in a high-dimensional context may come as a surprise, as it is well known that Ridge estimation can be very bad for say prediction. Nevertheless, our bias corrected Ridge procedure might not be optimal in terms of power, as indicated in Section 4.1. The main assumptions we make are the compatibility condition for the design, i.e., an identifiability condition, and knowledge of an upper bound of the sparsity (see Lemma 2). A related idea of using a linear estimator coupled with a bias correction for deriving confidence intervals has been earlier proposed by Zhang and Zhang 2011.

No tuning parameter. Our approach does not require the specification of a tuning parameter, except for the issue that we crudely bound the true sparsity as in (25); we always used $\xi=0.05$, and the Scaled Lasso initial estimator does not require the specification of a regularization parameter. All our numerical examples were run without tuning the method to a specific setting, and error control with our $p$-value approach is often conservative while the power seems reasonable. Furthermore, our method is generic which allows to test for any $H_{0,G}$ regardless whether the size of $G$ is small or large: we present in the Section .2 an additional simulation where $|G|$ is large. For multiple testing correction or for general hypotheses with sets $G$ where $|G|>1$, we rely on the power of simulation since analytical formulae for max-type statistics under dependence seem in-existing: yet, our simulation is extremely simple as we only need to generate dependent multivariate Gaussian random variables.

Small variance of Ridge estimator. As indicated before, it is surprising that corrected Ridge estimation performs rather well for statistical testing. Although the bias due to the projection $P_{\mathbf{X}}$ can be substantial, it is compensated by small variances $\sigma^{2}n^{-1}\Omega_{jj}$ of the Ridge estimator. It is not true that $\Omega_{jj}$’s become large as $p$ increases: that is, the Ridge estimator has small variance for an individual component when $p$ is very large, see Section 4.1. Therefore, the detection power of the method remains reasonably good as discussed in Section 4. Viewed from a different perspective, even though $|(P_{\mathbf{X}})_{jj}\beta^{0}_{j}|$ may be very small, the normalized version $a_{n,p;j}(\sigma)|(P_{\mathbf{X}})_{jj}\beta^{0}_{j}|$ can be sufficiently large for detection since $a_{n,p;j}(\sigma)$ may be very large (as the inverse of the square root of the variance). The values of $P_{\mathbf{X}}$ can be easily computed for a given problem: our analysis about sufficient conditions for detection in Section 4 could be made more complete by invoking random matrix theory for the projection $P_{\mathbf{X}}$ (assuming that $\mathbf{X}$ is a realization of i.i.d. row-vectors whose entries are potentially dependent). However, currently, most of the results on singular values and similar quantities of $\mathbf{X}$ are for the regime $p\leq n$ (Vershynin 2012), which leads in our context to the trivial projection $P_{\mathbf{X}}=I$, or for the regime $p/n\to C$ with $0\leq C<\infty$ (El Karoui 2008).

Extensions. Obvious but partially non-trivial model extensions include random design, non-Gaussian errors or generalized linear models. From a practical point of view, the second and third issue would be most valuable. Relaxing the fixed design assumption makes part of the mathematical arguments more complicated, yet a random design is better posed in terms of identifiability.

Appendix

Proof of Proposition 1 The statement about the bias is given in Shao and Deng 2012 (proof of their Theorem 1). The covariance matrix of $\hat{\beta}$ is

Proof of Proposition 3 (basis for proving Theorem 1) The bound from Proposition 1 for the estimation bias of the Ridge estimator leads to:

By using the representation from Proposition 2, invoking assumption (A′) and assuming that the null-hypothesis $H_{0,j}$ or $H_{0,G}$ holds, respectively, the proof is completed.

Proof of Theorem 1 Due to the choice of $\lambda=\lambda_{n}$ we have that $\|a_{n,p}b(\lambda_{n})\|_{\infty}=o(1)\ (n\to\infty)$. The proof then follows from Proposition 3 and invoking assumption (A) saying that the probabilities for the statements in Proposition 3 converge to 1 as $n\to\infty$.

where in the last inequality we used Proposition 2 and Taylor’s expansion. Thus, on $\mathcal{E}$:

Proof of Theorem 3 Throughout the proof, $\alpha_{n}\to 0$ is converging sufficiently slowly, possibly depending on the context of the different statements we prove. Regarding statement 1: it is sufficient that for $j\in S_{0}$,

From Proposition 2, we see that this can be enforced by requiring

Due to the choice of $\lambda=\lambda_{n}$ (as in Theorem 1), we have $a_{n,p;j}(\sigma)b_{j}(\lambda)\leq\|a_{n,p}(\sigma)b(\lambda)\|_{\infty}=o(1)$. Hence, (35) holds with probability converging to one if

For proving the second statement, we recall that

Using the union bound and the fact that $a_{n,p;j}(\sigma)|Z_{j}|\sim\mathcal{N}(0,1)$ (but dependent over different values of $j$), we have that

The argument is now analogous to the proof of the first statement above, using the representation from Proposition 2.

Regarding the third statement, we invoke the rough bound

with the non-truncated Bonferroni corrected $p$-value at the right-hand side. Hence,

Since this involves a standard Gaussian two-sided tail probability, the inequality can be enforced (for certain slowly converging $\alpha_{n}$) by

The argument is now analogous to the proof of the first statement above, using the representation from Proposition 2.

The fourth statement involves slight obvious modifications of the arguments above.

.2 PP-values for H0,GH_{0,G} with |G||G| large

We report here on a small simulation study for testing $H_{0,G}$ with $G=\{1,2,\ldots,100\}$. We consider model (M2) from Section 5.1 with 4 different configurations and we use the $p$-value from (27) with corresponding decision rule for rejection of $H_{0,G}$ if the $p$-value is smaller or equal to the nominal level 0.05. Table 2 describes the result based on 500 independent simulations (where the fixed design remains the same). The method works well with much better power than multiple testing of individual hypotheses but worse than average power for testing individual hypotheses without multiplicity adjustment (which is not a proper approach). This is largely in agreement with the theoretical results in Theorem 3. Furthermore, the type I error control is good.

.3 Number of false positives in simulated examples

We show in Table 3 the number of false positives $V=V_{0.05}$ in the simulated scenarios where the FWER (among individual hypotheses) was found too large. Although the FWER is larger than 0.05, the number of false positives is relatively small, except for the extreme model (M2), $p=2500$, $s=15$, $b=1$ which has a too large sparsity and a too strong signal strength. For the latter model, we would need to increase $\xi$ in (25) to achieve better error control.

.4 Further discussion about pp-values and bounds Δj\Delta_{j} in assumption (A)

The $p$-values in (26) and (27) are crucially based on the idea of correction with the bounds $\Delta_{j}$ in Section 2.4.1. The essential idea is contained in Proposition 2:

a correction with the bound $\Delta_{j}$ would not be necessary, but of course, it does not hurt in terms of type I error control. If

for some non-degenerate random variable $V$, the correction with the bound $\Delta_{j}$ is necessary and assuming that $\Delta_{j}$ is of the same order of magnitude as $V$, we have a balance between $\Delta_{j}$ and the stochastic term $a_{n,p;j}(\sigma)Z_{j}$. In the last case where

the bound $\Delta_{j}$ would be the dominating element in the $p$-value construction. We show in Figure 3 that there is empirical evidence that (38) applies most often.

Case (39) is comparable to a crude procedure which makes a hard decision about relevance of the underlying coefficients:

and the rejection would be “certain” corresponding to a $p$-value with value equal to $0$; and in case of a “$\leq$” relation, the corresponding $p$-value would be set to one. This is an analogue to the thresholding rule:

Acknowledgements

I would like to thank Cun-Hui Zhang for fruitful discussions and Stephanie Zhang for providing an R-program for the Scaled Lasso.

References