A significance test for the lasso

Introduction.

where $\lambda\geq 0$ is a tuning parameter, controlling the level of sparsity in $\hat{\beta}$. Here, we assume that the columns of $X$ are in general position in order to ensure uniqueness of the lasso solution [this is quite a weak condition, to be discussed again shortly; see also Tibshirani 2013].

There has been a considerable amount of recent work dedicated to the lasso problem, both in terms of computation and theory. A comprehensive summary of the literature in either category would be too long for our purposes here, so we instead give a short summary: for computational work, some relevant contributions are Friedman et al. 2007, Beck and Teboulle 2009, Friedman, Hastie and Tibshirani 2010, Becker, Bobin and Candès 2011, Boyd et al. 2011, Becker, Candès and Grant 2011; and for theoretical work see, for example, Greenshtein and Ritov 2004, Fuchs 2005, Donoho 2006, Candes and Tao 2006, Zhao and Yu 2006, Wainwright 2009, Candès and Plan 2009. Generally speaking, theory for the lasso

is focused on bounding the estimation error $\|X\hat{\beta}-X\beta^{*}\|_{2}^{2}$ or $\|\hat{\beta}-\beta^{*}\|_{2}^{2}$, or ensuring exact recovery of the underlying model, $\operatorname{supp}(\hat{\beta})=\operatorname{supp}(\beta^{*})$ [with $\operatorname{supp}(\cdot)$ denoting the support function]; favorable results in both respects can be shown under the right assumptions on the generative model (1) and the predictor matrix $X$. Strong theoretical backing, as well as fast algorithms, have made the lasso a highly popular tool.

Yet, there are still major gaps in our understanding of the lasso as an estimation procedure. In many real applications of the lasso, a practitioner will undoubtedly seek some sort of inferential guarantees for his or her computed lasso model—but, generically, the usual constructs like $p$-values, confidence intervals, etc., do not exist for lasso estimates. There is a small but growing literature dedicated to inference for the lasso, and important progress has certainly been made, with many methods being based on resampling or data splitting; we review this work in Section 2.5. The current paper focuses on a significance test for lasso models that does not employ resampling or data splitting, but instead uses the full data set as given, and proposes a test statistic that has a simple and exact asymptotic null distribution.

Section 2 defines the problem that we are trying to solve, and gives the details of our proposal—the covariance test statistic. Section 3 considers an orthogonal predictor matrix $X$, in which case the statistic greatly simplifies. Here, we derive its $\operatorname{Exp}(1)$ asymptotic distribution using relatively simple arguments from extreme value theory. Section 4 treats a general (nonorthogonal) $X$, and under some regularity conditions, derives an $\operatorname{Exp}(1)$ limiting distribution for the covariance test statistic, but through a different method of proof that relies on discrete-time Gaussian processes. Section 5 empirically verifies convergence of the null distribution to $\operatorname{Exp}(1)$ over a variety of problem setups. Up until this point, we have assumed that the error variance $\sigma^{2}$ is known; in Section 6, we discuss the case of unknown $\sigma^{2}$. Section 7 gives some real data examples. Section 8 covers extensions to the elastic net, generalized linear models, and the Cox model for survival data. We conclude with a discussion in Section 9.

Significance testing in linear modeling.

Classic theory for significance testing in linear regression operates on two fixed nested models. For example, if $M$ and $M\cup\{j\}$ are fixed subsets of $\{1,\ldots,p\}$, then to test the significance of the $j$th predictor in the model (with variables in) $M\cup\{j\}$, one naturally uses the chi-squared test, which computes the drop in residual sum of squares (RSS) from regression on $M\cup\{j\}$ and $M$,

and compares this to a $\chi_{1}^{2}$ distribution. (Here, $\sigma^{2}$ is assumed to be known; when $\sigma^{2}$ is unknown, we use the sample variance in its place, which results in the $F$-test, equivalent to the $t$-test, for testing the significance of variable $j$.)

Often, however, one would like to run the same test for $M$ and $M\cup\{j\}$ that are not fixed, but the outputs of an adaptive or greedy procedure. Unfortunately, adaptivity invalidates the use of a $\chi_{1}^{2}$ null distribution for the statistic (3). As a simple example, consider forward stepwise regression: starting with an empty model $M=\varnothing$, we enter predictors one at a time, at each step choosing the predictor $j$ that gives the largest drop in residual sum of squares. In other words, forward stepwise regression chooses $j$ at each step in order to maximize $R_{j}$ in (3), over all $j\notin M$. Since $R_{j}$ follows a $\chi_{1}^{2}$ distribution under the null hypothesis for each fixed $j$, the maximum possible $R_{j}$ will clearly be stochastically larger than $\chi_{1}^{2}$ under the null. Therefore, using a chi-squared test to evaluate the significance of a predictor entered by forward stepwise regression would be far too liberal (having type I error much larger than the nominal level). Figure 1(a) demonstrates this point by displaying the quantiles of $R_{1}$ in forward stepwise regression (the chi-squared statistic for the first predictor to enter) versus those of a $\chi_{1}^{2}$ variate, in the fully null case (when $\beta^{*}=0$). A test at the $5\%$ level, for example, using the $\chi^{2}_{1}$ cutoff of $3.84$, would have an actual type I error of about $39\%$.

The failure of standard testing methodology when applied to forward stepwise regression is not an anomaly—in general, there seems to be no direct way

to carry out the significance tests designed for fixed linear models in an adaptive setting. It is important to mention that a simple application of sample splitting can yield proper $p$-values for an adaptive procedure like forward stepwise: for example, run forward stepwise regression on one-half of the observations to construct a sequence of models, and use the other half to evaluate significance via the usual chi-squared test. Some of the related work mentioned in Section 2.5 does essentially this, but with more sophisticated splitting schemes. Our proposal uses the entire data set as given, and we do not consider sample splitting or resampling techniques. Aside from adding a layer of complexity, the use of sample splitting can result in a loss of power in significance testing. Our aim is hence to provide a (new) significance test for the predictor variables chosen adaptively by the lasso, which we describe next.

Before defining our statistic, we briefly review some properties of the lasso path.

The path $\hat{\beta}(\lambda)$ is a continuous and piecewise linear function of $\lambda$, with knots (changes in slope) at values $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{r}\geq 0$ (these knots depend on $y,X$).

At $\lambda=\infty$, the solution $\hat{\beta}(\infty)$ has no active variables (i.e., all variables have zero coefficients); for decreasing $\lambda$, each knot $\lambda_{k}$ marks the entry or removal of some variable from the current active set (i.e., its coefficient becomes nonzero or zero, resp.). Therefore, the active set, and also the signs of active coefficients, remain constant in between knots.

At any point $\lambda$ in the path, the corresponding active set $A=\operatorname{supp}(\hat{\beta}(\lambda))$ of the lasso solution indexes a linearly independent set of predictor variables, that is, $\operatorname{rank}(X_{A})=|A|$, where we use $X_{A}$ to denote the columns of $X$ in $A$.

For a matrix $X$ satisfying the positive cone condition (a restrictive condition that covers, e.g., orthogonal matrices), there are no variables removed from the active set as $\lambda$ decreases and, therefore, the number of knots is $p$.

We can now precisely define the problem that we are trying to solve: at a given step in the lasso path (i.e., at a given knot), we consider testing the significance of the variable that enters the active set. To this end, we propose a test statistic defined at the $k$th step of the path.

We propose the covariance test statistic defined by

Indeed, the natural choice for the tuning parameter in (5) is $\lambda=\lambda_{k+1}$: this allows the $j$th coefficient to have its fullest effect on the fit $X\hat{\beta}$ before the entry of the next variable at $\lambda_{k+1}$ (or possibly, the deletion of a variable from $A$ at $\lambda_{k+1}$).

that is, $T_{k}$ is asymptotically distributed as a standard exponential random variable, given reasonable assumptions on $X$ and the magnitudes of the nonzero true coefficients. [In some cases, e.g., when we have a strict inclusion $A\supsetneq\operatorname{supp}(\beta^{*})$, the use of an $\operatorname{Exp}(1)$ null distribution is actually conservative, because the limiting distribution of $T_{k}$ is stochastically smaller than $\operatorname{Exp}(1)$.] In the above limit, we are considering both $n,p\rightarrow\infty$; in Section 4, we allow for the possibility $p>n$, the high-dimensional case.

See Figure 1(b) for a quantile–quantile plot of $T_{1}$ versus an $\operatorname{Exp}(1)$ variate for the same fully null example ($\beta^{*}=0$) used in Figure 1(a); this shows that the weak convergence to $\operatorname{Exp}(1)$ can be quite fast, as the quantiles are decently matched even for $p=10$. Before proving this limiting distribution in Sections 3 (for an orthogonal $X$) and 4 (for a general $X$), we give an example of its application to real data, and discuss issues related to practical usage. We also derive useful alternative expressions for the statistic, present a connection to degrees of freedom, review related work, and finally, discuss the null hypothesis in more detail.

2 Prostate cancer data example and practical issues.

We consider a training set of 67 observations and 8 predictors, the goal being to predict log of the PSA level of men who had surgery for prostate cancer. For more details, see Hastie, Tibshirani and Friedman 2008 and the references therein. Table 1 shows the results of forward stepwise regression and the lasso. Both methods entered the same predictors in the same order. The forward stepwise $p$-values are smaller than the lasso $p$-values, and would enter four predictors at level $0.05$. The latter would enter only one or maybe two predictors. However, we know that the forward stepwise $p$-values are inaccurate, as they are based on a null distribution that does not account for the adaptive choice of predictors. We now make several remarks.

The above example implicitly assumed that one might stop entering variables into the model when the computed $p$-value rose above some threshold. More generally, our proposed test statistic and associated $p$-values could be used as the basis for multiple testing and false discovery rate control methods for this problem; we leave this to future work.

In the example, the lasso entered a predictor into the active set at each step. For a general $X$, however, a given predictor variable may enter the active set more than once along the lasso path, since it may leave the active set at some point. In this case, we treat each entry as a separate problem. Our test is specific to a step in the path, and not to a predictor variable at large.

One might imagine that such correlation would cause problems for our theory in Sections 3 and 4, which assumes i.i.d. normal errors in the model (1). However, a careful look at the arguments in these sections reveals that the only dependence on $y$ is through $X^{T}y$, the inner products of $y$ with the columns of $X$. Furthermore,

which is the same as it would have been without centering (here $\mathbh{1}\mathbh{1}^{T}$ is the matrix of all $1$s, and we used that the columns of $X$ are centered). Therefore, our arguments in Sections 3 and 4 apply equally well to centered data, and centering has no effect on the asymptotic distribution of $T_{k}$.

By design, the covariance test is applied in a sequential manner, estimating $p$-values for each predictor variable as it enters the model along the lasso path. A more difficult problem is to test the significance of any of the active predictors in a model fit by the lasso, at some arbitrary value of the tuning parameter $\lambda$. We discuss this problem briefly in Section 9.

3 Alternate expressions for the covariance statistic.

Here, we derive two alternate forms for the covariance statistic in (5). The first lends some insight into the role of shrinkage, and the second is helpful for the convergence results that we establish in Sections 3 and 4. We rely on some basic properties of lasso solutions; see, for example, Tibshirani and Taylor 2012, Tibshirani 2013. To remind the reader, we are assuming that $X$ has columns in general position.

For any fixed $\lambda$, if the lasso solution has active set $A=\operatorname{supp}(\hat{\beta}(\lambda))$ and signs $s_{A}=\operatorname{sign}(\hat{\beta}_{A}(\lambda))$, then it can be written explicitly (over active variables) as

In the above expression, the first term $(X_{A}^{T}X_{A})^{-1}X_{A}^{T}y$ simply gives the regression coefficients of $y$ on the active variables $X_{A}$, and the second term $-\lambda(X_{A}^{T}X_{A})^{-1}s_{A}$ can be thought of as a shrinkage term, shrinking the values of these coefficients toward zero. Further, the lasso fitted value at $\lambda$ is

where $P_{A}=X_{A}(X_{A}^{T}X_{A})^{-1}X_{A}^{T}$ denotes the projection onto the column space of $X_{A}$, and $(X_{A}^{T})^{+}=X_{A}(X_{A}^{T}X_{A})^{-1}$ is the (Moore–Penrose) pseudoinverse of $X_{A}^{T}$.

Using the representation (6) for the fitted values, we can derive our first alternate expression for the covariance statistic in (5). If $A$ and $s_{A}$ are the active set and signs just before the knot $\lambda_{k}$, and $j$ is the variable added to the active set at $\lambda_{k}$, with sign $s$ upon entry, then by (6),

and plugging the above two expressions into (5),

Note that the first term above is $y^{T}(P_{A\cup\{j\}}-P_{A})y/\sigma^{2}=(\|y-P_{A}y\|_{2}^{2}-\|y-P_{A\cup\{j\}}y\|_{2}^{2})/\sigma^{2}$, which is exactly the chi-squared statistic for testing the significance of variable $j$, as in (3). Hence, if $A,j$ were fixed, then without the second term, $T_{k}$ would have a $\chi_{1}^{2}$ distribution under the null. But of course $A,j$ are not fixed, and so much like we saw previously with forward stepwise regression, the first term in (2.3) will be generically larger than $\chi_{1}^{2}$, because $j$ is chosen adaptively based on its inner product with the current lasso residual vector. Interestingly, the second term in (2.3) adjusts for this adaptivity: with this term, which is composed of the shrinkage factors in the solutions of the two relevant lasso problems (on $X$ and $X_{A}$), we prove in the coming sections that $T_{k}$ has an asymptotic $\operatorname{Exp}(1)$ null distribution. Therefore, the presence of the second term restores the (asymptotic) mean of $T_{k}$ to $1$, which is what it would have been if $A,j$ were fixed and the second term were missing. In short, adaptivity and shrinkage balance each other out.

This insight aside, the form (2.3) of the covariance statistic leads to a second representation that will be useful for the theoretical work in Sections 3 and 4. We call this the knot form of the covariance statistic, described in the next lemma.

Let $A$ be the active set just before the $k$th step in the lasso path, that is, $A=\operatorname{supp}(\hat{\beta}(\lambda_{k}))$, with $\lambda_{k}$ being the $k$th knot. Also, let $s_{A}$ denote

the signs of the active coefficients, $s_{A}=\operatorname{sign}(\hat{\beta}_{A}(\lambda_{k}))$, $j$ be the predictor that enters the active set at $\lambda_{k}$, and $s$ be its sign upon entry. Then, assuming that

or in other words, all coefficients are active in the reduced lasso problem (4) at $\lambda_{k+1}$ and have signs $s_{A}$, we have

and $s_{A\cup\{j\}}$ is the concatenation of $s_{A}$ and $s$.

The proof starts with expression (2.3), and arrives at (9) through simple algebraic manipulations. We defer it until Appendix .1.

When $X$ satisfies the positive cone condition (which includes $X$ orthogonal), because no variables ever leave the active set in this case. In fact, for $X$ orthogonal, it is straightforward to check that $C(A,s_{A},j,s)=1$, so $T_{k}=\lambda_{k}(\lambda_{k}-\lambda_{k+1})/\sigma^{2}$.

When $k=1$ (we are testing the first variable to enter), as a variable cannot leave the active set right after it has entered. If $k=1$ and $X$ has unit normed columns, $\|X_{i}\|_{2}=1$ for $i=1,\ldots,p$, then we again have $C(A,s_{A},j,s)=1$ (note that $A=\varnothing$), so $T_{1}=\lambda_{1}(\lambda_{1}-\lambda_{2})/\sigma^{2}$.

When $s_{A}=\operatorname{sign}((X_{A})^{+}y)$, that is, $s_{A}$ contains the signs of the least squares coefficients on $X_{A}$, because the same active set and signs cannot appear at two different knots in the lasso path (applied here to the reduced lasso problem on $X_{A}$).

The first and second scenarios are considered in Sections 3 and 4.1, respectively. The third scenario is actually somewhat general and occurs, for example, when $s_{A}=\operatorname{sign}((X_{A})^{+}y)=\operatorname{sign}(\beta_{A}^{*})$; in this case, both the lasso and least squares on $X_{A}$ recover the signs of the true coefficients. Section 4.2 studies the general $X$ and $k\geq 1$ case, wherein this third scenario is important.

4 Connection to degrees of freedom.

There is an interesting connection between the covariance statistic in (5) and the degrees of freedom of a fitting procedure. In the regression setting (1), for an estimate $\hat{y}$ [which we think of as a fitting procedure $\hat{y}=\hat{y}(y)$], its degrees of freedom is typically defined [Efron 1986] as

In words, $\operatorname{df}(\hat{y})$ sums the covariances of each observation $y_{i}$ with its fitted value $\hat{y}_{i}$. Hence, the more adaptive a fitting procedure, the higher this covariance, and the greater its degrees of freedom. The covariance test evaluates the significance of adding the $j$th predictor via something loosely like a sample version of degrees of freedom, across two models: that fit on $A\cup\{j\}$, and that on $A$. This was more or less the inspiration for the current work.

Using the definition (10), one can reason [and confirm by simulation, just as in Figure 1(a)] that with $k$ predictors entered into the model, forward stepwise regression had used substantially more than $k$ degrees of freedom. But something quite remarkable happens when we consider the lasso: for a model containing $k$ nonzero coefficients, the degrees of freedom of the lasso fit is equal to $k$ (either exactly or in expectation, depending on the assumptions) [Efron et al. 2004, Zou, Hastie and Tibshirani 2007, Tibshirani and Taylor 2012]. Why does this happen? Roughly speaking, it is the same adaptivity versus shrinkage phenomenon at play. [Recall our discussion in the last section following the expression (2.3) for the covariance statistic.] The lasso adaptively chooses the active predictors, which costs extra degrees of freedom; but it also shrinks the nonzero coefficients (relative to the usual least squares estimates), which decreases the degrees of freedom just the right amount, so that the total is simply $k$.

5 Related work.

There is quite a lot of recent work related to the proposal of this paper. Wasserman and Roeder 2009 propose a procedure for variable selection and $p$-value estimation in high-dimensional linear models based on sample splitting, and this idea was extended by Meinshausen, Meier and Bühlmann 2009. Meinshausen and Bühlmann 2010 propose a generic method using resampling called “stability selection,” which controls the expected number of false positive variable selections. Minnier, Tian and Cai 2011 use perturbation resampling-based procedures to approximate the distribution of a general class of penalized parameter estimates. One big difference with the work here: we propose a statistic that utilizes the data as given and does not employ any resampling or sample splitting.

Zhang and Zhang 2014 derive confidence intervals for contrasts of high-dimensional regression coefficients, by replacing the usual score vector with the residual from a relaxed projection (i.e., the residual from sparse linear regression). Bühlmann 2013 constructs $p$-values for coefficients in high-dimensional regression models, starting with ridge estimation and then employing a bias correction term that uses the lasso. Even more recently, van de Geer and Bühlmann 2013, Javanmard and Montanari 2013a (Javanmard and Montanari 2013a; Javanmard and Montanari 2013b) all present approaches for debiasing the lasso estimate based on estimates of the inverse covariance matrix of the predictors. (The latter work focuses on the special case of a predictor matrix $X$ with i.i.d. Gaussian rows; the first two consider a general matrix $X$.) These debiased lasso estimates are asymptotically normal, which allows one to compute $p$-values both marginally for an individual coefficient, and simultaneously for a group of coefficients. All of the work mentioned in the present paragraph provides a way to make inferential statements about preconceived predictor variables of interest (or preconceived groups of interest); this is in contrast to our work, which instead deals directly with variables that have been adaptively selected by the lasso procedure. We discuss this next.

6 What precisely is the null hypothesis?

The referees of a preliminary version of this manuscript expressed some confusion with regard to the null distribution considered by the covariance test. Given a fixed number of steps $k\geq 1$ along the lasso path, the covariance test examines the set of variables $A$ selected by the lasso before the $k$th step (i.e., $A$ is the current active set not including the variable to be added at the $k$th step). In particular, the null distribution being tested is

where $\beta^{*}$ is the true underlying coefficient vector in the model (1). For $k=1$, we have $A=\varnothing$ (no variables are selected before the first step), so this reduces to a test of the global null hypothesis: $\beta^{*}=0$. For $k>1$, the set $A$ is random (it depends on $y$), and hence the null hypothesis in (11) is itself a random event. This makes the covariance test a conditional hypothesis test beyond the first step in the path, as the null hypothesis that it considers is indeed a function of the observed data. Statements about its null distribution must therefore be made conditional on the event that $A\supseteq\operatorname{supp}(\beta^{*})$, which is precisely what is done in Sections 3.2 and 4.2.

Compare the null hypothesis in (11) to a null hypothesis of the form

where $S\subseteq\{1,\ldots,p\}$ is a fixed subset. The latter hypothesis, in (12), describes the setup considered by Zhang and Zhang 2014, Bühlmann 2013, van de Geer and Bühlmann 2013, Javanmard and Montanari 2013a (Javanmard and Montanari 2013a; Javanmard and Montanari 2013b). At face value, the hypotheses (11) and (12) may appear similar [the test in (11) looks just like that in (12) with $S=\{1,\ldots,p\}\setminus A$], but they are fundamentally very different. The difference is that the null hypothesis in (11) is random, whereas that in (12) is fixed; this makes the covariance test a conditional hypothesis test, while the tests constructed in all of the aforementioned work are traditional (unconditional) hypothesis tests. It should be made clear that the goal of our work and these works also differ. Our test examines an adaptive subset of variables $A$ deemed interesting by the lasso procedure; for such a goal, it seems necessary to consider a random null hypothesis, as theory designed for tests of fixed hypotheses would not be valid here. In principle, fixed hypothesis tests can be used along with the appropriate correction for multiple comparisons in order to test a random null hypotheses. Aside from being conservative, it is unclear how to efficiently carry out such a procedure when the random null hypothesis consists of a group of coefficients (as opposed to a single one). The main goal of Zhang and Zhang 2014, Bühlmann 2013, van de Geer and Bühlmann 2013, Javanmard and Montanari 2013a (Javanmard and Montanari 2013a; Javanmard and Montanari 2013b), it appears, is to construct a new set of variables, say $\widetilde{A}$, based on testing the hypotheses in (12) with $S=\{j\}$ for $j=1,\ldots,p$. Though the construction of this new set $\widetilde{A}$ may have started from a lasso estimate, it need not be true that $\widetilde{A}$ matches the lasso active set $A$, and ultimately it is this new set $\widetilde{A}$ (and inferential statements concerning $\widetilde{A}$) that these authors consider the point of interest.

An orthogonal predictor matrix XX.

We examine the special case of an orthogonal predictor matrix $X$, that is, one that satisfies $X^{T}X=I$. Even though the results here can be seen as special cases of those for a general $X$ in Section 4, the arguments in the current orthogonal $X$ case rely on relatively straightforward extreme value theory and are hence much simpler than their general $X$ counterparts (which analyze the knots in the lasso path via Gaussian process theory). Furthermore, the $\operatorname{Exp}(1)$ limiting distribution for the covariance statistic translates in the orthogonal case to a few interesting and previously unknown (as far as we can tell) results on the order statistics of independent standard $\chi_{1}$ variates. For these reasons, we discuss the orthogonal $X$ case in detail.

As noted in the discussion following Lemma 1 (see the first point), for an orthogonal $X$, we know that the covariance statistic for testing the entry of the variable at step $k$ in the lasso path is

Again using orthogonality, we rewrite $\|y-X\beta\|_{2}^{2}=\|X^{T}y-\beta\|_{2}^{2}+C$ for a constant $C$ (not depending on $\beta$) in the criterion in (2), and then we can see that the lasso solution at any given value of $\lambda$ has the closed-form:

Letting $U_{j}=X_{j}^{T}y$, $j=1,\ldots,p$, the knots in the lasso path are simply the values of $\lambda$ at which the coefficients become nonzero (i.e., cease to be thresholded),

where $|U_{(1)}|\geq|U_{(2)}|\geq\cdots\geq|U_{(p)}|$ are the order statistics of $|U_{1}|,\ldots,|U_{p}|$ (somewhat of an abuse of notation). Therefore,

Next, we study the special case $k=1$, the test for the first predictor to enter the active set along the lasso path. We then examine the case $k\geq 1$, the test at a general step in the lasso path.

Consider the covariance test statistic for the first predictor to enter the active set, that is, for $k=1$,

We are interested in the distribution of $T_{1}$ under the null hypothesis; since we are testing the first predictor to enter, this is

Under the null, $U_{1},\ldots,U_{p}$ are i.i.d., $U_{j}\sim N(0,\sigma^{2})$, and so $|U_{1}|/\sigma,\ldots,|U_{p}|/\sigma$ follow a $\chi_{1}$ distribution (absolute value of a standard Gaussian). That $T_{1}$ has an asymptotic $\operatorname{Exp}(1)$ null distribution is now given by the next result.

Let $V_{1}\geq V_{2}\geq\cdots\geq V_{p}$ be the order statistics of an independent sample of $\chi_{1}$ variates (i.e., they are the sorted absolute values of an independent sample of standard Gaussian variates). Then

This lemma reveals a remarkably simple limiting distribution for the largest of independent $\chi_{1}$ random variables times the gap between the largest two; we skip its proof, as it is a special case of the following generalization.

If $V_{1}\geq V_{2}\geq\cdots\geq V_{p}$ are the order statistics of an independent sample of $\chi_{1}$ variates, then for any fixed $k\geq 1$,

where $\Phi$ is the standard normal CDF. We first compute

the last equality using Mills’ ratio. Theorem 2.2.1 in de Haan and Ferreira 2006 then implies that, for constants $a_{p}=F^{-1}(1-1/p)$ and $b_{p}=pF^{\prime}(a_{p})$,

where $E_{0}$ is a standard exponential variate, so $-\log{E_{0}}$ has the standard (or type I) extreme value distribution. Hence, according to Theorem 3 in Weissman 1978, for any fixed $k\geq 1$, the random variables $W_{0}=b_{p}(V_{k+1}-a_{p})$ and $W_{i}=b_{p}(V_{i}-V_{i+1})$, $i=1,\ldots,k$, converge jointly:

where $G_{0},E_{1},\ldots,E_{k}$ are independent, $G_{0}$ is Gamma distributed with scale parameter 1 and shape parameter $k$, and $E_{1},\ldots,E_{k}$ are standard exponentials. Now note that

We claim that $a_{p}/b_{p}\rightarrow 1$; this would give the desired result as the second term converges to zero, using $b_{p}\rightarrow\infty$. Writing $a_{p},b_{p}$ more explicitly, we see that $1-1/p=2\Phi(a_{p})-1$, that is, $1-\Phi(a_{p})=1/(2p)$, and $b_{p}=2p\phi(a_{p})$. Using Mills’ inequalities,

Since $a_{p}\rightarrow\infty$, this means that $b_{p}/a_{p}\rightarrow 1$, completing the proof.

Practically, Lemma 3 tells us that under the global null hypothesis $y\sim N(0,\sigma^{2})$, comparing the covariance statistic $T_{k}$ at the $k$th step of the lasso path to an $\operatorname{Exp}(1)$ distribution is increasingly conservative [at the first step, $T_{1}$ is asymptotically $\operatorname{Exp}(1)$, at the second step, $T_{2}$ is asymptotically $\operatorname{Exp}(1/2)$, at the third step, $T_{3}$ is asymptotically $\operatorname{Exp}(1/3)$, and so forth]. This progressive conservatism is favorable, if we place importance on parsimony in the fitted model: we are less and less likely to incur a false rejection of the null hypothesis as the size of the model grows. Moreover, we know that the test statistics $T_{1},T_{2},\ldots$ at successive steps are independent, and hence so are the corresponding $p$-values; from the point of view of multiple testing corrections, this is nearly an ideal scenario.

Of real interest is the distribution of $T_{k}$, $k\geq 1$, not under the global null hypothesis, but rather, under the weaker null hypothesis that all variables excluded from the current lasso model are truly inactive (i.e., they have zero coefficients in the true model). We study this in next section.

2 A general step, k≥1k\geq 1.

We suppose that exactly $k_{0}$ components of the true coefficient vector $\beta^{*}$ are nonzero, and consider testing the entry of the predictor at step $k=k_{0}+1$. Let $A^{*}=\operatorname{supp}(\beta^{*})$ denote the true active set (so $k_{0}=|A^{*}|$), and let $B$ denote the event that all truly active variables are added at steps $1,\ldots,k_{0}$,

We show that under the null hypothesis (i.e., conditional on $B$), the test statistic $T_{k_{0}+1}$ is asymptotically $\operatorname{Exp}(1)$, and further, the test statistic $T_{k_{0}+d}$ at a future step $k=k_{0}+d$ is asymptotically $\operatorname{Exp}(1/d)$.

The basic idea behind our argument is as follows: if we assume that the nonzero components of $\beta^{*}$ are large enough in magnitude, then it is not hard to show (relying on orthogonality, here) that the truly active predictors are added to the model along the first $k_{0}$ steps of the lasso path, with probability tending to one. The test statistic at the $(k_{0}+1)$st step and beyond would therefore depend on the order statistics of $|U_{i}|$ for truly inactive variables $i$, subject to the constraint that the largest of these values is smaller than the smallest $|U_{j}|$ for truly active variables $j$. But with our strong signal assumption, that is, that the nonzero entries of $\beta^{*}$ are large in absolute value, this constraint has essentially no effect, and we are back to studying the order statistics from a $\chi_{1}$ distribution, as in the last section. This is made precise below.

The same convergence in distribution holds conditionally on $B$.

Note that $U_{j}\sim N(\beta^{*}_{j},\sigma^{2})$, independently for $j=1,\ldots,p$. For $j\in A^{*}$,

Hence, we are essentially back in the setting of the last section, and the desired convergence result follows from the same arguments as those for Lemma 3.

A general predictor matrix XX.

In this section, we consider a general predictor matrix $X$, with columns in general position. Recall that our proposed covariance test statistic (5) is closely intertwined with the knots $\lambda_{1}\geq\cdots\geq\lambda_{r}$ in the lasso path, as it was defined in terms of difference between fitted values at successive knots. Moreover, Lemma 1 showed that (provided there are no sign changes in the reduced lasso problem over $[\lambda_{k+1},\lambda_{k}]$) this test statistic can be expressed even more explicitly in terms of the values of these knots. As was the case in the last section, this knot form is quite important for our analysis here. Therefore, it is helpful to recall [Efron et al. 2004; Tibshirani 2013] the precise formulae for the knots in the lasso path. If $A$ denotes the active set and $s_{A}$ denotes the signs of active coefficients at a knot $\lambda_{k}$,

then the next knot $\lambda_{k+1}$ is given by

where recall $P_{A}=X_{A}(X_{A}^{T}X_{A})^{-1}X_{A}^{T}$, and $(X_{A}^{T})^{+}=X_{A}(X_{A}^{T}X_{A})^{-1}$; and

As we did in Section 3 with the orthogonal $X$ case, we begin by studying the asymptotic distribution of the covariance statistic in the special case $k=1$ (i.e., the first model along the path), wherein the expressions for the next knot (14), (15), (16) greatly simplify. Following this, we study the more difficult case $k\geq 1$. For the sake of readability, we defer the proofs and most technical details until the Appendix.

We assume here that $X$ has unit normed columns: $\|X_{i}\|_{2}=1$, for $i=1,\ldots,p$; we do this mostly for simplicity of presentation, and the generalization to a matrix $X$ whose columns are not unit normed is given in the next section (though the exponential limit is now a conservative upper bound). As per our discussion following Lemma 1 (see the second point), we know that the first predictor to enter the active set along the lasso path cannot leave at the next step, so the constant sign condition (8) holds, and by Lemma 1 the covariance statistic for testing the entry of the first variable can be written as

(the leading factor $C$ being equal to one since we assumed that $X$ has unit normed columns). Now let $U_{j}=X_{j}^{T}y$, $j=1,\ldots,p$, and $R=X^{T}X$. With $\lambda_{0}=\infty$, we have $A=\varnothing$, and trivially, no variables can leave the active set. The first knot is hence given by (15), which can be expressed as

Letting $j_{1},s_{1}$ be the first variable to enter and its sign (i.e., they achieve the maximum in the above expression), and recalling that $j_{1}$ cannot leave the active set immediately after it has entered, the second knot is again given by (15), written as

The general position assumption on $X$ implies that $|R_{j,j_{1}}|<1$, and so $1-ss_{1}R_{j,j_{1}}>0$, all $j\neq j_{1}$, $s\in\{-1,1\}$. It is easy to show then that the indicator inside the maximum above can be dropped, and hence

Our goal now is to calculate the asymptotic distribution of $T_{1}=\lambda_{1}(\lambda_{1}-\lambda_{2})/\sigma^{2}$, with $\lambda_{1}$ and $\lambda_{2}$ as above, under the null hypothesis; to be clear, since we are testing the significance of the first variable to enter along the lasso path, the null hypothesis is

The strategy that we use here for the general $X$ case—which differs from our extreme value theory approach for the orthogonal $X$ case—is to treat the quantities inside the maxima in expressions (17), (18) for $\lambda_{1},\lambda_{2}$ as discrete-time Gaussian processes. First, we consider the zero mean Gaussian process

We can easily compute the covariance function of this process:

where the expectation is taken over the null distribution in (19). From (17), we know that the first knot is simply

In addition to (20), we consider the process

An important property: for fixed $j_{1},s_{1}$, the entire process $h^{(j_{1},s_{1})}(j,s)$ is independent of $g(j_{1},s_{1})$. This can be seen by verifying that

and noting that $g(j_{1},s_{1})$ and $h^{(j_{1},s_{1})}(j,s)$, all $j\neq j_{1}$, $s\in\{-1,1\}$, are jointly normal. Now define

and from the above we know that for fixed $j_{1},s_{1}$, $M(j_{1},s_{1})$ is independent of $g(j_{1},s_{1})$. If $j_{1},s_{1}$ are instead treated as random variables that maximize $g(j,s)$ (the argument maximizers being almost surely unique), then from (18) we see that the second knot is $\lambda_{2}=M(j_{1},s_{1})$. Therefore, to study the distribution of $T_{1}=\lambda_{1}(\lambda_{1}-\lambda_{2})/\sigma^{2}$, we are interested in the random variable

It turns out that this event, which concerns the argument maximizers of $g$, can be rewritten as an event concerning only the relative values of $g$ and $M$ [see Taylor, Takemura and Adler 2005 for the analogous result for continuous-time processes].

With $g,M$ as defined in (20), (21), (22), we have

This is an important realization because the dual representation $\{g(j_{1},s_{1})>M(j_{1},s_{1})\}$ is more tractable, once we partition the space over the possible argument minimizers $j_{1},s_{1}$, and use the fact that $M(j_{1},s_{1})$ is independent of $g(j_{1},s_{1})$ for fixed $j_{1},s_{1}$. In this vein, we express the distribution of $T_{1}=\lambda_{1}(\lambda_{1}-\lambda_{2})/\sigma^{2}$ in terms of the sum

The terms in the above sum can be simplified: dropping for notational convenience the dependence on $j_{1},s_{1}$, we have

where $u(a,b)=(b+\sqrt{b^{2}+4a})/2$, which follows by simply solving for $g$ in the quadratic equation $g(g-M)/\sigma^{2}=t$. Therefore,

We now examine the term inside the braces in (4.1), the difference between a ratio of normal survival functions and $e^{-t}$; our next lemma shows that this term vanishes as $m\rightarrow\infty$.

Hence, loosely speaking, if each $M(j_{1},s_{1})\rightarrow\infty$ fast enough as $p\rightarrow\infty$, then the right-hand side in (4.1) converges to zero, and $T_{1}$ converges weakly to $\operatorname{Exp}(1)$. This is made precise below.

Consider $M(j_{1},s_{1})$ defined in (21), (22) over $j_{1}=1,\ldots,p$ and $s_{1}\in\{-1,1\}$. If for any fixed $m_{0}>0$

The assumption in (25) is written in terms of random variables whose distributions are induced by the steps along the lasso path; to make our assumptions more transparent, we show that (25) is implied by a conditional variance bound involving the predictor matrix $X$ alone, and arrive at the main result of this section.

and the size of $S$ growing faster than $\log{p}$,

An example of a matrix $X$ that does not satisfy (26) and (27) is one with fixed rank as $p$ grows. (This, of course, would also not satisfy the general position assumption.) In this case, we would not be able to find a subset of the variables $U_{i}=X_{i}^{T}y$, $i=1,\ldots,p$, that is both linearly independent and has size larger than $r=\operatorname{rank}(X)$, which violates the conditions. We note that in general, since $|S|\leq\operatorname{rank}(X)\leq n$, and $|S|/\log{p}\rightarrow\infty$, conditions (26) and (27) require that $n/\log{p}\rightarrow\infty$.

2 A general step, k≥1k\geq 1.

In this section, we no longer assume that $X$ has unit normed columns (in any case, this provides no simplification in deriving the null distribution of the test statistic at a general step in the lasso path). Our arguments here have more or less the same form as they did in the last section, but overall the calculations are more complicated.

Fix an integer $k_{0}\geq 0$, subset $A_{0}\subseteq\{1,\ldots,p\}$ containing the true active set $A_{0}\supseteq A^{*}=\operatorname{supp}(\beta^{*})$, and sign vector $s_{A_{0}}\in\{-1,1\}^{|A_{0}|}$. Consider the event

First note that on $B$, we have $s_{A}=\operatorname{sign}((X_{A})^{+}y)$, and as discussed in the third point following Lemma 1, this implies that the solution of the reduced problem (4) on $X_{A}$ cannot incur any sign changes over the interval $[\lambda_{k},\lambda_{k+1}]$. Hence, we can apply Lemma 1 to write the covariance statistic on $B$ as

$C(A,s_{A},j_{k},s_{k})=\|(X_{A\cup\{j_{k}\}}^{T})^{+}s_{A\cup\{j_{k}\}}-(X_{A}^{T})^{+}s_{A}\|_{2}^{2}$, $A$ and $s_{A}$ are the active set and signs at step $k-1$, and $j_{k}$ is the variable added to the active set at step $k$, with sign $s_{k}$. Now, analogous to our definition in the last section, we define the discrete-time Gaussian process

For any fixed $A,s_{A}$, the above process has mean zero provided that $A\supseteq A^{*}$. Additionally, for any such fixed $A,s_{A}$, we can compute its covariance function

Note that on the event $B$, the $k$th knot in the lasso path is

For fixed $j_{k},s_{k}$, we also consider the process

(above, $s_{A\cup\{j_{k}\}}$ is the concatenation of $s_{A}$ and $s_{k}$) and its achieved maximum value, subject to being less than the maximum of $g^{(A,s_{A})}$,

If $j_{k},s_{k}$ indeed maximize $g^{(A,s_{A})}$, that is, they correspond to the variable added to the active set at $\lambda_{k}$ and its sign (note that these are almost surely unique), then on $B$, we have $\lambda_{k+1}=M^{(A,s_{A})}(j_{k},s_{k})$. To study the distribution of $T_{k}$ on $B$, we are therefore interested in the random variable

where we have replaced all instances of $A$ and $s_{A}$ on the right-hand side above with the fixed subset $A_{0}$ and sign vector $s_{A_{0}}$. This is a helpful simplification, because in what follows we may now take $A=A_{0}$ and $s_{A}=s_{A_{0}}$ as fixed, and consider the distribution of the random processes $g^{(A_{0},s_{A_{0}})}$ and $M^{(A_{0},s_{A_{0}})}$. With $A=A_{0}$ and $s_{A}=s_{A_{0}}$ fixed, we drop the notational dependence on them and write these processes as $g$ and $M$. We also write the scaling factor $C(A_{0},s_{A_{0}},j_{k},s_{k})$ as $C(j_{k},s_{k})$.

The setup in (4.2) looks very much like the one in the last section [and to draw an even sharper parallel, the scaling factor $C(j_{k},s_{k})$ is actually equal to one over the variance of $g(j_{k},s_{k})$, meaning that $\sqrt{C(j_{k},s_{k})}\cdot g(j_{k},s_{k})$ is standard normal for fixed $j_{k},s_{k}$, a fact that we will use later in the proof of Lemma 8]. However, a major complication is that $g(j_{k},s_{k})$ and $M(j_{k},s_{k})$ are no longer independent for fixed $j_{k},s_{k}$. Next, we derive a dual representation for the event (34) (analogous to Lemma 4 in the last section), introducing a triplet of random variables $M^{+},M^{-},M^{0}$—it turns out that $g$ is independent of this triplet, for fixed $j_{k},s_{k}$.

Let $g$ be as defined in (29) (with $A,s_{A}$ fixed at $A_{0},s_{A_{0}}$). Let $\Sigma_{j,j^{\prime}}$ denote the covariance function of $g$ [short form for the expression in (30)]. To be perfectly clear, here $\Sigma_{j,j^{\prime}}$ actually depends on $s,s^{\prime}$, but our notation suppresses this dependence for brevity. Define

the event $E(j_{k},s_{k})$ in (34), that $j_{k},s_{k}$ maximize $g$, can be written as an intersection of events involving $M^{+},M^{-},M^{0}$:

As a result of Lemma 7, continuing from (4.2), we can decompose the tail probability of $T_{k}$ as

A key point here is that, for fixed $j_{k},s_{k}$, the triplet $M^{+}(j_{k},s_{k})$, $M^{-}(j_{k},s_{k})$, $M^{0}(j_{k},s_{k})$ is independent of $g(j_{k},s_{k})$, which is true because

and $g(j_{k},s_{k})$, along with $g(j,s)-(\Sigma_{j_{k},j}/\Sigma_{j_{k},j_{k}})g(j_{k},s_{k})$, for all $j,s$, form a jointly Gaussian collection of random variables. If we were to now replace $M$ by $M^{+}$ in the first line of (4.2), and define a modified statistic $\widetilde{T}_{k}$ via its tail probability,

Consider $g$ as defined in (29) (with $A,s_{A}$ fixed at $A_{0},s_{A_{0}}$), and $M^{+},M^{-},M^{0}$ as defined in (), (), (). Assume that for any fixed $m_{0}$,

Consider $g,M$ as defined in (29), (4.2), (4.2) (with $A,s_{A}$ fixed at $A_{0},s_{A_{0}}$), and consider $M^{+}$ as defined in (). Then for any fixed $j_{k},s_{k}$, on the event $E(j_{k},s_{k})$ in (34), we have

Though Lemma 9 establishes a (conservative) exponential limit for the covariance statistic $T_{k}$, it does so by enforcing assumption (42), which is phrased in terms of the tail distribution of a random process defined at the $k$th step in the lasso path. We translate this into an explicit condition on the covariance structure in (30), to make the stated assumptions for exponential convergence more concrete.

Assume that the diagonal elements in $R$ are all of the same order, that is, $R_{ii}/R_{jj}\leq C$ for all $i,j$ and some constant $C>0$. Finally assume that, for each fixed $j\notin A_{0}$, there is a set $S\subseteq\{1,\ldots,p\}\setminus(A_{0}\cup\{j\})$ such that for all $i\in S$,

where $\delta>0$ is a constant (not depending on $j$), and the size of $S$ grows faster than $\log{p}$,

Some readers will likely recognize condition (43) as that of mutual incoherence or strong irrepresentability, commonly used in the lasso literature on exact support recovery [see, e.g., Wainwright 2009, Zhao and Yu 2006]. This condition, in addition to a lower bound on the magnitudes of the true coefficients, is sufficient for the lasso solution to recover the true active set $A^{*}$ with probability tending to one, at a carefully chosen value of $\lambda$. It is important to point out that we do not place any requirements on the magnitudes of the true nonzero coefficients; instead, we assume directly that the lasso converges (with probability approaching one) to some fixed model defined by $A_{0},s_{A_{0}}$ at the $(k_{0})$th step in the path. Here, $A_{0}$ is large enough that it contains the true support, $A_{0}\supseteq A^{*}$, and the signs $s_{A_{0}}$ are arbitrary—they may or may not match the signs of the true coefficients over $A_{0}$. In a setting in which the nonzero coefficients in $\beta^{*}$ are well separated from zero, a condition quite similar to the irrepresentable condition can be used to show that the lasso converges to the model with support $A_{0}=A^{*}$ and

signs $s_{A_{0}}=\operatorname{sign}(\beta_{A_{0}}^{*})$, at step $k_{0}=|A_{0}|$ of the path. Our result extends beyond this case, and allows for situations in which the lasso model converges to a possibly larger set of “screened” variables $A_{0}$, and fixed signs $s_{A_{0}}$.

In fact, one can modify the above arguments to account for the case that $A_{0}$ does not contain the entire set $A^{*}$ of truly nonzero coefficients, but rather, only the “strong” coefficients. While “strong” is rather vague, a more precise way of stating this is to assume that $\beta^{*}$ has nonzero coefficients both large and small in magnitude, and with $A_{0}$ corresponding to the set of large coefficients, we assume that the (left-out) small coefficients must be small enough that the mean of the process $g$ in (29) (with $A=A_{0}$ and $s_{A}=s_{A_{0}}$) grows much faster than $M^{+}$. The details, though not the main ideas, of the arguments would change, and the result would still be a conservative exponential limit for the covariance statistic $T_{k}$ at step $k=k_{0}+1$. We may pursue this extension in future work.

Simulation of the null distribution.

We investigate the null distribution of the covariance statistic through simulations, starting with an orthogonal predictor matrix $X$, and then considering more general forms of $X$.

Figure 3 shows the results for testing the 5th, 6th and 7th predictors to enter the lasso model. An $\operatorname{Exp}(1)$-based test will now be conservative: at a nominal $5\%$ level, the actual type I errors are about $1\%$, $0.2\%$ and $0.0\%$, respectively. The solid line has slope 1, and the broken lines have slopes $1/2,1/3,1/4$, as predicted by Theorem 1.

2 General predictor matrix.

In Table 2, we simulated null data (i.e., $\beta^{*}=0$), and examined the distribution of the covariance test statistic $T_{1}$ for the first predictor to enter. We varied the numbers of predictors $p$, correlation parameter $\rho$, and structure of the predictor correlation matrix. In the first two correlation setups, the correlation between each pair of predictors was $\rho$, in the data and population, respectively. In the $\mathit{AR}(1)$ setup, the correlation between predictors $j$ and $j^{\prime}$ is $\rho^{|j-j^{\prime}|}$. Finally, in the block diagonal setup, the correlation matrix has two equal-sized blocks, with population correlation $\rho$ in each block. We computed the mean, variance and tail probability of the covariance statistic $T_{1}$ over $500$ simulated data sets for each setup. We see that the $\operatorname{Exp}(1)$ distribution is a reasonably good approximation throughout.

In Table 3, the setup was the same as in Table 2, except that we set the first $k$ coefficients of the true coefficient vector equal to 4, and the rest zero, for $k=1,2,3$. The dimensions were also fixed at $n=100$ and $p=50$. We computed the mean, variance, and tail probability of the covariance statistic $T_{k+1}$ for entering the next (truly inactive) $(k+1)$st predictor, discarding those simulations in which a truly inactive predictor was selected in the first $k$ steps. (This occurred $1.7\%$, $4.0\%$ and $7.0\%$ of the time, resp.) Again, we see that the $\operatorname{Exp}(1)$ approximation is reasonably accurate throughout.

In Figure 4, we estimate the power curves for significance testing via the drop in RSS test for forward stepwise regression, and the covariance test for the lasso. In the former, we use simulation-derived cutpoints, and in the latter we use the theoretically-based $\operatorname{Exp}(1)$ cutpoints, to control the type I error at the 5% level. We find that the tests have similar power, though the cutpoints for forward stepwise would not be typically available in practice. For more details, see the figure caption.

The case of unknown σ2\sigma^{2}.

Up until now, we have assumed that the error variance $\sigma^{2}$ is known; in practice it will typically be unknown. In this case, provided that $n>p$, we can easily estimate it and proceed by analogy to standard linear model theory. In particular,

This follows because $F_{k}=T_{k}/(\hat{\sigma}^{2}/\sigma^{2})$, the numerator $T_{k}$ being asymptotically $\operatorname{Exp}(1)=\chi_{2}^{2}/2$, the denominator $\hat{\sigma}^{2}/\sigma^{2}$ being asymptotically $\chi_{n-p}^{2}/\penalty(n-p)$, and we claim that the two are independent. Why? Note that the lasso solution path is unchanged if we replace $y$ by $P_{X}y$, so the lasso fitted values in $T_{k}$ are functions of $P_{X}y$; meanwhile, $\hat{\sigma}^{2}$ is a function of $(I-P_{X})y$. The quantities $P_{X}y$ and $(I-P_{X})y$ are uncorrelated, and hence independent (recalling normality of $y$), so $T_{k}$ and $\hat{\sigma}^{2}$ are functions of independent quantities and, therefore, independent.

As an example, consider one of the setups from Table 2, with $n=100$, $p=80$ and predictor correlation of the $\mathit{AR}(1)$ form $\rho^{|j-j^{\prime}|}$. The true model is null, and we test the first predictor to enter along the lasso path. (We choose $n,p$ of roughly equal sizes here to expose the differences between the $\sigma^{2}$ known and unknown cases.) Table 4 shows the results of 1000 simulations from each of the $\rho=0$ and $\rho=0.8$ scenarios. We see that with $\sigma^{2}$ estimated, the $F_{2,n-p}$ distribution provides a more accurate finite-sample approximation than does $\operatorname{Exp}(1)$.

When $p\geq n$, estimation of $\sigma^{2}$ is not nearly as straightforward; one idea is to estimate $\sigma^{2}$ from the least squares fit on the support of the model selected by cross-validation. One would then hope that the resulting statistic, with this plug-in estimate of $\sigma^{2}$, is close in distribution to $F_{2,n-r}$ under the null, where $r$ is the size of the model chosen by cross-validation. This is by analogy to the low-dimensional $n>p$ case in (48), but is not supported by rigorous theory. Simulations (withheld for brevity) show that this approximation is not too far off, but that the variance of the observed statistic is sometimes inflated compared that of an $F_{2,n-r}$ distribution (this unaccounted variability is likely due to the model selection process via cross-validation). Other authors have argued that using cross-validation to estimate $\sigma^{2}$ when $p\gg n$ is not necessarily a good approach, as it can be anti-conservative; see, for example, Fan, Guo and Hao 2012, Sun and Zhang 2012 for alternative techniques. In future work, we will address the important issue of estimating $\sigma^{2}$ in the context of the covariance statistic, when $p\geq n$.

Real data examples.

We demonstrate the use of covariance test with some real data examples. As mentioned previously, in any serious application of significance testing over many variables (many steps of the lasso path), we would need to consider the issue of multiple comparisons, which we do not here. This is a topic for future work.

Table 5 shows the results for the wine quality data taken from the UCI database. There are $p=11$ predictors, and $n=1599$ observations, which we split randomly into approximately equal-sized training and test sets. The outcome is a wine quality rating, on a scale between 0 and 10. The table shows the training set $p$-values from forward stepwise regression (with the chi-squared test) and the lasso (with the covariance test). Forward stepwise enters 6 predictors at the 0.05 level, while the lasso enters only 3.

In the left panel of Figure 5, we repeated this $p$-value computation over 500 random splits into training test sets. The right panel shows the corresponding test set prediction error for the models of each size. The lasso test error decreases sharply once the 3rd predictor is added, but then somewhat flattens out from the 4th predictor onward; this is in general qualitative agreement with the lasso $p$-values in the left panel, the first 3 being very small, and the 4th $p$-value being about 0.2. This also echoes the well-known difference between hypothesis testing and minimizing prediction error. For example, the $C_{p}$ statistic stops entering variables when the $p$-value is larger than about 0.16.

2 HIV data.

Rhee et al. (Rhee et al. 2003) study six nucleotide reverse transcriptase inhibitors (NRTIs) that are used to treat HIV-1. The target of these drugs can become resistant through mutation, and they compare a collection of models for predicting the (log) susceptibility of the drugs, a measure of drug resistance, based on the location of mutations. We focused on the first drug (3TC), for which there are $p=217$ sites and $n=1057$ samples. To examine the behavior of the covariance test in the $p>n$ setting, we divided the data at random into training and test sets of size 150 and 907, respectively, a total of 50 times. Figure 6 shows the results, in the same format as Figure 5. We used the model chosen by cross-validation to estimate $\sigma^{2}$. The covariance test for the lasso suggests that there are only one or two important predictors (in marked contrast to the chi-squared test for forward stepwise), and this is confirmed by the test error plot in the right panel.

Extensions.

We discuss some extensions of the covariance statistic, beyond significance testing for the lasso. The proposals here are supported by simulations [in terms of having an $\operatorname{Exp}(1)$ null distribution], but we do not offer any theory. This may be a direction for future work.

The elastic net estimate [Zou and Hastie 2005] is defined as

with $U_{j}=X_{j}^{T}y$, $j=1,\ldots,p$. This means that for an orthogonal $X$, under the null,

and one is tempted to use this approximation beyond the orthogonal setting as well. In Figure 7, we evaluated the distribution of $(1+\gamma)T_{1}$ (for the first predictor to enter), for orthogonal and correlated scenarios, and for three different values of $\gamma$. Here, $n=100$, $p=10$ and the true model was null. It seems to be reasonably close to $\operatorname{Exp}(1)$ in all cases.

2 Generalized linear models and the Cox model.

This is the implicit concept used by Efron 1986 in his definition of the “optimism” of the training error. The same idea could be used to define degrees of freedom for the penalized estimate in (50) for any $\lambda>0$, and this motivates the definition of the covariance statistic, as follows. If the tuning parameter value $\lambda=\lambda_{k}$ marks the entry of a new predictor into the active set $A$, then we define the covariance statistic

we can analogously define the covariance test statistic at a knot $\lambda_{k}$, marking the entry of a predictor into the active set $A$, as

Discussion.

We proposed a simple covariance statistic for testing the significance of predictor variables as they enter the active set, along the lasso solution path. We showed that the distribution of this statistic is asymptotically $\operatorname{Exp}(1)$, under the null hypothesis that all truly active predictors are contained in the current active set. [See Theorems 1, 2 and 3; the conditions required for this convergence result vary depending on the step $k$ along the path that we are considering, and the covariance structure of the predictor matrix $X$; the $\operatorname{Exp}(1)$ limiting distribution is in some cases a conservative upper bound under the null.] Such a result accounts for the adaptive nature of the lasso procedure, which is not true for the usual chi-squared test (or $F$-test) applied to, for example, forward stepwise regression.

We feel that our work has shed light not only on the lasso path (as given by LARS), but also, at a high level, on forward stepwise regression. Both the lasso and forward stepwise start by entering the predictor variable most correlated with the outcome (thinking of standardized predictors), but the two differ in what they do next. Forward stepwise is greedy, and once it enters this first variable, it proceeds to fit the first coefficient fully, ignoring the effects of other predictors. The lasso, on the other hand, increases (or decreases) the coefficient of the first variable only as long as its correlation with the residual is larger than that of the inactive predictors. Subsequent steps follow similarly. Intuitively, it seems that forward stepwise regression inflates coefficients unfairly, while the lasso takes more appropriately sized steps. This intuition is confirmed in one sense by looking at degrees of freedom (recall Section 2.4). The covariance test and its simple asymptotic null distribution reveal another way in which the step sizes used by the lasso are “just right.”

The problem of assessing significance in an adaptive linear model fit by the lasso is a difficult one, and what we have presented in this paper by no means a complete solution. We describe some current work and ideas for future projects below.

Significance test for generic lasso models. A natural direction to consider is the generic lasso testing problem: given a lasso model computed at some fixed value of $\lambda$, how do we carry out a significance test for each predictor in the active set? Work on this is in progress.

Nonasymptotic null distributions. A geometric characterization of the first knot in the lasso path provides an alternative test for the global null hypothesis, $\beta^{*}=0$. When all predictors have unit norm, $\|X_{i}\|_{2}=1$, for $i=1,\ldots,p$, this test has the form

Remarkably, this above result is exact (nonasymptotic), valid for any $n$ and $p$, requiring (essentially) only Gaussianity of the errors, and no real assumptions about the matrix $X$. For most reasonably behaved predictor matrices $X$, the $\operatorname{Exp}(1)$ approximation agrees closely with this test. Details are in Taylor, Loftus and Tibshirani 2013. Work to extend this formula to subsequent steps along the solution path, that is, to test a hypothesis beyond the global null, is underway.

Generalizations to other penalties and models. The manuscript of Taylor, Loftus and Tibshirani 2013 applies to a regularized regression setting with a general seminorm penalty, and derives explicit results for the group lasso and nuclear norm penalties (in addition to the lasso penalty). The nuclear norm result yields a test for principal components and matrix completion. The recent work of Grazier G’Sell, Taylor and Tibshirani 2013 studies the covariance test for graphical models, based on a sparse estimate of the inverse covariance matrix.

Sequential procedures with false discovery rate control. It is also interesting to consider how the sequence of covariance test $p$-values can be used to construct a sequential test with good power properties, and a guaranteed bound on its false discovery rate. A number of such approaches are proposed in Grazier G’Sell et al. 2013.

Proper $p$-values for forward stepwise. Perhaps surprisingly, a test analogous to the covariance test can be used in forward stepwise regression, to construct valid $p$-values for this greedy procedure. This work is in progress.

Other related problems include: estimation of $\sigma^{2}$ when $p\geq n$, in the context of the covariance test; power calculations and confidence interval estimation; theory for linear models having strong and weak signals (large and small true coefficients); theory for the elastic net, generalized linear models, and the Cox model.

As is clear from the above discussion, the covariance test work has created much excitement and activity among our close collaborators and students. It is our hope that the current paper will also broadly stimulate other researchers’ interest in this area, and that at some point, the joint efforts of the community will yield a full set of inferential tools for the lasso and other commonly used adaptive procedures.

Appendix

By continuity of the lasso solution path at $\lambda_{k}$,

From this, we can obtain two identities: the first is

obtained by squaring both sides in (56) (more precisely, taking the inner product of the left-hand side with itself and the right-hand side with itself), and noting that $(P_{A\cup\{j\}}-P_{A})^{2}=P_{A\cup\{j\}}-P_{A}$; the second is

obtained by taking the inner product of both sides in (56) with $y$, and then using (57). Plugging (57) and (.1) in for the first and second terms in (2.3), respectively, then gives the result in (9).

.2 Proof of Lemma 4.

the first step following since $1-ss_{1}R_{j,j_{1}}>0$, and the second step following from the definition of $h^{(j_{1},s_{1})}$. The intersection of the right-hand side above, over all $(j,s)\neq(j_{1},s_{1})$, is equivalent to

But the former inequality is the same as $g(j_{1},s_{1})>0$, because $g(j_{1},s_{1})$ and $g(j_{1},-s_{1})$ have opposite signs. Further, the inequality $g(j_{1},s_{1})>0$ is redundant, as $M(j_{1},s_{1})\geq 0$. This gives the result.

.3 Proof of Lemma 5.

where $\phi$ is the standard normal density. First, note that

Also, a straightforward calculation shows

where in the last step we used the fact that $(1-\sqrt{1+4t/m^{2}})/(2/m^{2})\rightarrow-t/2$, again by l’Hôpital’s rule. Therefore, $\phi(u(t,m))/\phi(m)\rightarrow e^{-t}$, which completes the proof.

.4 Proof of Lemma 6.

Fix $\varepsilon>0$, and choose $m_{0}$ large enough that

Above, the term multiplying $\varepsilon$ is equal to 1, and the second term can be made arbitrarily small (say, less than $\varepsilon$) by taking $p$ sufficiently large.

.5 Proof of Theorem 2.

We will show that for any fixed $m_{0}>0$ and $j_{1},s_{1}$,

where $S\subseteq\{1,\ldots,p\}\setminus\{j_{1}\}$ is as in the theorem for $j=j_{1}$, with size $|S|\geq d_{p}$, and $c<1$ is a constant (not depending on $j_{1}$). This would imply that

where we used the fact that $d_{p}/\log{p}\rightarrow\infty$ by (27). The above sum tending to zero now implies the desired convergence result by Lemma 6, and hence it suffices to show (59). To this end, consider

where in both inequalities above we used the fact that $|R_{j,j_{1}}|<1$. We can therefore use the bound

where we define $V_{j}=(U_{j}-R_{j,j_{1}}U_{j_{1}})/2$ for $j\in S$. Let $r=|S|$, and without a loss of generality, let $S=\{1,\ldots,r\}$. We will show that