Exact Post Model Selection Inference for Marginal Screening

Jason D Lee, Jonathan E Taylor

Introduction

where μ(x)\mu(x) is an arbitrary function, and xi∈Rpx_{i}\in\mathbf{R}^{p}. Our goal is to perform inference on (XTX)−1XTμ(X^{T}X)^{-1}X^{T}\mu, which is the best linear predictor of μ\mu. In the classical setting of n>pn>p , the least squares estimator

is a commonly used estimator for (XTX)−1XTμ(X^{T}X)^{-1}X^{T}\mu. Under the linear model assumption μ=Xβ0\mu=X\beta^{0}, the exact distribution of β^\hat{\beta} is

Using the normal distribution, we can test the hypothesis H0:βj0=0H_{0}:\beta^{0}_{j}=0 and form confidence intervals for βj0\beta^{0}_{j} using the z-test.

However in the high-dimensional p>np>n setting, the least squares estimator is an underdetermined problem, and the predominant approach is to perform variable selection or model selection . There are many approaches to variable selection including AIC/BIC, greedy algorithms such as forward stepwise regression, orthogonal matching pursuit, and regularization methods such as the Lasso. The focus of this paper will be on the model selection procedure known as marginal screening, which selects the kk most correlated features xjx_{j} with the response yy.

Marginal screening is the simplest and most commonly used of the variable selection procedures . Marginal screening requires only O(np)O(np) computation and is several orders of magnitude faster than regularization methods such as the Lasso; it is extremely suitable for extremely large datasets where the Lasso may be computationally intractable to apply. Furthermore, the selection properties are comparable to the Lasso . In the ultrahigh dimensional setting p=O(enk)p=O(e^{n^{k}}), marginal screening is shown to have the SURE screening property, P(S⊂S^P(S\subset\hat{S}), that is marginal screening selects a superset of the truly relevant variables . Marginal screening can also be combined with a second variable selection procedure such as the Lasso to further reduce the dimensionality; our statistical inference methods extend to the Marginal Screening+Lasso method.

Since marginal screening utilizes the response variable yy, the confidence intervals and statistical tests based on the distribution in (3) are not valid; confidence intervals with nominal 1−α1-\alpha coverage may no longer cover at the advertised level:

Several authors have previously noted this problem including recent work in . A major line of work has described the difficulty of inference post model selection: the distribution of post model selection estimates is complicated and cannot be approximated in a uniform sense by their asymptotic counterparts.

In this paper, we describe how to form exact confidence intervals for linear regression coefficients post model selection. We assume the model (1), and operate under the fixed design matrix XX setting. The linear regression coefficients constrained to a subset of variables SS is linear in μ\mu, ejT(XSTXS)−1XSTμ=ηTμe_{j}^{T}(X_{S}^{T}X_{S})^{-1}X_{S}^{T}\mu=\eta^{T}\mu for some η\eta. We derive the conditional distribution of ηTy\eta^{T}y for any vector η\eta, so we are able to form confidence intervals and test regression coefficients.

In Section 2 we discuss related work on high-dimensional statistical inference, and Section 3 introduces the marginal screening algorithm and shows how z intervals may fail to have the correct coverage properties. Section 4 and 5 show how to represent the marginal screening selection event as constraints on yy, and construct pivotal quantities for the truncated Gaussian. Section 6 uses these tools to develop valid hypothesis tests and confidence intervals.

Although the focus of this paper is on marginal screening, the “condition on selection” framework, first proposed for the Lasso in , is much more general; we use marginal screening as a simple and clean illustration of the applicability of this framework. In Section 7, we discuss several extensions including how to apply the framework to other variable/model selection procedures and to nonlinear regression problems. Section 7 covers

marginal screening+Lasso, a screen and clean procedure that first uses marginal screening and cleans with the Lasso,

Related Work

Most of the theoretical work on high-dimensional linear models focuses on consistency. Such results establish, under restrictive assumptions on XX, the Lasso β^\hat{\beta} is close to the unknown β0\beta^{0} and selects the correct model . We refer to the reader to for a comprehensive discussion about the theoretical properties of the Lasso.

There is also recent work on obtaining confidence intervals and significance testing for penalized M-estimators such as the Lasso. One class of methods uses sample splitting or subsampling to obtain confidence intervals and p-values . In the post model selection literature, the recent work of proposed the POSI approach, a correction to the usual t-test confidence intervals by controlling the familywise error rate for all parameters in any possible submodel. The POSI approach will produce valid confidence intervals for any possible model selection procedure; however for a given model selection procedure such as marginal regression, it will be conservative. In addition, the POSI methodology is extremely computationally intensive and currently only applicable for p≤30p\leq 30.

A separate line of work establishes the asymptotic normality of a corrected estimator obtained by “inverting” the KKT conditions . The corrected estimator b^\hat{b} has the form b^=β^+λΘ^z^,\hat{b}=\hat{\beta}+\lambda\hat{\Theta}\hat{z}, where z^\hat{z} is a subgradient of the penalty at β^\hat{\beta} and Θ^\hat{\Theta} is an approximate inverse to the Gram matrix XTXX^{T}X. The two main drawbacks to this approach are 1) the confidence intervals are valid only when the M-estimator is consistent, and thus require restricted eigenvalue conditions on XX, 2) obtaining Θ^\hat{\Theta} is usually much more expensive than obtaining β^\hat{\beta}, and 3) the method is specific to regularized estimators, and does not extend to marginal screening, forward stepwise, and other variable selection methods.

Most closely related to our work is the “condition on selection” framework laid out in for the Lasso. Our work extends this methodology to other variable selection methods such as marginal screening, marginal screening followed by the Lasso (marginal screening+Lasso), orthogonal matching pursuit, and non-negative least squares. The primary contribution of this work is the observation that many model selection methods, including marginal screening and Lasso, lead to “selection events” that can be represented as a set of constraints on the response variable yy. By conditioning on the selection event, we can characterize the exact distribution of ηTy\eta^{T}y. This paper focuses on marginal screening, since it is the simplest of variable selection methods, and thus the applicability of the “conditioning on selection event” framework is most transparent. However, this framework is not limited to marginal screening and can be applied to a wide a class of model selection procedures including greedy algorithms such as matching pursuit and orthogonal matching pursuit. We discuss some of these possible extensions in Section 7, but leave a thorough investigation to future work.

A remarkable aspect of our work is that we only assume XX is in general position, and the test is exact, meaning the distributional results are true even under finite samples. By extension, we do not make any assumptions on nn and pp, which is unusual in high-dimensional statistics . Furthermore, the computational requirements of our test are negligible compared to computing the linear regression coefficients.

Our test assumes that the noise variance σ2\sigma^{2} is known. However, there are many methods for estimating σ2\sigma^{2} in high dimensions. A data splitting technique is used in , while proposes a method that computes the regression estimate and an estimate of the variance simultaneously. We refer the reader to for a survey and comparison of the various methods, and assume σ2\sigma^{2} is known for the remainder of the paper.

Marginal Screening

We will assume that XX is in general position and has unit norm columns. The algorithm estimates β^\hat{\beta} via Algorithm 1.

The marginal screening algorithm chooses the kk variables with highest absolute dot product with yy, and then fits a linear model over those kk variables. We will assume k≤min⁡(n,p)k\leq\min(n,p). For any fixed subset of variables SS, the distribution of β^S=(XSTXS)−1XSTy\hat{\beta}_{S}=(X_{S}^{T}X_{S})^{-1}X_{S}^{T}y is

We will use the notation βj∈S⋆:=(βS⋆)j\beta^{\star}_{j\in S}:=\left(\beta^{\star}_{S}\right)_{j}, where jj is indexing a variable in the set SS. The z-test intervals for a regression coefficient are

and each interval has 1−α1-\alpha coverage, meaning Pr⁡(βj∈S⋆∈C(α,j,S))=1−α\Pr\left(\beta^{\star}_{j\in S}\in C(\alpha,j,S)\right)=1-\alpha. However if S^\hat{S} is chosen using a model selection procedure that depends on yy, the distributional result (5) no longer holds and the z-test intervals will not cover at the 1−α1-\alpha level. It is possible that

Similarly, the test of the hypothesis H0:βj∈S^⋆=0H_{0}:\beta^{\star}_{j\in\hat{S}}=0 will not control type I error at level α\alpha, meaning \Pr\left(\text{rejectH_{0}}|H_{0}\right)>\alpha.

We will illustrate empirically that the z-test intervals do not cover at 1−α1-\alpha when S^\hat{S} is chosen by marginal screening in Algorithm 1.

For this experiment we generated XX from a standard normal with n=20n=20 and p=200p=200. The signal vector is 22 sparse with β10,β20=SNR\beta^{0}_{1},\beta^{0}_{2}=\text{SNR}, y=Xβ0+ϵy=X\beta^{0}+\epsilon, and ϵ∼N(0,1)\epsilon\sim N(0,1). The confidence intervals were constructed for the k=2k=2 variables selected by the marginal screening algorithm. The z-test intervals were constructed via (6) with α=.1\alpha=.1, and the adjusted intervals were constructed using Algorithm 3. The results are described in Figure 1. The y-axis plots the coverage proportion or the fraction of times the true parameter value fell in the confidence interval. Each point represents 500500 independent trials. The x-axis varies the SNR parameter over the values 0.1,.2,.5,1,2,5,100.1,.2,.5,1,2,5,10. From the figure, we see that the z intervals can have coverage proportion drastically less than the nominal level of 1−α=.91-\alpha=.9, and only for SNR=1010 does the coverage tend to .9.9. This motivates the need for intervals that have the correct coverage proportion after model selection.

Representing the selection event

Since Equation (5) does not hold for a selected S^\hat{S} when the selection procedure depends on yy, the z-test intervals are not valid. Our strategy will be to understand the conditional distribution of yy and contrasts (linear functions of yy) ηTy\eta^{T}y, then construct inference conditional on the selection event E^\hat{E}. We will use E^(y)\hat{E}(y) to represent a random variable, and EE to represent an element of the range of E^(y)\hat{E}(y). In the case of marginal screening, the selection event E^(y)\hat{E}(y) corresponds to the set of selected variables S^\hat{S} and signs ss:

for some matrix A(S^,s^)A(\hat{S},\hat{s}) and vector b(S^,s^)b(\hat{S},\hat{s})bb can be taken to be for marginal screening, but this extra generality is needed for other model selection methods. We will use the selection event E^\hat{E} and the selected variables/signs pair (S^,s^)(\hat{S},\hat{s}) interchangeably since they are in bijection.

The space Rn\mathbf{R}^{n} is partitioned by the selection events,

The vector yy can be decomposed with respect to the partition as follows

The previous equation establishes that yy is a different constrained Gaussian for each element of the partition, where the partition is specified by a possible subset of variables and signs (S,s)(S,s). The above discussion can be summarized in the following theorem.

The distribution of yy conditional on the selection event is a constrained Gaussian,

The event EE is in bijection with a pair (S,s)(S,s), and yy is unconditionally Gaussian. Thus the conditional y\big{|}\{A(S,s)y\leq b(S,s)\} is a Gaussian constrained to the set {A(S,s)y≤b(S,s)}\{A(S,s)y\leq b(S,s)\}. ∎

Truncated Gaussian test

This section summarizes the recent tools developed in for testing contrastsA contrast of yy is a linear function of the form ηTy\eta^{T}y. ηTy\eta^{T}y of a constrained Gaussian yy. The results are stated without proof and the proofs can be found in .

The conditioning set can be rewritten in terms of ηTy\eta^{T}y as follows:

Moreover, (V+,V−,V0)({\cal V}^{+},{\cal V}^{-},{\cal V}^{0}) are independent of ηTy\eta^{T}y.

The geometric picture gives more intuition as to why V+{\cal V}^{+} and V−{\cal V}^{-} are independent of ηTy\eta^{T}y. Without loss of generality, we assume ∣∣η∣∣2=1||\eta||_{2}=1 and y∼N(μ,I)y\sim N(\mu,I) (otherwise we could replace yy by Σ−12y\Sigma^{-\frac{1}{2}}y). Now we can decompose yy into two independent components, a 1-dimensional component ηTy\eta^{T}y and an (n−1)(n-1)-dimensional component orthogonal to η\eta:

The case of n=2n=2 is illustrated in Figure 2. Since the two components are independent, the distribution of ηTy\eta^{T}y is the same as ηTy∣{Pη⊥y}\eta^{T}y|\{P_{\eta^{\perp}}y\}. If we condition on Pη⊥yP_{\eta^{\perp}}y, it is clear from Figure 2 that in order for yy to lie in the set, it is necessary for V−≤ηTy≤V+{\cal V}^{-}\leq\eta^{T}y\leq{\cal V}^{+}, where V−{\cal V}^{-} and V+{\cal V}^{+} are functions of Pη⊥yP_{\eta^{\perp}}y.

The distribution of ηTy\eta^{T}y conditioned on {Ay≤b,V+(y)=v+,V−(y)=v−}\{Ay\leq b,{\cal V}^{+}(y)=v^{+},{\cal V}^{-}(y)=v^{-}\} is a (univariate) Gaussian truncated to fall between V−{\cal V}^{-} and V+{\cal V}^{+}, i.e.

where W∼TN(ηTμ,ηTΣη,v−,v+)W\sim TN(\eta^{T}\mu,\eta^{T}\Sigma\eta,v^{-},v^{+}). TN(μ,σ,a,b)TN(\mu,\sigma,a,b) is the normal distribution truncated to lie between aa and bb.

In Figure 3, we plot the density of the truncated Gaussian, noting that its shape depends on the location of μ\mu relative to [a,b][a,b] as well as the width relative to σ\sigma.

The following pivotal quantityThe distribution of a pivotal quantity does not depend on unobserved parameters. follows from Corollary 5.2 via the probability integral transform.

Let Φ(x)\Phi(x) denote the CDF of a N(0,1)N(0,1) random variable, and let Fμ,σ2[a,b]F_{\mu,\sigma^{2}}^{[a,b]} denote the CDF of TN(μ,σ,a,b)TN(\mu,\sigma,a,b), i.e.:

Then FηTμ, ηTΣη[V−,V+](ηTy)F_{\eta^{T}\mu,\ \eta^{T}\Sigma\eta}^{[{\cal V}^{-},{\cal V}^{+}]}(\eta^{T}y) is a pivotal quantity, conditional on {Ay≤b}\{Ay\leq b\}:

where V−{\cal V}^{-} and V+{\cal V}^{+} are defined in (10) and (11).

Inference for marginal screening

In this section, we apply the theory summarized in Sections 4 and 5 to marginal screening. In particular, we will construct confidence intervals for the selected variables.

To summarize the developments so far, recall that our model (1) says that y∼N(μ,σ2I)y\sim N(\mu,\sigma^{2}I). The distribution of interest is y∣{E^(y)=E}y|\{\hat{E}(y)=E\}, and by Theorem 4.1, this is equivalent to y∣{A(S,s)z≤b(S,s)}y|{\{A(S,s)z\leq b(S,s)\}}, where y∼N(μ,σ2I)y\sim N(\mu,\sigma^{2}I). By applying Theorem 5.3, we obtain the pivotal quantity

for any η\eta, where V−{\cal V}^{-} and V+{\cal V}^{+} are defined in (10) and (11).

In this section, we describe how to form confidence intervals for the components of βS^⋆=(XS^TXS^)−1XS^Tμ\beta^{\star}_{\hat{S}}=(X_{\hat{S}}^{T}X_{\hat{S}})^{-1}X_{\hat{S}}^{T}\mu. The best linear predictor of μ\mu that uses only the selected variables is βS^⋆\beta^{\star}_{\hat{S}} , and β^S^=(XS^TXS^)−1XS^Ty\hat{\beta}_{\hat{S}}=(X_{\hat{S}}^{T}X_{\hat{S}})^{-1}X_{\hat{S}}^{T}y is an unbiased estimate of βS^⋆\beta^{\star}_{\hat{S}}. In this section, we propose hypothesis tests and confidence intervals for βS^⋆\beta^{\star}_{\hat{S}}. If we choose

then ηjTμ=βj∈S^⋆\eta_{j}^{T}\mu=\beta_{j\in\hat{S}}^{\star}, so the above framework provides a method for inference about the jthj^{\text{th}} variable in the model S^\hat{S}. This choice of η\eta is not fixed before marginal screening selects S^\hat{S}, but it is measurable with respect to the σ\sigma-algebra generated by the partition. Since it is measurable, η\eta is constant on each partition, so the pivot is uniformly distributed on each element of the partition, and thus uniformly distributed for all yy.

If we assume the linear model μ=Xβ0\mu=X\beta^{0} for some β0∈Rp\beta^{0}\in\mathbf{R}^{p}, S0:=support(β0)⊂S^S^{0}:=\text{support}(\beta^{0})\subset\hat{S}, and XS^X_{\hat{S}} is full rank, then by the following computation βS^⋆=βS^0\beta^{\star}_{\hat{S}}=\beta_{\hat{S}}^{0}:

In , the screening property S0⊂S^S^{0}\subset\hat{S} for the marginal screening algorithm is established under mild conditions. Thus under the screening property, our method provides hypothesis tests and confidence intervals for βS^0\beta_{\hat{S}}^{0}.

By applying Theorem 5.3, we obtain the following (conditional) pivot for βj∈S^⋆\beta^{\star}_{j\in\hat{S}}:

The quantities jj and ηj\eta_{j} are both random through E^\hat{E}, a quantity which is fixed after conditioning, therefore Theorem 5.3 holds even for this choice of η\eta.

Consider testing the hypothesis H0:βj∈S^⋆=βjH_{0}:\beta^{\star}_{j\in\hat{S}}=\beta_{j}. A valid test statistic is given by Fβj, σ2∣∣ηj∣∣2[V−,V+](ηjTy)F_{\beta_{j},\ \sigma^{2}||\eta_{j}||^{2}}^{[{\cal V}^{-},{\cal V}^{+}]}(\eta_{j}^{T}y), which is uniformly distributed under the null hypothesis and y∣{E^(y)=E}y|\{\hat{E}(y)=E\}. Thus, this test would reject when Fβj, σ2∣∣ηj∣∣2[V−,V+](ηjTy)>1−α2F_{\beta_{j},\ \sigma^{2}||\eta_{j}||^{2}}^{[{\cal V}^{-},{\cal V}^{+}]}(\eta_{j}^{T}y)>1-\frac{\alpha}{2} or Fβj, σ2∣∣ηj∣∣2[V−,V+](ηjTy)<α2F_{\beta_{j},\ \sigma^{2}||\eta_{j}||^{2}}^{[{\cal V}^{-},{\cal V}^{+}]}(\eta_{j}^{T}y)<\frac{\alpha}{2}.

The test of H0:βj∈S^⋆=βjH_{0}:\beta^{\star}_{j\in\hat{S}}=\beta_{j} that accepts when

Under H0H_{0}, we have βj∈S^⋆=βj\beta^{\star}_{j\in\hat{S}}=\beta_{j}, so by (17) F_{\beta_{j},\ \sigma^{2}||\eta_{j}||^{2}}^{[{\cal V}^{-},{\cal V}^{+}]}(\eta_{j}^{T}y)\big{|}\{\hat{E}(y)=E\} is uniformly distributed. Thus

and the type 1 error is exactly α\alpha. Under H0H_{0}, but not conditional on selection event E^\hat{E}, we have

For each element of the partition EE, the conditional (on selection) hypothesis test is level 1−α1-\alpha, so by summing over the partition the unconditional test is level 1−α1-\alpha. ∎

Our hypothesis test is not conservative, in the sense that the type 1 error is exactly α\alpha; also, it is non-asymptotic, since the statement holds for fixed nn and pp. We summarize the hypothesis test in this section in the following algorithm.

2 Confidence intervals for selected variables

Next, we discuss how to obtain confidence intervals for βj∈S^⋆\beta^{\star}_{j\in\hat{S}}. The standard way to obtain an interval is to invert a pivotal quantity . In other words, since

one can define a (1−α)(1-\alpha) (conditional) confidence interval for βj,E^⋆\beta_{j,\hat{E}}^{\star} as

In fact, FF is monotone decreasing in xx, so to find its endpoints, one need only solve for the root of a smooth one-dimensional function. The monotonicity is a consequence of the fact that the truncated Gaussian distribution is a natural exponential family and hence has monotone likelihood ratio in μ\mu .

We now formalize the above observations in the following result, an immediate consequence of Theorem 5.3.

Let ηj\eta_{j} be defined as in (16), and let Lα=Lα(ηj,(S^,s^))L_{\alpha}=L_{\alpha}(\eta_{j},(\hat{S},\hat{s})) and Uα=Uα(ηj,(S^,s^))U_{\alpha}=U_{\alpha}(\eta_{j},(\hat{S},\hat{s})) be the (unique) values satisfying

Then [Lα,Uα][L_{\alpha},U_{\alpha}] is a (1−α)(1-\alpha) confidence interval for βj∈S^⋆\beta^{\star}_{j\in\hat{S}}, conditional on E^\hat{E}:

The confidence region of βj∈S^⋆\beta^{\star}_{j\in\hat{S}} is the set of βj\beta_{j} such that the test of H0:βj∈S^⋆H_{0}:\beta^{\star}_{j\in\hat{S}} accepts at the 1−α1-\alpha level. The function Fx, σ2∣∣ηj∣∣2[V−,V+](ηjTy)F_{x,\ \sigma^{2}||\eta_{j}||^{2}}^{[{\cal V}^{-},{\cal V}^{+}]}(\eta_{j}^{T}y) is monotone in xx, so solving for LαL_{\alpha} and UαU_{\alpha} identify the most extreme values where H0H_{0} is still accepted. This gives a 1−α1-\alpha confidence interval. ∎

In relation to the literature on False Coverage Rate (FCR) , our procedure also controls the FCR.

Furthermore, the FCR of the intervals {[Lαj,Uαj]}j∈E^\left\{[L_{\alpha}^{j},U_{\alpha}^{j}]\right\}_{j\in\hat{E}} is α\alpha.

By (20), the conditional coverage of the confidence intervals are 1−α1-\alpha. The coverage holds for every element of the partition {E^(y)=E}\{\hat{E}(y)=E\}, so

We summarize the algorithm for selecting and constructing confidence intervals below.

3 Experiments on Diabetes dataset

In Figure 1, we have already seen that the confidence intervals constructed using Algorithm 3 have exactly 1−α1-\alpha coverage proportion. In this section, we perform an experiment on real data where the linear model does not hold, the noise is not Gaussian, and the noise variance is unknown. The diabetes dataset contains n=442n=442 diabetes patients measured on p=10p=10 baseline variables . The baseline variables are age, sex, body mass index, average blood pressure, and six blood serum measurements, and the response yy is a quantitative measure of disease progression measured one year after the baseline. The goal is to use the baseline variables to predict yy, the measure of disease progression after one year, and determine which baseline variables are statistically significant for predicting yy.

Extensions

The purpose of this section is to illustrate the broad applicability of the condition on selection framework. This framework was first proposed in to form valid hypothesis tests and confidence intervals after model selection via the Lasso. However, the framework is not restricted to the Lasso, and we have shown how to apply it to marginal screening. For expository purposes, we focused the paper on marginal screening where the framework is particularly easy to understand. In the rest of this section, we show how to apply the framework to marginal screening+Lasso, orthogonal matching pursuit, and non-negative least squares. This is a non-exhaustive list of selection procedures where the condition on selection framework is applicable, but we hope this incomplete list emphasizes the ease of constructing tests and confidence intervals post-model selection via conditioning.

The marginal screening+Lasso procedure was introduced in as a variable selection method for the ultra-high dimensional setting of p=O(enk)p=O(e^{n^{k}}). Fan et al. recommend applying the marginal screening algorithm with k=n−1k=n-1, followed by the Lasso on the selected variables. This is a two-stage procedure, so to properly account for the selection we must encode the selection event of marginal screening followed by Lasso. This can be done by representing the two stage selection as a single event. Let (S^m,s^m)(\hat{S}_{m},\hat{s}_{m}) be the variables and signs selected by marginal screening, and the (S^L,z^L)(\hat{S}_{L},\hat{z}_{L}) be the variables and signs selected by Lasso . In Proposition 2.2 of , it is shown how to encode the Lasso selection event (S^L,z^L)(\hat{S}_{L},\hat{z}_{L}) as a set of constraints {ALy≤bL}\{A_{L}y\leq b_{L}\} The Lasso selection event is with respect to the Lasso optimization problem after marginal screening., and in Section 4 we showed how to encode the marginal screening selection event (S^m,s^m)(\hat{S}_{m},\hat{s}_{m}) as a set of constraints {Amy≤bm}\{A_{m}y\leq b_{m}\}. Thus the selection event of marginal screening+Lasso can be encoded as {ALy≤bL,Amy≤bm}\{A_{L}y\leq b_{L},A_{m}y\leq b_{m}\}. Using these constraints, the hypothesis test and confidence intervals described in Algorithms 2 and 3 are valid for marginal screening+Lasso.

2 Orthogonal Matching Pursuit

Orthogonal matching pursuit (OMP) is a commonly used variable selection method. At each iteration, OMP selects the variable most correlated with the residual rr, and then recomputes the residual using the residual of least squares using the selected variables. The description of the OMP algorithm is given in Algorithm 4.

Similar to Section 4, we can represent the OMP selection event as a set of linear constraints on yy.

The selection event encodes that OMP selected a certain variable and the sign of the correlation of that variable with the residual, at steps 11 to kk. The primary difference between the OMP selection event and the marginal screening selection event is that the OMP event also describes the order at which the variables were chosen. The marginal screening event only describes that the variable was among the top kk most correlated, and not whether a variable was the most correlated or kthkth most correlated.

Since the selection event can be represented as constraints on yy, the hypothesis test and confidence intervals described in Algorithms 2 and 3 are valid for OMP selected β^S^\hat{\beta}_{\hat{S}}.

3 Nonnegative Least Squares

Non-negative least squares (NNLS) is a simple modification of the linear regression estimator with non-negative constraints on β\beta:

Under a positive eigenvalue conditions on XX, several authors have shown that NNLS is comprable to the Lasso in terms of prediction and estimation errors. The NNLS estimator also does not have any tuning parameters, since the sign constraint provides a natural form of regularization. NNLS has found applications when modeling non-negative data such as prices, incomes, count data. Non-negativity constraints arise naturally in non-negative matrix factorization, signal deconvolution, spectral analysis, and network tomography; we refer to for a comprehensive survey of the applications of NNLS.

We show how our framework can be used to form exact hypothesis tests and confidence intervals for NNLS estimated coefficients. The primal dual solution pair (β^,λ^)(\hat{\beta},\hat{\lambda}) is a solution iff the KKT conditions are satisfied,

Let S^={i:−xiT(y−Xβ^)=0}\hat{S}=\{i:-x_{i}^{T}(y-X\hat{\beta})=0\}. By complementary slackness β^−S^=0\hat{\beta}_{-\hat{S}}=0, where −S^-\hat{S} is the complement to the “active” variables S^\hat{S} chosen by NNLS. Given the active set we can solve the KKT equation for the value of β^S^\hat{\beta}_{\hat{S}},

which is a linear contrast of yy. The NNLS selection event is

The selection event encodes that for a given yy the NNLS optimization program will select a subset of variables S^(y)\hat{S}(y). Similar to the case in OMP and marginal screening, we can use Algorithms 2 and 3, since the selection event is represented by a set of linear constraints {y:A(S^)y≤0}\{y:A(\hat{S})y\leq 0\}.

Conclusion

Due to the increasing size of datasets, marginal screening has become an important method for fast variable selection. However, the standard hypothesis tests and confidence intervals used in linear regression are invalid after using marginal screening to select important variables. We have described a method to perform hypothesis and form confidence intervals after marginal screening. The conditional on selection framework is not restricted to marginal screening, and also applies to OMP, marginal screening + Lasso, and NNLS.

Acknowledgements

Jonathan Taylor was supported in part by NSF grant DMS 1208857 and AFOSR grant 113039. Jason Lee was supported by a NSF graduate fellowship, and a Stanford Graduate Fellowship.

References