Exact post-selection inference, with application to the lasso

Introduction

As a statistical technique, linear regression is both simple and powerful. Not only does it provide estimates of the “effect” of each variable, but it also quantifies the uncertainty in those estimates, allowing inferences to be made about the effects. However, in many applications, a practitioner starts with a large pool of candidate variables, such as genes or demographic features, and does not know a priori which are relevant. This is especially problematic when there are more variables than observations, since then the model is unidentifiable (at least in the setting where the predictors are assumed fixed).

In such settings, it is tempting to let the data decide which variables to include in the model. For example, one common approach when the number of variables is not too large is to fit a linear model with all variables included, observe which ones are significant at level $\alpha$, and then refit the linear model with only those variables included. The problem with this is that the $p$-values can no longer be trusted, since the variables that are selected will tend to be those that are significant. Intuitively, we are “overfitting” to a particular realization of the data.

where $X_{M}^{+}\equiv(X_{M}^{T}X_{M})^{-1}X_{M}^{T}$ is the pseudo-inverse of $X_{M}$. Notice that (2) implies that the targets $\beta^{M}_{j}$ and $\beta^{M^{\prime}}_{j}$ in different models $M\neq M^{\prime}$ are in general different. This is simply a restatement of the well-known fact that a regression coefficient describes the effect of a predictor, adjusting for the other predictors in the model. In general, the coefficient of a predictor cannot be compared across different models.

Thus, “inference after selection” is ambiguous in linear regression because the target of inference changes with the selected model [Berk et al. (2013)]. In the next section, we discuss several ways to resolve this ambiguity.

Post-selection inference in linear regression

At first blush, the fact that the target ${\bolds\beta}^{M}$ changes with the model is deeply troubling, since it seems to imply that the parameters are random. However, the randomness is actually in the choice of which parameters to consider, not in the parameters themselves. Imagine that there are a priori $p2^{p-1}$ well-defined population parameters, one for each coefficient in all $2^{p}$ possible models:

We only ever form inferences for the parameters $\beta^{\hat{M}}_{j}$ in the model $\hat{M}$ we select. This adaptive choice of which parameters to consider can lead to inferences with undesirable frequency properties, as noted by Benjamini and Yekutieli (2005) and Benjamini, Heller and Yekutieli (2009).

To be concrete, suppose we want a confidence interval $C^{\hat{M}}_{j}$ for a parameter $\beta^{\hat{M}}_{j}$. What frequency properties should $C^{\hat{M}}_{j}$ have? By analogy to the classical setting, we might require that

but the event inside the probability is not well-defined because $\beta^{M}_{j}$ is undefined when $j\notin M$. Two ways around this issue are suggested by Berk et al. (2013): {longlist}[2.]

Conditional coverage: Since we form an interval for $\beta^{M}_{j}$ if and only if model $M$ is selected, that is, $\hat{M}=M$, it makes sense to condition on this event. Hence, we might require that our confidence interval $C^{M}_{j}$ satisfy

The benefit of this approach is that we avoid ever having to compare coefficients across two different models $M\neq M^{\prime}$.

Another way to understand conditioning on the model is to consider data splitting [Cox (1975)], an approach to post-selection inference that most statisticians would agree is valid. In data splitting, the data is divided into two halves, with one half used to select the model and the other used to conduct inference. Fithian, Sun and Taylor (2014) argues that inferences obtained by data splitting are only valid conditional on the model that was selected on the first half of the data. Therefore, conditional coverage is a reasonable frequency property to require of a post-selection confidence interval.

Simultaneous coverage: It also makes sense to talk about events that are defined simultaneously over all $j\in\hat{M}$. Berk et al. (2013) propose controlling the familywise error rate

but this is very stringent when many predictors are involved.

Instead of controlling the probability of making any error, we can control the expected proportion of errors—although “proportion of errors” is ambiguous in the event that we select zero variables. Benjamini and Yekutieli (2005) simply declare the error to be zero when $|\hat{M}|=0$:

while Storey (2003) suggests conditioning on $|\hat{M}|>0$:

Consider a family of intervals $\{C^{\hat{M}}_{j}\}_{j\in\hat{M}}$ that each have conditional $(1-\alpha)$ coverage:

Condition on $\hat{M}$ and iterate expectations:

Theorem 2 in Weinstein, Fithian and Benjamini (2013) proves a special case of Lemma 2.1 for a particular selection procedure, and Proposition 11 in Fithian, Sun and Taylor (2014) provides a more general result, but this result is sufficient for our purposes: to establish that conditional coverage is a sensible criterion to consider in post-selection inference.

Although the criterion is easy to state, how do we construct an interval with conditional coverage? This requires that we understand the conditional distribution

One of the main contributions of this paper is to show that this distribution is indeed possible to characterize, making valid post-selection inference feasible in the context of linear regression.

Outline of our approach

We have argued that post-selection intervals for regression coefficients should have $1-\alpha$ coverage conditional on the selected model:

both because this criterion is interesting in its own right and because it implies FCR control. To obtain an interval with this property, we study the conditional distribution

which will allow, more generally, conditional inference for parameters of the form ${\bolds\eta}_{M}^{T}{\bolds\mu}$. In particular, the regression coefficients $\beta^{M}_{j}={\mathbf{e}}_{j}^{T}X_{M}^{+}{\bolds\mu}$ can be written in this form, as can many other linear contrasts.

Our paper focuses on the specific case where the lasso is used to select the model $\hat{M}$. We begin in Section 4 by characterizing the event $\{\hat{M}=M\}$ for the lasso. As it turns out, this event is a union of polyhedra. More precisely, the event $\{\hat{M}=M,\hat{\mathbf{s}}_{M}={\mathbf{s}}_{M}\}$, that specifies the model and the signs of the selected variables, is a polyhedron of the form

Therefore, if we condition on both the model and the signs, then we only need to study

We do this in Section 5. It turns out that this conditional distribution is essentially a (univariate) truncated Gaussian. We use this to derive a statistic $F^{\mathbf{z}}({\bolds\eta}^{T}{\mathbf{y}})$ whose distribution given $\{A{\mathbf{y}}\leq{\mathbf{b}}\}$ is $\operatorname{Unif}(0,1)$.

The resulting post-selection test has a similar structure to the pathwise significance tests of Lockhart et al. (2014) and Taylor et al. (2014), which also are conditional tests. However, the intended application of our test is different. While their significance tests are specifically intended for the path context, our framework allows more general questions about the model the lasso selects: we can test the model at any value of $\lambda$ or form confidence intervals for an individual coefficient in the model.

There is also a parallel literature on confidence intervals for coefficients in high-dimensional linear models based on the lasso estimator [van de Geer et al. (2013); Zhang and Zhang (2014); Javanmard and Montanari (2013)]. The difference between their work and ours is that they do not address post-selection inference; their target is ${\bolds\beta}^{0}$, the coefficients in the true model, rather than ${\bolds\beta}^{\hat{M}}$, the coefficients in the selected model. The two will not be the same unless $\hat{M}$ happens to contain all nonzero coefficients of ${\bolds\beta}^{0}$. Although inference for ${\bolds\beta}^{0}$ is appealing, it requires assumptions about correctness of the linear model and sparsity of ${\bolds\beta}^{0}$. Pötscher and Schneider (2010) consider confidence intervals for the hard-thresholding and soft-thresholding estimators in the case of orthogonal design. Our approach instead regards the selected model as a linear approximation to the truth, a view shared by Berk et al. (2013) and Miller (2002).

The idea of post-selection inference conditional on the selected model appears in Pötscher (1991), although the notion of inference conditional on certain relevant subsets dates back to Fisher (1956); see also Robinson (1979). Leeb and Pötscher (2005; 2006) obtained a number of negative results about estimating the distribution of a post-selection estimator, although they note their results do not necessarily preclude the possibility of post-selection inference. Benjamini and Yekutieli (2005) also consider conditioning on the selection event, although they argue that this is too conservative. To the contrary, we show that conditioning on the selected model can produce reasonable confidence intervals in a wide variety of situations.

Inference conditional on selection has also appeared in literature on the winner’s curse: Sampson and Sill (2005); Sill and Sampson (2009); Zhong and Prentice (2008); Zollner and Pritchard (2007). These works are not really associated with model selection in linear regression, though they employ a similar approach to inference.

The lasso and its selection event

Because the lasso produces sparse solutions, we can define model “selected” by the lasso to be simply the set of predictors with nonzero coefficients:

Then post-selection inference seeks to make inferences about ${\bolds\beta}^{M}$, given $\{\hat{M}=M\}$, as defined in (2).

The rest of this section focuses on characterizing this event $\{\hat{M}=M\}$. We begin by noting that in order for a vector of coefficients $\hat{\bolds\beta}$ and a vector of signs $\hat{\mathbf{s}}$ to be solutions to the lasso problem (9), it is necessary and sufficient that they satisfy the Karush–Kuhn–Tucker (KKT) conditions:

Following Tibshirani (2013), we consider the equicorrelation set

Notice that we have implicitly identified the model $\hat{M}$ with the equicorrelation set. Since $|\hat{s}_{i}|=1$ for any $\hat{\beta}_{i}\neq 0$, the equicorrelation set does in fact contain all predictors with nonzero coefficients, although it may also include some predictors with zero coefficients. However, for almost every $\lambda$, the equicorrelation set is precisely the set of predictors with nonzero coefficients.

It turns out that it is easier to first characterize $\{(\hat{M},\hat{\mathbf{s}})=(M,{\mathbf{s}})\}$ and obtain $\{\hat{M}=M\}$ as a corollary by taking a union over the possible signs. The next result is an important first step.

Assume the columns of $X$ are in general position [Tibshirani (2013)]. Let $M\subset\{1,\dots,p\}$ and ${\mathbf{s}}\in\{-1,1\}^{|M|}$ be a candidate set of variables and their signs, respectively. Define the random variables

where $P_{M}\equiv X_{M}(X_{M}^{T}X_{M})^{-1}X_{M}$ is projection onto the column span of $X_{M}$. Then the selection procedure can be rewritten in terms of ${\mathbf{w}}$ and ${\mathbf{u}}$ as

First, we rewrite the KKT conditions (10) by partitioning them according to the equicorrelation set $\hat{M}$, adopting the convention that $-\hat{M}$ means “variables not in $\hat{M}$”:

Since the KKT conditions are necessary and sufficient for a solution, we obtain that ${\{(\hat{M},\hat{\mathbf{s}})=(M,{\mathbf{s}})\}}$ if and only if there exist ${\mathbf{w}}$ and ${\mathbf{u}}$ satisfying

We can solve the first two equations for ${\mathbf{w}}$ and ${\mathbf{u}}$ to obtain the equivalent set of conditions

where the first two are the definitions of ${\mathbf{w}}$ and ${\mathbf{u}}$ given in (13) and (14), and the last two are the conditions on ${\mathbf{w}}$ and ${\mathbf{u}}$ given in (15).

Lemma 4.1 is remarkable because it says that the event ${\{(\hat{M},\hat{\mathbf{s}})=(M,{\mathbf{s}})\}}$ can be rewritten as affine constraints on ${\mathbf{y}}$. This is because ${\mathbf{w}}$ and ${\mathbf{u}}$ are already affine functions of ${\mathbf{y}}$, and the constraints $\operatorname{sign}(\cdot)={\mathbf{s}}$ and $\|\cdot\|_{\infty}<1$ can also be rewritten in terms of affine constraints. The following proposition makes this explicit.

Let ${\mathbf{w}}$ and ${\mathbf{u}}$ be defined as in (13) and (14). Then

where $A_{0},{\mathbf{b}}_{0}$ encode the “inactive” constraints $\{\|{\mathbf{u}}\|_{\infty}<1\}$, and $A_{1},{\mathbf{b}}_{1}$ encode the “active” constraints $\{\operatorname{sign}({\mathbf{w}})={\mathbf{s}}\}$. These matrices have the explicit forms

First, substituting expression (13) for ${\mathbf{w}}$, we rewrite the “active” constraints as

Next, substituting expression (14) for ${\mathbf{u}}$, we rewrite the “inactive” constraints as

Combining Lemma 4.1 with Proposition 4.2, we obtain the following.

Let $A(M,{\mathbf{s}})={A_{0}(M,{\mathbf{s}})\choose A_{1}(M,{\mathbf{s}})}$ and $b(M,{\mathbf{s}})={{\mathbf{b}}_{0}(M,{\mathbf{s}})\choose{\mathbf{b}}_{1}(M,{\mathbf{s}})}$, where $A_{i}$ and $b_{i}$ are defined in Proposition 4.2. Then

As a corollary, $\{\hat{M}=M\}$ is simply the union of the above events over all possible sign patterns.

$\{\hat{M}=M\}=\bigcup_{{\mathbf{s}}\in\{-1,1\}^{|M|}}\{A(M,{\mathbf{s}}){\mathbf{y}}\leq{\mathbf{b}}(M,{\mathbf{s}})\}$.

Polyhedral conditioning sets

In order to obtain inference conditional on the model, we need to understand the distribution of

However, as we saw in the previous section, $\{\hat{M}=M\}$ is a union of polyhedra, so it is easier to condition on both the model and the signs,

since the conditioning event is a single polyhedron $\{A(M,{\mathbf{s}}){\mathbf{y}}\leq{\mathbf{b}}(M,{\mathbf{s}})\}$. Notice that inferences that are valid conditional on this finer event will also be valid conditional on $\{\hat{M}=M\}$. For example, if a confidence interval $C^{M}_{j}$ for $\beta^{M}_{j}$ has $(1-\alpha)$ coverage conditional on the model and signs

it will also have $(1-\alpha)$ coverage conditional only on the model by the Law of Total Probability:

This section is divided into two subsections. First, we study how to condition on a single polyhedron; this will allow us to condition on $\{\hat{M}=M,\hat{\mathbf{s}}={\mathbf{s}}\}$. Then we extend the framework to condition on a union of polyhedra, which will allow us to condition only on the model $\{\hat{M}=M\}$. The inferences obtained by conditioning on the model will in general be more efficient (i.e., narrower intervals, more powerful tests), at the price of more computation.

we rewrite $\{A{\mathbf{y}}\leq{\mathbf{b}}\}$ in terms of ${\bolds\eta}^{T}{\mathbf{y}}$ and a component ${\mathbf{z}}$ which is independent of ${\bolds\eta}^{T}{\mathbf{y}}$. That component is

It is easy to verify that ${\mathbf{z}}$ is uncorrelated with, and hence independent of, ${\bolds\eta}^{T}{\mathbf{y}}$. Notice that in the case where $\Sigma=\sigma^{2}I_{n}$, ${\mathbf{z}}$ is simply the residual $(I_{n}-P_{{\bolds\eta}}){\mathbf{y}}$ from projecting ${\mathbf{y}}$ onto ${\bolds\eta}$.

We can now rewrite $\{A{\mathbf{y}}\leq{\mathbf{b}}\}$ in terms of ${\bolds\eta}^{T}{\mathbf{y}}$ and ${\mathbf{z}}$.

Let ${\mathbf{z}}$ be defined as in (18) and ${\mathbf{c}}$ as in (19). Then the conditioning set can be rewritten as follows:

Note that ${\mathcal{V}}^{-}$, ${\mathcal{V}}^{+}$, and ${\mathcal{V}}^{0}$ refer to functions. Since they are functions of ${\mathbf{z}}$ only, (20)–(22) are independent of ${\bolds\eta}^{T}{\mathbf{y}}$.

We can decompose ${\mathbf{y}}={\mathbf{c}}({\bolds\eta}^{T}{\mathbf{y}})+{\mathbf{z}}$ and rewrite the polyhedron as

where in the last step, we have divided the components into three categories depending on whether $(A{\mathbf{c}})_{j}\gtreqless 0$, since this affects the direction of the inequality (or whether we can divide at all). Since ${\bolds\eta}^{T}{\mathbf{y}}$ is the same quantity for all $j$, it must be at least the maximum of the lower bounds, which is ${\mathcal{V}}^{-}({\mathbf{z}})$, and no more than the minimum of the upper bounds, which is ${\mathcal{V}}^{+}({\mathbf{z}})$.

Since ${\mathcal{V}}^{+}({\mathbf{z}}),{\mathcal{V}}^{-}({\mathbf{z}}),{\mathcal{V}}^{0}({\mathbf{z}})$ are independent of ${\bolds\eta}^{T}{\mathbf{y}}$, they behave as “fixed” quantities. Thus, ${\bolds\eta}^{T}{\mathbf{y}}$ is conditionally like a normal random variable, truncated to be between ${\mathcal{V}}^{-}({\mathbf{z}})$ and ${\mathcal{V}}^{+}({\mathbf{z}})$. We would like to be able to say

but this is technically incorrect, since the distribution on the right-hand side changes with ${\mathbf{z}}$. By conditioning on the value of ${\mathbf{z}}$, ${\bolds\eta}^{T}{\mathbf{y}}|\{A{\mathbf{y}}\leq{\mathbf{b}},{\mathbf{z}}={\mathbf{z}}_{0}\}$ is a truncated normal. We can then use the probability integral transform to obtain a statistic $F^{\mathbf{z}}({\bolds\eta}^{T}{\mathbf{y}})$ that has a $\operatorname{Unif}(0,1)$ distribution for any value of ${\mathbf{z}}$. Hence, $F^{\mathbf{z}}({\bolds\eta}^{T}{\mathbf{y}})$ will also have a $\operatorname{Unif}(0,1)$ distribution marginally over ${\mathbf{z}}$. We make this precise in the next theorem.

Let $F_{\mu,\sigma^{2}}^{[a,b]}$ denote the CDF of a $N(\mu,\sigma^{2})$ random variable truncated to the interval $[a,b]$, that is,

where $\Phi$ is the CDF of a $N(0,1)$ random variable. Then

where ${\mathcal{V}}^{-}$ and ${\mathcal{V}}^{+}$ are defined in (20) and (21). Furthermore,

The only random quantities left are ${\bolds\eta}^{T}{\mathbf{y}}$ and ${\mathbf{z}}$. Now we can eliminate ${\mathbf{z}}={\mathbf{z}}_{0}$ from the condition using independence:

Letting $F^{\mathbf{z}}({\bolds\eta}^{T}{\mathbf{y}})\equiv F_{{\bolds\eta}^{T}{\bolds\mu},{\bolds\eta}^{T}\Sigma{\bolds\eta}}^{[{\mathcal{V}}^{-}({\mathbf{z}}),{\mathcal{V}}^{+}({\mathbf{z}})]}({\bolds\eta}^{T}{\mathbf{y}})$, we can apply the probability integral transform to the above result to obtain

If we let $p_{X}$ denote the density of a random variable $X$ given $\{A{\mathbf{y}}\leq{\mathbf{b}}\}$, what we have just shown is that

for any ${\mathbf{z}}_{0}$. The desired result now follows by integrating over ${\mathbf{z}}_{0}$:

2 Conditioning on a union of polyhedra

We have just characterized the distribution of ${\bolds\eta}^{T}{\mathbf{y}}$, conditional on ${\mathbf{y}}$ falling into a single polyhedron $\{A{\mathbf{y}}\leq{\mathbf{b}}\}$. We obtain such a polyhedron if we condition on both the model and the signs $\{\hat{M}=M,\hat{\mathbf{s}}={\mathbf{s}}\}$. If we want to only condition on the model $\{\hat{M}=M\}$, then we will have to understand the distribution of ${\bolds\eta}^{T}{\mathbf{y}}$, conditional on ${\mathbf{y}}$ falling into a union of such polyhedra, that is,

As Figure 3 makes clear, the argument proceeds exactly as before, except that ${\bolds\eta}^{T}{\mathbf{y}}$ is now truncated to a union of intervals, instead of a single interval. There is a ${\mathcal{V}}^{-}$ and a ${\mathcal{V}}^{+}$ for each possible sign pattern ${\mathbf{s}}$, so we index the intervals by the signs. This leads immediately to the next theorem, whose proof is essentially the same as that of Theorem 5.2.

Let $F_{\mu,\sigma^{2}}^{S}$ denote the CDF of a $N(\mu,\sigma^{2})$ random variable truncated to the set $S$. Then

where ${\mathcal{V}}_{{\mathbf{s}}}^{-}({\mathbf{z}})$ and ${\mathcal{V}}_{{\mathbf{s}}}^{+}({\mathbf{z}})$ are defined in (20) and (21) and $A=A_{\mathbf{s}}$ and $b=b_{\mathbf{s}}$.

Post-selection intervals for regression coefficients

In this section, we combine the characterization of the lasso selection event in Section 4 with the results about the distribution of a Gaussian truncated to a polyhedron (or union of polyhedra) in Section 5 to form post-selection intervals for lasso-selected regression coefficients. The key link is that the lasso selection event can be expressed as a union of polyhedra:

where $A(M,{\mathbf{s}})$ and ${\mathbf{b}}(M,{\mathbf{s}})$ are defined in Theorem 4.3. Therefore, conditioning on selection is the same as conditioning on a union of polyhedra, so we can apply the framework of Section 5.

Recall that our goal is to form confidence intervals for $\beta^{M}_{j}={\mathbf{e}}_{j}^{T}X_{M}^{+}{\bolds\mu}$, with $(1-\alpha)$-coverage conditional on $\{\hat{M}=M\}$. Taking ${\bolds\eta}=(X_{M}^{+})^{T}{\mathbf{e}}_{j}$, we can use Theorem 5.3 to obtain

This gives us a test statistic for testing any hypothesized value of ${\bolds\beta}^{M}_{j}$. We can invert this test to obtain a confidence set

In fact, the set $C^{M}_{j}$ is an interval, as formalized in the next result.

Let ${\bolds\eta}=(X_{M}^{+})^{T}{\mathbf{e}}_{j}$. Let $L$ and $U$ be the (unique) values satisfying

Then $[L,U]$ is a $(1-\alpha)$ confidence interval for $\beta^{M}_{j}$, conditional on $\{\hat{M}=M\}$, that is,

Alternatively, we could have conditioned on the signs, in addition to the model, so that we would only have to condition on a single polyhedron. We also showed in Section 5 that

Inverting this statistic will produce intervals that have $(1-\alpha)$ coverage conditional on $\{\hat{M}=M,\hat{\mathbf{s}}={\mathbf{s}}\}$, and hence $(1-\alpha)$ coverage conditional on $\{\hat{M}=M\}$. However, these intervals will be less efficient; they will in general be wider. However, one may be willing to sacrifice statistical efficiency for computational efficiency. Notice that the main cost in computing intervals according to Theorem 6.1 is determining the intervals $[{\mathcal{V}}^{-}_{\mathbf{s}}({\mathbf{z}}),{\mathcal{V}}^{+}_{\mathbf{s}}({\mathbf{z}})]$ for each ${\mathbf{s}}\in\{-1,1\}^{|M|}$. The number of such sign patterns is $2^{|M|}$. While this might be feasible when $|M|$ is small, it is not feasible when we select hundreds of variables. Conditioning on the signs means that we only have to compute the interval $[{\mathcal{V}}^{-}_{\mathbf{s}}({\mathbf{z}}),{\mathcal{V}}^{+}_{\mathbf{s}}({\mathbf{z}})]$ for the sign pattern ${\mathbf{s}}$ that was actually observed.

Figure 4 shows the tradeoff in statistical efficiency. When the signal is strong, as in the left-hand plot, there is virtually no difference between the intervals obtained by conditioning on just the model, or the model and signs. On the other hand, in the right-hand plot, we see that we can obtain very wide intervals when the signal is weak. The widest intervals are for actual noise variables, as expected.

To understand why post-selection intervals are sometimes very wide, notice that when a truncated Gaussian random variable $Z$ is close to the endpoints of the truncation interval $[a,b]$, there are many means $\mu$ that would be consistent with that observation—hence, the wide intervals. Figure 5 shows confidence intervals for $\mu$ as a function of $Z$. When $Z$ is far from the endpoints of the truncation interval, we basically recover the nominal OLS intervals (i.e., not adjusted for selection).

The implications are clear. When the signal is strong, ${\bolds\eta}^{T}{\mathbf{y}}$ will be far from the endpoints of the truncation region, so we obtain the nominal OLS intervals. On the other hand, when a variable just barely entered the model, then ${\bolds\eta}^{T}{\mathbf{y}}$ will be close to the edge of the truncation region, and the interval will be wide.

We have derived a confidence interval $C^{M}_{j}$ whose conditional coverage, given $\{\hat{M}=M\}$, is at least $1-\alpha$. The fact that we have found such an interval is not remarkable, since many such intervals have this property. However, given two intervals with the same coverage, we generally prefer the shorter one. This problem is considered in Fithian, Sun and Taylor (2014) where it is shown that $C^{M}_{j}$ is, with one small tweak, the shortest interval among all unbiased intervals with $1-\alpha$ coverage.

An unbiased interval $C$ for a parameter $\theta$ is one which covers no other parameter $\theta^{\prime}$ with probability more than $1-\alpha$, that i,

Unbiasedness is a common restriction to ensure the existence of an optimal interval [Lehmann and Romano (2005)]. The shortest unbiased interval for $\beta^{M}_{j}$, among all intervals with conditional $1-\alpha$ coverage, resembles to the interval $[L,U]$ in Theorem 6.1. There, the critical values $L$ and $U$ were chosen symmetrically so that the pivot has $\alpha/2$ area in either tail. However, it may be possible to obtain a shorter interval on average by allocating the a probability unequally between the two tails. Theorem 5 of Fithian, Sun and Taylor (2014) provides a general formula for obtaining shortest unbiased intervals in exponential families.

Data example

We apply our post-selection intervals to the diabetes data set from Efron et al. (2004). Since $p<n$ for this data set, we can estimate $\sigma^{2}$ using the residual sum of squares from the full regression model with all $p$ predictors. After standardizing all variables, we chose $\lambda$ according to the strategy in Negahban et al. (2012), $\lambda=2\mathbf{E}(\|X^{T}\varepsilon\|_{\infty})$. This expectation was computed by simulation, where $\varepsilon\sim N(0,\hat{\sigma}^{2})$, resulting in $\lambda\approx 190$. The lasso selected four variables: BMI, BP, S3 and S5.

The post-selection intervals are shown in Figure 6, alongside the nominal confidence intervals produced by fitting OLS to the four selected variables, ignoring selection. The nominal intervals do not have $(1-\alpha)$ coverage conditional on the model and are not valid post-selection intervals. Also depicted are the confidence intervals obtained by data splitting, as discussed in Section 2. This is a competitor method that also produces valid confidence intervals conditional on the model. The lasso selected the same four variables on half of the data, and then nominal intervals for these four variables using OLS on the other half of the data.

We can make two observations from Figure 6. {longlist}[2.]

The adjusted intervals provided by our method essentially reproduces the OLS intervals for the strong effects, whereas data splitting intervals are wider by a factor of $\sqrt{2}$ (since only $n/2$ observations are used in the inference). For this dataset, the POSI intervals are $1.36$ times wider than the OLS intervals. For all the variables, our method produces the shortest intervals among the methods that control selective type 1 error.

One variable, S3 which would have been deemed significant using the OLS intervals, is no longer significant after accounting for selection. Data splitting, our selection-adjusted intervals, and POSI intervals conclude that S3 is not significant. This demonstrates that taking model selection into account can have substantive impacts on the conclusions.

Extensions

The above results rely on knowing $\sigma^{2}$ or at least having a good estimate of it. If $n>p$, then the variance $\hat{\sigma}^{2}$ of the residuals from fitting the full model is a consistent estimator and in general can be substituted for $\sigma^{2}$ to yield asymptotically valid confidence intervals. Formally, the condition is that the pivot is smooth with respect to $\sigma$. Geometrically speaking, the upper and lower truncation limits ${\mathcal{V}}^{+}$ and ${\mathcal{V}}^{-}$ must be well-separated (with high probability). We refer the interested reader to Section 2.3 in Tian and Taylor (2015) for details.

In the setting where $p>n$, obtaining an estimate of $\sigma^{2}$ is more challenging, but if the pivot satisfies a monotonicity property, plugging in an overestimate of the variance gives conservative confidence intervals. We refer the reader to Theorem 11 in Tibshirani et al. (2015) for details.

2 Elastic net

These four conditions differ from those of Lemma 4.1 in only one respect: $X_{M}^{T}X_{M}$ in the first expression is replaced by $X_{M}^{T}X_{M}+\gamma I$. Continuing the argument of Section 4, we see that the selection event can be rewritten

Now that we have rewritten the selection event in the form $\{A{\mathbf{y}}\leq{\mathbf{b}}\}$, we can once again apply the framework in Section 5 to obtain a test for the elastic net conditional on this event.

Conclusion

Model selection and inference have long been regarded as conflicting goals in linear regression. Following the lead of Berk et al. (2013), we have proposed a framework for post-selection inference that conditions on which model was selected, that is, the event $\{\hat{M}=M\}$. We characterize this event for the lasso and derive optimal and exact confidence intervals for linear contrasts ${\bolds\eta}^{T}{\bolds\mu}$, conditional on $\{\hat{M}=M\}$. With this general framework, we can form post-selection intervals for regression coefficients, equipping practitioners with a way to obtain “valid” intervals even after model selection.

Appendix: Monotonicity of F𝐹F

Let $F_{\mu}(x):=F_{\mu,\sigma^{2}}^{[a,b]}(x)$ denote the cumulative distribution function of a truncated Gaussian random variable, as defined as in (24). Then $F_{\mu}(x)$ is monotone decreasing in $\mu$.

First, the truncated Gaussian distribution with CDF $F_{\mu}:=F_{\mu,\sigma^{2}}^{[a,b]}$ is a natural exponential family in $\mu$, since it is just a Gaussian with a different base measure. Therefore, it has monotone likelihood ratio in $\mu$. That is, for all $\mu_{1}>\mu_{0}$ and $x_{1}>x_{0}$:

where $f_{\mu_{i}}:=dF_{\mu_{i}}$ denotes the density. (Instead of appealing to properties of exponential families, this property can also be directly verified.)

Therefore, the inequality is preserved if we integrate both sides with respect to $x_{0}$ on $(-\infty,x)$ for $x<x_{1}$. This yields

Now we integrate both sides with respect to $x_{1}$ on $(x,\infty)$ to obtain

which establishes $F_{\mu_{0}}(x)>F_{\mu_{1}}(x)$ for all $\mu_{1}>\mu_{0}$.

Acknowledgements

We thank Will Fithian, Sam Gross and Josh Loftus for helpful comments and discussions. In particular, Will Fithian provided insights that led to the geometric intuition of our procedure shown in Figure 2.

References