High-dimensional variable selection
The Annals of Statistics
2009, Vol. 37, No. 5A, 2178–2201
DOI: 10.1214/08-AOS646
c© Institute of Mathematical Statistics, 2009
HIGH-DIMENSIONAL VARIABLE SELECTION
By Larry Wasserman and Kathryn Roeder1
Carnegie Mellon University
This paper explores the following question: what kind of statis-
tical guarantees can be given when doing variable selection in high-
dimensional models? In particular, we look at the error rates and
power of some multi-stage regression methods. In the first stage we
fit a set of candidate models. In the second stage we select one model
by cross-validation. In the third stage we use hypothesis testing to
eliminate some variables. We refer to the first two stages as “screen-
ing” and the last stage as “cleaning.” We consider three screening
methods: the lasso, marginal regression, and forward stepwise regres-
sion. Our method gives consistent variable selection under certain
conditions.
1. Introduction. Several methods have been developed lately for high-
dimensional linear regression such as the lasso [Tibshirani (1996)], Lars
[Efron et al. (2004)] and boosting [Bühlmann (2006)]. There are at least
two different goals when using these methods. The first is to find models
with good prediction error. The second is to estimate the true “sparsity
pattern,” that is, the set of covariates with nonzero regression coefficients.
These goals are quite different and this paper will deal with the second goal.
(Some discussion of prediction is in the Appendix.) Other papers on this
topic include Meinshausen and Bühlmann (2006), Candes and Tao (2007),
Wainwright (2006), Zhao and Yu (2006), Zou (2006), Fan and Lv (2008),
Meinshausen and Yu (2008), Tropp (2004, 2006), Donoho (2006) and Zhang
and Huang (2006). In particular, the current paper builds on ideas in Mein-
shausen and Yu (2008) and Meinshausen (2007).
Let (X1, Y1), . . . , (Xn, Yn) be i.i.d. observations from the regression model
Yi =X
i β + εi,(1)
Received June 2008; revised August 2008.
1Supported by NIH Grant MH057881.
AMS 2000 subject classifications. Primary 62J05; secondary 62J07.
Key words and phrases. Lasso, stepwise regression, sparsity.
This is an electronic reprint of the original article published by the
Institute of Mathematical Statistics in The Annals of Statistics,
2009, Vol. 37, No. 5A, 2178–2201. This reprint differs from the original in
pagination and typographic detail.
http://arxiv.org/abs/0704.1139v2
http://www.imstat.org/aos/
http://dx.doi.org/10.1214/08-AOS646
http://www.imstat.org
http://www.ams.org/msc/
http://www.imstat.org
http://www.imstat.org/aos/
http://dx.doi.org/10.1214/08-AOS646
2 L. WASSERMAN AND K. ROEDER
where ε ∼ N(0, σ2), Xi = (Xi1, . . . ,Xip)T ∈ Rp and p = pn > n. Let X be
the n × p design matrix with jth column X•j = (X1j , . . . ,Xnj)T and let
Y = (Y1, . . . , Yn)
T . Let
D= {j :βj 6= 0}
be the set of covariates with nonzero regression coefficients. Without loss
of generality, assume that D = {1, . . . , s} for some s. A variable selection
procedure D̂n maps the data into subsets of S = {1, . . . , p}.
The main goal of this paper is to derive a procedure D̂n such that
lim sup
P(D̂n ⊂D)≥ 1− α,(2)
that is, the asymptotic type I error is no more than α. Note that throughout
the paper we use ⊂ to denote nonstrict set-inclusion. Moreover, we want
D̂n to have nontrivial power. Meinshausen and Bühlmann (2006) control a
different error measure. Their method guarantees lim supn→∞P(D̂n ∩ V 6=
∅)≤ α where V is the set of variables not connected to Y by any path in
an undirected graph.
Our procedure involves three stages. In stage I we fit a suite of candidate
models, each model depending on a tuning parameter λ,
S = {Ŝn(λ) :λ ∈ Λ}.
In stage II we select one of those models Ŝn using cross-validation to select λ̂.
In stage III we eliminate some variables by hypothesis testing. Schematically,
stage I−→ S stage II−→︸ ︷︷ ︸
screen
stage III−→ D̂n︸ ︷︷ ︸
clean
Genetic epidemiology provides a natural setting for applying screen and
clean. Typically, the number of subjects, n, is in the thousands, while p
ranges from tens of thousands to hundereds of thousands of genetic fea-
tures. The number of genes exhibiting a detectable association with a trait
is extremely small. Indeed, for type I diabetes only ten genes have exhibited
a reproducible signal [Wellcome Trust (2007)]. Hence, it is natural to assume
that the true model is sparse. A common experimental design involves a 2-
stage sampling of data, with stages 1 and 2 corresponding to the screening
and cleaning processes, respectively.
In stage 1 of a genetic association study, n1 subjects are sampled and one
or more traits such as bone mineral density are recorded. Each subject is
also measured at p locations on the chromosomes. These genetic covariates
usually have two forms in the population due to variability at a single nu-
cleotide and hence are called single nucleotide polymorphisms (SNPs). The
distinct forms are called alleles. Each covariate takes on a value (0, 1 or 2)
VARIABLE SELECTION 3
indicating the number of copies of the less common allele observed. For a
well-designed genetic study, individual SNPs are nearly uncorrelated unless
they are physically located in very close proximity. This feature makes it
much easier to draw causal inferences about the relationship between SNPs
and quantitative traits. It is standard in the field to infer that an association
discovered between a SNP and a quantitative trait implies a causal genetic
variant is physically located near the one exhibiting association. In stage 2,
n2 subjects are sampled at a subset of the SNPs assessed in stage 1. SNPs
measured in stage 2 are often those that achieved a test statistic that ex-
ceeded a predetermined threshold of significance in stage 1. In essence, the
two stage design pairs naturally with a screen and clean procedure.
For the screen and clean procedure, it is essential that Ŝn has two prop-
erties as n→∞ as follows:
P(D ⊂ Ŝn)→ 1(3)
|Ŝn|= oP (n),(4)
where |M | denotes the number of elements in a setM . Condition (3) ensures
the validity of the test in stage III while condition (4) ensures that the power
of the test is not too small. Without condition (3), the hypothesis test in
stage III would be biased. We will see that the power goes to 1, so taking
α = αn→ 0 implies consistency: P(D̂n =D)→ 1. For fixed α, the method
also produces a confidence sandwich for D, namely,
lim inf
P(D̂n ⊂D ⊂ Ŝn)≥ 1−α.
To fit the suite of candidate models, we consider three methods. In method
Ŝn(λ) = {j : β̃j(λ) 6= 0},
where β̃j(λ) is the lasso estimator, the value of β that minimizes
(Yi −XTi β)2 + λ
|βj |.
In method 2, take Ŝn(λ) to be the set of variables chosen by forward stepwise
regression after λ steps. In method 3, marginal regression, we take
Ŝn = {j : |µ̂j|> λ},
where µ̂j is the marginal regression coefficient from regressing Y onXj . (This
is equivalent to ordering by the absolute t-statistics since we will assume that
the covariates are standardized.) These three methods are very similar to
basis pursuit, orthogonal matching pursuit and thresholding [see, e.g., Tropp
(2004, 2006) and Donoho (2006)].
4 L. WASSERMAN AND K. ROEDER
Notation. Let ψ =minj∈D |βj |. Define the loss of any estimator β̂ by
L(β̂) =
(β̂ − β)TXTX(β̂ − β) = (β̂ − β)T Σ̂n(β̂ − β),(5)
where Σ̂n = n
−1XTX . For convenience, when β̂ ≡ β̂(λ) depends on λ we
write L(λ) instead of L(β̂(λ)). If M ⊂ S, let XM be the design matrix
with columns (X•j : j ∈M) and let β̂M = (XTMXM )−1XTMY denote the least-
squares estimator, assuming it is well defined. Note that our use ofX•j differs
from standard ANOVA notation. Write Xλ instead of XM whenM = Ŝn(λ).
When convenient, we extend β̂M to length p by setting β̂M (j) = 0, for j /∈M .
We use the norms
v2j , ‖v‖1 =
|vj | and ‖v‖∞ =max
|vj |.
If C is any square matrix, let φ(C) and Φ(C) denote the smallest and
largest eigenvalues of C. Also, if k is an integer, define
φn(k) = min
M : |M |=k
XTMXM
and Φn(k) = max
M : |M |=k
XTMXM
We will write zu for the upper quantile of a standard normal, so that P(Z >
zu) = u where Z ∼N(0,1).
Our method will involve splitting the data randomly into three groups
D1, D2 and D3. For ease of notation, assume the total sample size is 3n and
that the sample size of each group is n.
Summary of assumptions. We will use the following assumptions through-
out except in Section 8:
(A1) Yi =X
i β + εi where εi ∼N(0, σ2), for i= 1, . . . , n.
(A2) The dimension pn of X satisfies pn→∞ and pn ≤ c1en
, for some
c1 > 0 and 0≤ c2 < 1.
(A3) s≡ |{j :βj 6= 0}|=O(1) and ψ =min{|βj | :βj 6= 0}> 0.
(A4) There exist positive constants C0,C1 and κ such that
P(lim supn→∞Φn(n)≤C0) = 1 and P(lim infn→∞ φn(C1 logn)≥ κ) = 1. Also,
P(φn(n)> 0) = 1, for all n.
(A5) The covariates are standardized: E(Xij) = 0 and E(X
ij) = 1. Also,
there exists 0<B <∞ such that P(|Xjk| ≤B) = 1.
For simplicity, we include no intercepts in the regressions. The assump-
tions can be weakened at the expense of more complicated proofs. In par-
ticular, we can let s increase with n and ψ decrease with n. Similarly, the
normality and constant variance assumptions can be relaxed.
VARIABLE SELECTION 5
2. Error control. Define the type I error rate q(D̂n) = P(D̂n ∩ Dc 6=
∅) and the asymptotic error rate lim supn→∞ q(D̂n). We define the power
π(D̂n) = P(D ⊂ D̂n) and the average power
πav =
P(j ∈ D̂n).
It is well known that controlling the error rate is difficult for at least
three reasons: correlation of covariates, high-dimensionality of the covariate
and unfaithfulness (cancellations of correlations due to confounding). Let us
briefly review these issues.
It is easy to construct examples where q(D̂n)≤ α implies that π(D̂n)≈ α.
Consider the following two models for random variables Z = (Y,X1,X2):
Model 1
X1 ∼N(0,1),
Y = ψX1 +N(0,1),
X2 = ρX1 +N(0, τ
Model 2
X2 ∼N(0,1),
Y = ψX2 +N(0,1),
X1 = ρX2 +N(0, τ
Under models 1 and 2, the marginal distribution of Z is P1 =N(0,Σ1)
and P2 =N(0,Σ2), where
ψ2 + 1 ψ ρψ
ψ 1 ρ
ρψ ρ ρ2 + τ2
 , Σ2 =
ψ2 + 1 ρψ ψ
ρψ ρ2 + τ2 ρ
ψ ρ 1
Given any ε > 0, we can choose ρ sufficiently close to 1 and τ sufficiently close
to 0 such that Σ1 and Σ2 are as close as we like, and hence, d(P
1 , P
2 )< ε
where d is total variation distance. It follows that
P2(2 /∈ D̂)≥ P1(2 /∈ D̂)− ε≥ 1− α− ε.
Thus, if q ≤ α, then the power is less than α+ ε.
Dimensionality is less of an issue thanks to recent methods. Most methods,
including those in this paper, allow pn to grow exponentially. But all the
methods require some restrictions on the number s of nonzero βj ’s. In other
words, some sparsity assumption is required. In this paper we take s fixed
and allow pn to grow.
False negatives can occur during screening due to cancellations of cor-
relations. For example, the correlation between Y and X1 can be 0 even
when β1 is huge. This problem is called unfaithfulness in the causality liter-
ature [see Spirtes, Glymour and Scheines (2001) and Robins et al. (2003)].
False negatives during screening can lead to false positives during the second
stage.
6 L. WASSERMAN AND K. ROEDER
Let µ̂j denote the regression coefficient from regressing Y on Xj . Fix j ≤ s
and note that
µj ≡ E(µ̂j) = βj +
k 6=j
1≤k≤s
βkρkj,
where ρkj = corr(Xk,Xj). If
k 6=j
1≤k≤s
βkρkj ≈−βj,
then µj ≈ 0 no matter how large βj is. This problem can occur even when
n is large and p is small.
For example, suppose that β = (10,−10,0,0) and that ρ(Xi,Xj) = 0 ex-
cept that ρ(X1,X2) = ρ(X1,X3) = ρ(X2,X4) = 1− ε, where ε > 0 is small.
Then,
β = (10,−10,0,0), but µ≈ (0,0,10,−10).
Marginal regression is extremely susceptible to unfaithfulness. The lasso
and forward stepwise, less so. However, unobserved covariates can induce
unfaithfulness in all the methods.
3. Loss and cross-validation. Let Xλ = (X•j : j ∈ Ŝn(λ)) denote the de-
sign matrix corresponding to the covariates in Ŝn(λ) and let β̂(λ) be the
least-squares estimator for the regression restricted to Ŝn(λ), assuming the
estimator is well defined. Hence, β̂(λ) = (XTλXλ)
−1XTλ Y . More generally,
β̂M is the least-squares estimator for any subset of variables M . When con-
venient, we extend β̂(λ) to length p by setting β̂j(λ) = 0, for j /∈ Ŝn(λ).
3.1. Loss. Now we record some properties of the loss function. The first
part of the following lemma is essentially Lemma 3 of Meinshausen and Yu
(2008).
Lemma 3.1. Let M+m = {M ⊂ S : |M | ≤m,D⊂M}. Then,
L(β̂M )≤
4m log p
nφn(m)
→ 1.(6)
Let M−m = {M ⊂ S : |M | ≤m,D 6⊂M}. Then,
L(β̂M )≥ ψ2φn(m+ s)
→ 1.(7)
VARIABLE SELECTION 7
3.2. Cross-validation. Recall that the data have been split into groups
D1, D2 and D3 each of size n. Construct β̂(λ) from D1 and let
L̂(λ) =
Xi∈D2
(Yi −XTi β̂(λ))
We would like L̂(λ) to order the models the same way as the true loss
L(λ) [defined after (5)]. This requires that, asymptotically, L̂(λ)−L(λ)≈ δn,
where δn does not involve λ. The following bounds will be useful. Note that
L(λ) and L̂(λ) are both step functions that only change value when a variable
enters or leaves the model.
Theorem 3.2. Suppose that maxλ∈Λn |Ŝn(λ)| ≤ kn. Then, there exists
a sequence of random variables δn =OP (1) that do not depend on λ or X,
such that, with probability tending to 1,
|L(λ)− L̂(λ)− δn|=OP
n1−c2
4. Multi-stage methods. Themulti-stage methods use the following steps.
As mentioned earlier, we randomly split the data into three parts, D1, D2
and D3, which we take to be of equal size:
1. Stage I. Use D1 to find Ŝn(λ), for each λ.
2. Stage II. Use D2 to find λ̂ by cross-validation, and let Ŝn = Ŝn(λ̂).
3. Stage III. Use D3 to find the least-squares estimate β̂ for the model Ŝn.
D̂n = {j ∈ Ŝn : |Tj |> cn},
where Tj is the usual t-statistic, cn = zα/2m and m= |Ŝn|.
4.1. The lasso. The lasso estimator [Tibshirani (1996)] β̃(λ) minimizes
Mλ(λ) =
(Yi −XTi β)2 + λ
|βj |
and let Ŝn(λ) = {j : β̃j(λ) 6= 0}. Recall that β̂(λ) is the least-squares estima-
tor using the covariates in Ŝn(λ).
Let kn =A logn where A> 0 is a positive constant.
Theorem 4.1. Assume that (A1)–(A5) hold. Let Λn = {λ : |Ŝn(λ)| ≤
kn}. Then:
1. The true loss overfits: P(D ⊂ Ŝn(λ∗))→ 1 where λ∗ = argminλ∈Λn L(λ).
8 L. WASSERMAN AND K. ROEDER
2. Cross-validation also overfits: P(D⊂ Ŝn(λ̂))→ 1 where λ̂= argminλ∈Λn L̂(λ).
3. Type I error is controlled: lim supn→∞ P(D
c ∩ D̂n 6=∅) ≤ α.
If we let α= αn→ 0, then D̂n is consistent for variable selection.
Theorem 4.2. Assume that (A1)–(A5) hold. Let αn→ 0 and
∞. Then, the multi-stage lasso is consistent,
P(D̂n =D)→ 1.(10)
The next result follows directly. The proof is thus omitted.
Theorem 4.3. Assume that (A1)–(A5) hold. Let α be fixed. Then,
(D̂n, Ŝn) forms a confidence sandwich
lim inf
P(D̂n ⊂D ⊂ Ŝn)≥ 1−α.(11)
Remark 4.4. This confidence sandwich is expected to be conservative
in the sense that the coverage can be much larger than 1− α.
4.2. Stepwise regression. Let kn = A logn for some A > 0. The version
of stepwise regression we consider is as follows:
1. Initialize: Res = Y , λ= 0, Ŷ = 0 and Ŝn(λ) =∅.
2. Let λ← λ+1. Compute µ̂j = n−1〈Xj ,Res〉 for j = 1, . . . , p.
3. Let J = argmaxj |µ̂j |. Set Ŝn(λ) = {Ŝn(λ−1), J}. Set Ŷ =Xλβ̂(λ) where
β̂λ = (X
−1XTλ Y , and let Res = Y − Ŷ .
4. If λ= kn, stop. Otherwise, go to step 2.
For technical reasons, we assume that the final estimator xT β̂ is truncated
to be no larger than B. Note that λ is discrete and Λn = {0,1, . . . , kn}.
Theorem 4.5. With Ŝn(λ) defined as above, the statements of Theorems
4.1, 4.2 and 4.3 hold.
4.3. Marginal regression. This is probably the oldest, simplest and most
common method. It is quite popular in gene expression analysis. It used to
be regarded with some derision but has enjoyed a revival. A version appears
in a recent paper by Fan and Lv (2008). Let Ŝn(λ) = {j : |µ̂j | ≥ λ} where
µ̂j = n
−1〈Y,X•j〉.
Let µj = E(µ̂j), and let µ(j) denote the value of µ ordered by their absolute
values,
|µ(1)| ≥ |µ(2)| ≥ · · · .
VARIABLE SELECTION 9
Theorem 4.6. Let kn →∞ with kn = o(
n). Let Λn = {λ : |Ŝn(λ)| ≤
kn}. Assume that
|µj|> |µ(kn)|.(12)
Then, the statements of Theorems 4.1, 4.2 and 4.3 hold.
Assumption (12) limits the degree of unfaithfulness (small partial corre-
lations induced by cancellation of parameters). Large values of kn weaken
assumption (12), thus making the method more robust to unfaithfulness, but
at the expense of lower power. Fan and Lv (2008) make similar assumptions.
They assume that there is a C > 0 such that |µj| ≥ C|βj | for all j, which
rules out unfaithfulness. However, they do not explicitly relate the values of
µj for j ∈D to the values outside D as we have done. On the other hand,
they assume that Z =Σ−1/2X has a spherically symmetric distribution. Un-
der this assumption and their faithfulness assumption, they deduce that the
µj ’s outside D cannot strongly dominate the µj ’s within D. We prefer to
simply make this an explicit assumption without placing distributional as-
sumptions on X . At any rate, any method that uses marginal regressions as
a starting point must make some sort of faithfulness assumptions to succeed.
4.4. Modifications. Let us now discuss a few modifications of the basic
method. First, consider splitting the data only into two groups, D1 and D2.
Then do the following steps:
1. Stage I. Find Ŝn(λ) for λ ∈Λn, where |Ŝn(λ)| ≤ kn for each λ ∈ Λn using
2. Stage II. Find λ̂ by cross-validation, and let Ŝn = Ŝn(λ̂) using D2.
3. Stage III. Find the least-squares estimate β̂
using D2. Let D̂n = {j ∈
Ŝn : |Tj |> cn}, where Tj is the usual t-statistic.
Theorem 4.7. Choosing
log logn
2kn log(2pn)
controls asymptotic type I error.
The critical value in (13) is hopelessly large and it does not appear it can
be substantially reduced. We present this mainly to show the value of the
extra data-splitting step. It is tempting to use the same critical value as in
the tri-split case, namely cn = zα/2m where m = |Ŝn|, but we suspect this
will not work in general. However, it may work under extra conditions.
10 L. WASSERMAN AND K. ROEDER
5. Application. As an example, we illustrate an analysis based on part
of the osteoporotic fractures in men study [MrOS, Orwoll et al. (2005)]. A
sample of 860 men were measured at a large number of genes and outcome
measures. We consider only 296 SNPs which span 30 candidate genes for
bone mineral density. An aim of the study was to identify genes associated
with bone mineral density that could help in understanding the genetic basis
of osteoporosis in men. Initial analyses of this subset of the data revealed
no SNPs with a clear pattern of association with the phenotype; however,
three SNPs, numbered (67), (277) and (289), exhibited some association in
the screening of the data. To further explore the effacacy of the lasso screen
and clean procedure, we modified the phenotype to enhance this weak signal
and then reanalyzed the data to see if we could detect this planted signal.
We were interested in testing for main effects and pairwise interactions
in these data; however, including all interactions results in a model with
43,660 additional terms, which is not practical for this sample size. As a
compromise, we selected 2 SNPs per gene to model potential interaction
effects. This resulted in a model with a total of 2066 potential coefficients,
including 296 main effects and 1770 interaction terms. With this model, our
initial screen detected 10 terms, including the 3 enhanced signals, 2 other
main effects and 5 interactions. After cleaning, the final model detected the
3 enhanced signals and no other terms.
6. Simulations. To further explore the screen and clean procedures, we
conducted simulation experiments with four models. For each model Yi =
XTi β + εi where the measurement errors, εi and ε
ij , are i.i.d. normal(0,1)
and the covariates Xij ’s are normal(0,1) (except for model D). Models differ
in how Yi is linked to Xi and the dependence structure of the Xi’s. Models
A, B and C explore scenarios with moderate and large p, while model D
focuses on confounding and unfaithfullness, as follows:
(A) Null model: β = (0, . . . ,0) and the Xij ’s are i.i.d.
(B) Triangle model: βj = δ(10− j), j = 1, . . . ,10, βj = 0, j > 10 and Xij ’s are
i.i.d.
(C) Correlated Triangle model: as B, but withXij(+1) = ρXij+(1−ρ2)1/2ε∗ij ,
for j > 1, and ρ= 0.5.
(D) Unfaithful model: Yi = β1Xi1 + β2Xi2 + εi, for β1 = −β2 = 10, where
the Xij ’s are i.i.d. for j = {1,5,6,7,8,9,10}, but Xi2 = ρXi1 + τε∗i2,
Xi3 = ρXi1 + τε
i10, and Xi4 = ρXi2 + τε
i11, for τ = 0.01 and ρ= 0.95.
We used a maximum model size of kn = n
1/2 which technically goes be-
yond the theory but works well in practice. Prior to analysis, the covariates
are scaled so that each has mean 0 and variance 1. The tests were initially
performed using a third of the data for each of the 3 stages of the proce-
dure (Table 1, top half, 3 splits). For models A, B and C, each approach
VARIABLE SELECTION 11
Table 1
Size and power of screen and clean procedures using lasso, stepwise and marginal
regression for the screening step. For all procedures α= 0.05. For p= 100, δ = 0.5 and
for p= 1000, δ = 1.5. Reported power is πav. The top 8 rows of simulations were
conducted using three stages as described in Section 4, with a third of the data used for
each stage. The bottom 8 rows of simulations were conducted splitting the data in half,
using the first portion with leave-one-out cross validation for stages 1 and 2 and the
second portion for cleaning
Size Power
Splits n p Model Lasso Step Marg Lasso Step Marg
2 100 100 A 0.005 0.001 0.004 0 0 0
2 100 100 B 0.01 0.02 0.03 0.62 0.62 0.31
2 100 100 C 0.001 0.01 0.01 0.77 0.57 0.21
2 100 10 D 0.291 0.283 0.143 0.08 0.08 0.04
2 100 1000 A 0.001 0.002 0.010 0 0 0
2 100 1000 B 0.002 0.020 0.010 0.17 0.09 0.11
2 100 1000 C 0.02 0.14 0.01 0.27 0.15 0.11
2 1000 10 D 0.291 0.283 0.143 0.08 0.08 0.04
3 100 100 A 0.040 0.050 0.030 0 0 0
3 100 100 B 0.02 0.01 0.02 0.91 0.90 0.56
3 100 100 C 0.03 0.04 0.03 0.91 0.88 0.41
3 100 10 D 0.382 0.343 0.183 0.16 0.18 0.09
3 100 1000 A 0.035 0.045 0.040 0 0 0
3 100 1000 B 0.045 0.020 0.035 0.57 0.66 0.29
3 100 1000 C 0.06 0.070 0.020 0.74 0.65 0.19
3 1000 10 D 0.481 0.486 0.187 0.17 0.17 0.13
has type I error less than α, except the stepwise procedure which has trou-
ble with model C when n = p = 100. We also calculated the false positive
rate and found it to be very low (about 10−4 when p= 100 and 10−5 when
p= 1000) indicating that even when a type I error occurs, only a very small
number of terms are included erroneously. The lasso screening procedure
exhibited a slight power advantage over the stepwise procedure. Both meth-
ods dominated the marginal approach. The Markov dependence structure
in model C clearly challenged the marginal approach. For model D, none of
the approaches controlled the type I error.
To determine the sensitivity of the approach to using distinct data for each
stage of the analysis, simulations were conducted screening on the first half
of the data and cleaning on the second half (2 splits). The tuning parameter
was selected using leave-one-out cross validation (Table 1, bottom half). As
expected, this approach lead to a dramatic increase in the power of all the
procedures. More surprising is the fact that the type I error was near α or
12 L. WASSERMAN AND K. ROEDER
below for models A, B and C. Clearly this approach has advantages over
data splitting and merits further investigation.
A natural competitor to screen and clean procedure is a two-stage adap-
tive lasso [Zou (2006)]. In our implementation, we split the data and used
one half for each stage of the analysis. At stage one, leave-one-out cross val-
idation lasso screens the data. In stage two, the adaptive lasso, with weights
wj = |β̂j |−1, cleans the data. The tuning parameter for the lasso was again
chosen using leave-one-out cross validation. Table 2 provides the size, power
and false positive rate (FPR) for this procedure. Naturally, the adaptive
lasso does not control the size of the test, but the FPR is small. The power
of the test is greater than we found for our lasso screen and clean procedure,
but this extra power comes at the cost of a much higher type I error rate.
7. Proofs. Recall that if A is a square matrix, then φ(A) and Φ(A)
denote the smallest and largest eigenvalues of A. Throughout the proofs we
make use of the following fact. If v is a vector and A is a square matrix,
φ(A)‖v‖2 ≤ vTAv ≤Φ(A)‖v‖2.(14)
We use the following standard tail bound: if Z ∼N(0,1), then P(|Z|> t)≤
t−1e−t
2/2. We will also use the following results about the lasso from Mein-
shausen and Yu (2008). Their results are stated and proved for fixed X but,
under the conditions (A1)–(A5), it is easy to see that their conditions hold
with probability tending to one and so their results hold for random X as
well.
Theorem 7.1 [Meinshausen and Yu (2008)]. Let β̃(λ) be the lasso esti-
mator.
Table 2
Size, power and false positive rate (FPR) of two-stage adaptive lasso procedure
n p Model Size Power FPR
100 100 A 0.93 0 0.032
100 100 B 0.84 0.97 0.034
100 100 C 0.81 0.96 0.031
100 10 D 0.67 0.21 0.114
100 1000 A 0.96 0 0.004
100 1000 B 0.89 0.65 0.004
100 1000 C 0.76 0.77 0.002
1000 10 D 0.73 0.24 0.013
VARIABLE SELECTION 13
1. The squared error satisfies
‖β̃(λ)− β‖22 ≤
cm log pn
nφ2n(m)
→ 1,(15)
where m= |Ŝn(λ)| and c > 0 is a constant.
2. The size of Ŝn(λ) satisfies
|Ŝn(λ)| ≤
τ2Cn2
→ 1,(16)
where τ2 = E(Y 2i ).
Proof of Lemma 3.1. Let D ⊂M and φ= φ(n−1XTMXM ). Then,
L(β̂M ) =
εTXM (X
MXM )
−1XTMε≤
‖XTMε‖2 =
Z2j ,
where Zj = n
−1/2XT•jε. Conditional on X , Zi ∼N(0, a2j ) where a2j = n−1 ×∑n
ij . Let A
n =max1≤j≤pn a
j . By Hoeffding’s inequality, (A2) and (A5),
P(En)→ 1 where En = {An ≤
2}. So
1≤j≤pn
|Zj |>
4 log pn
1≤j≤pn
|Zj |>
4 log pn,En
1≤j≤pn
|Zj |>
4 log pn,E
1≤j≤pn
|Zj |>
4 log pn,En
+ P(Ecn)
An max
1≤j≤pn
|Zj |
4 log pn,En
+ o(1)
1≤j≤pn
|Zj |
2 log pn
+ o(1)
1≤j≤pn
|Zj |
2 log pn
)∣∣∣X
+ o(1)
2 log pn
+ o(1) = o(1).
j∈M Z
j ≤mmax1≤j≤pn Z2j and (6) follows.
Now we lower bound L(β̂M ). Let M be such that D 6⊂ M . Let A =
{j : β̂M (j) 6= 0} ∪D. Then, |A| ≤m + s. Therefore, with probability tend-
ing to 1,
L(β̂M ) =
(β̂M − β)TXTX(β̂M − β) =
(β̂M − β)TXTAXA(β̂M − β)
14 L. WASSERMAN AND K. ROEDER
≥ φn(m+ s)‖β̂M − β‖2 = φn(m+ s)
(β̂M (j)− β(j))2
≥ φn(m+ s)
j∈D∩Mc
(0− β(j))2 ≥ φn(m+ s)ψ2.
Proof of Theorem 3.2. Let Ỹ denote the responses, and X̃ the design
matrix, for the second half of the data. Then, Ỹ = X̃β + ε̃. Now
L(λ) =
(β̂(λ)− β)TXTX(β̂(λ)− β) = (β̂(λ)− β)T Σ̂n(β̂(λ)− β)
L̂(λ) = n−1‖Ỹ −X̃β̂(λ)‖2 = (β̂(λ)−β)T Σ̃n(β̂(λ)−β)+δn+
〈ε̃, X̃(β̂(λ)−β)〉,
where δn = ‖ε̃‖2/n, and Σ̂n = n−11 XTX and Σ̃n = n−1X̃T X̃ . By Hoeffding’s
inequality
P(|Σ̂n(j, k)− Σ̃n(j, k)|> ε)≤ e−ncε
for some c > 0, and so
|Σ̂n(j, k)− Σ̃n(j, k)|> ε
≤ p2ne−ncε
Choose εn = 4/(cn
1−c2). It follows that
|Σ̂n(j, k)− Σ̃n(j, k)|>
cn1−c2
≤ e−2nc2 → 0.
Note that
|{j : β̂j(λ) 6= 0} ∪ {j :βj 6= 0}| ≤ kn + s.
Hence, with probability tending to 1
|L(λ)− L̂(λ)− δn| ≤
cn1−c2
‖β̂(λ)− β‖21 + 2ξn(λ)
for all λ ∈Λn, where
ξn(λ) =
ε̃iµi(λ)
and µi(λ) = X̃
i (β̂(λ)−β). Now ‖β̂(λ)−β‖21 =OP ((kn+s)2) since ‖β̂(λ)‖2 =
OP (kn/φ(kn)). Thus, ‖β̂(λ)−β‖1 ≤C(kn+s) with probability tending to 1,
for some C > 0. Also, |µi(λ)| ≤B‖β̂(λ)−β‖1 ≤BC(kn+ s) with probability
tending to 1. Let W ∼N(0,1). Conditional on D1,
|ξn(λ)|
µ2i (λ)|W | ≤
BC(kn + s)|W |,
VARIABLE SELECTION 15
so supλ∈Λn |ξn(λ)|=OP (kn/
n). �
Proof of Theorem 4.1. 1. Let λn = τn
C/kn,M = Ŝn(λn) and m=
|M |. Then, P(m≤ kn)→ 1 due to (16). Hence, P(λn ∈ Λn)→ 1. From (15),
‖β̃(λn)− β‖22 ≤O
kn log pn
= oP (1).
Hence, ‖β̃(λn)− β‖2∞ = oP (1). So, for each j ∈D,
|β̃j(λn)| ≥ |βj | − |β̃j(λn)− βj | ≥ ψ+ oP (1)
and hence, P(minj∈D |β̃j(λn)|> 0)→ 1. Therefore, Γn = {λ ∈Λn :D ⊂ Ŝn(λ)}
is nonempty. By Lemma 3.1,
L(λn)≤ cm log pn/(nφ(m)) =OP (kn log pn/n).(17)
On the other hand, from Lemma 3.1,
λ∈Λn∩Γcn
L(β̂λ)>ψ
2φ(kn)
→ 1.(18)
Now, nφn(kn)/(kn log pn)→∞, and so (17) and (18) imply that
λ∈Λn∩Γcn
L(β̂λ)>L(λn)
Thus, if λ∗ denotes the minimizer of L(λ) over Λn, we conclude that P(λ∗ ∈
Γn)→ 1, and hence, P(D ⊂ Ŝn(λ∗))→ 1.
2. This follows from part 1 and Theorem 3.2.
3. Let A= Ŝn ∩Dc. We want to show that
|Tj|> cn
≤ α+ o(1).
|Tj |> cn
|Tj |> cn,D ⊂ Ŝn
|Tj |> cn,D 6⊂ Ŝn
|Tj |> cn,D ⊂ Ŝn
+ P(D 6⊂ Ŝn)
|Tj |> cn,D ⊂ Ŝn
+ o(1).
Conditional on (D1,D2), β̂A is normally distributed with mean 0 and vari-
ance matrix σ2(XTAXA)
−1 when D⊂ Ŝn. Recall that
Tj(M) =
eTj (X
MXM )
−1XTMY
eTj (X
MXM )
β̂M,j
16 L. WASSERMAN AND K. ROEDER
whereM = Ŝn, s
j = σ̂
2eTj (X
MXM )
−1ej and ej = (0, . . . ,0,1,0, . . . ,0)
T , where
the 1 is in the jth coordinate. When D ⊂ Ŝn, each Tj , for j ∈ A, has a t-
distribution with n−m degrees of freedom wherem= |Ŝn|. Also, cn/tα/2m→
1 where tu denotes the upper tail critical value for the t-distribution. Hence,
|Tj |> cn,D ⊂ Ŝn|D1,D2
|Tj |> tα/2m,D ⊂ Ŝn|D1,D2
≤ α+ an,
where an = o(1), since |A| ≤m. It follows that
|Tj |> cn,D ⊂ Ŝn
≤ α+ o(1).
Proof of Theorem 4.2. From Theorem 4.1, P(D̂n ∩Dc 6= ∅) ≤ αn
and so P(D̂n ∩Dc 6=∅)→ 0. Hence, P(D̂n ⊂D)→ 1. It remains to be shown
P(D⊂ D̂n)→ 1.(19)
The test statistic for testing βj = 0 when Ŝn =M is
Tj(M) =
eTj (X
MXM )
−1XTMY
eTj (X
MXM )
For simplicity in the proof, let us take σ̂ = σ, the extension to unknown σ
being straightforward. Let j ∈D,M= {M : |M | ≤ kn,D⊂M}. Then,
P(j /∈ D̂n) = P(j /∈ D̂n,D ⊂ Ŝn) + P(j /∈ D̂n,D 6⊂ Ŝn)
≤ P(j /∈ D̂n,D ⊂ Ŝn) + P(D 6⊂ Ŝn)
= P(j /∈ D̂n,D ⊂ Ŝn) + o(1)
P(j /∈ D̂n, Ŝn =M) + o(1)
P(|Tj(M)|< cn, Ŝn =M) + o(1)
P(|Tj(M)|< cn) + o(1).
Conditional on D1 ∪D2, for each M ∈M, Tj(M) = (βj/sj) +Z, where Z ∼
N(0,1). Without loss of generality, assume that βj > 0. Hence,
P(|Tj(M)|< cn|D1 ∪D2) = P
−cn −
<Z < cn −
VARIABLE SELECTION 17
Fix a small ε > 0. Note that s2j ≤ σ2/(nκ). It follows that, for all large n,
cn − βj/sj <−ε
n. So,
P(|Tj(M)|< cn|D1 ∪D2)≤ P(Z <−ε
n)≤ e−nε2/2.
The number of models inM is
pn − s
j − s
pn − s
kn − s
(pn − s)e
kn − s
)kn−s
≤ knpknn ,
where we used the inequality
P(|Tj(M)|< cn|D1 ∪D2)≤ knpknn e−nε
2 → 0
by (A2). We have thus shown that P(j /∈ D̂n)→ 0, for each j ∈D. Since |D|
is finite, it follows that P(j /∈ D̂n for some j ∈D)→ 0 and hence (19). �
Proof of Theorem 4.5. A simple modification of Theorem 3.1 of
Barron et al. (2008) shows that
L(kn) =
‖Ŷkn −Xβ‖2 = oP (1).
[The modification is needed because Barron et al. (2008) require Y to be
bounded while we have assumed that Y is normal. By a truncation argument,
we can still derive the bound on L(kn).] So
‖β̂kn − β‖2 ≤
L(kn)
φn(kn + s)
L(kn)
= oP (1).
Hence, for any ε > 0, with probability tending to 1, ‖β̂(kn)−β‖2 < ε so that
|β̂j | > ψ/2 > 0, for all j ∈D. Thus, P(D ⊂ Ŝn(kn))→ 1. The remainder of
the proof of part 1 is the same as in Theorem 4.1. Part 2 follows from the
previous result together with Theorem 3.2. The proof of part 3 is the same
as for Theorem 4.1. �
Proof of Theorem 4.6. Note that µ̂j − µj = n−1
i=1Xijεi. Hence,
µ̂j − µj ∼N(0,1/n). So, for any δ > 0,
|µ̂j − µj |> δ
P(|µ̂j − µj|> δ)
2/2 ≤ c1e
2/2→ 0.
18 L. WASSERMAN AND K. ROEDER
By (12), conclude that D ⊂ Ŝn(λ) when λ = µ̂(kn). The remainder of the
proof is the same as the proof of Theorem 4.5. �
Proof of Theorem 4.7. Let A= Ŝn ∩Dc. We want to show that
|Tj|> cn
≤ α+ o(1).
For fixed A, β̂A is normal with mean 0 but this is not true for random A.
Instead we need to bound Tj . Recall that
Tj(M) =
eTj (X
MXM )
−1XTMY
eTj (X
MXM )
β̂M,j
whereM = Ŝn, s
j = σ̂
2eTj (X
MXM )
−1ej and ej = (0, . . . ,0,1,0, . . . ,0)
T where
the 1 is in the jth coordinate. The probabilities that follow are conditional
on D1 but this is supressed for notational convenience. First, write
|Tj |> cn
|Tj |> cn,D ⊂ Ŝn
|Tj |> cn,D 6⊂ Ŝn
|Tj |> cn,D ⊂ Ŝn
+ P(D 6⊂ Ŝn)
|Tj |> cn,D ⊂ Ŝn
+ o(1).
When D ⊂ Ŝn,
)−1 1
where Q
= ((1/n)XT
)−1, γ
= n−1XT
ε, and β
(j) = 0, for j ∈A.
Now, s2j ≥ σ̂2/(nC) so that
|Tj | ≤
nC|β̂
Ŝn,j
n log logn|β̂
Ŝn,j
for j ∈ Ŝn. Therefore,
|Tj |> cn,D ⊂ Ŝn
Ŝn,j
σ̂cn√
,D ⊂ Ŝn
Let γ = n−1XT ε. Then,
‖β̂A‖2 ≤ γT
knmax1≤j≤pn γ
VARIABLE SELECTION 19
It follows that
Ŝn,j
knmax1≤j≤pn |γj |
kn log logn max
1≤j≤pn
|γj |,
since κ > 0. So,
Ŝn,j
|> σ̂cn√
n log logn
,D ⊂ Ŝn
1≤j≤pn
|γj |>
log logn
Note that γj ∼N(0, σ2/n), and hence,
|γj |
2σ2 log(2pn)
There exists εn → 0 such that P(Bn)→ 1 where Bn = {(1 − εn) ≤ σ̂/σ ≤
(1 + ε)}. So,
1≤j≤pn
|γj|>
log logn
1≤j≤pn
|γj |>
σcn(1− εn)
log logn
σ(1− εn)cn
log logn
|γj |
≤ α+ o(1). �
8. Discussion. The multi-stage method presented in this paper success-
fully controls type I error while giving reasonable power. The lasso and step-
wise have similar performance. Although theoretical results assume indepen-
dent data for each of the three stages, simulations suggest that leave-one-out
cross-validation leads to valid type I error rates and greater power. Screening
the data in one phase of the experiment and cleaning in a followup phase
leads to an efficient experimental design. Certainly this approach deserves
further theoretical investigation. In particular, the question of optimality is
an open question.
The literature on high-dimensional variable selection is growing quickly.
The most important deficiency in much of this work, including this paper,
is the assumption that the model Y = XTβ + ε is correct. In reality, the
model is at best an approximation. It is possible to study linear procedures
when the linear model is not assumed to hold as in Greenshtein and Ritov
(2004). We discuss this point in the Appendix. Nevertheless, it seems useful
to study the problem under the assumption of linearity to gain insight into
these methods. Future work should be directed at exploring the robustness
of the results when the model is wrong.
Other possible extensions include: dropping the normality of the errors,
permitting nonconstant variance, investigating the optimal sample sizes for
each stage and considering other screening methods besides cross-validation.
20 L. WASSERMAN AND K. ROEDER
Finally, let us note that the example involving unfaithfulness, that is,
cancellations of parameters to make the marginal correlation much different
than the regression coefficient, pose a challenge for all the methods and
deserve more attention even in cases of small p.
APPENDIX: PREDICTION
Realistically, there is little reason to believe that the linear model is cor-
rect. Even if we drop the assumption that the linear model is correct, sparse
methods like the lasso can still have good properties as shown in Greenshtein
and Ritov (2004). In particular, they showed that the lasso satisfies a risk
consistency property. In this appendix we show that this property continues
to hold if λ is chosen by cross-validation.
The lasso estimator is the minimizer of
i=1(Yi −XTi β)2 + λ‖β‖1. This
is equivalent to minimizing
i=1(Yi−XTi β)2 subject to ‖β‖1 ≤Ω, for some
Ω. (More precisely, the set of estimators as λ varies is the same as the set of
estimators as Ω varies.) We use this second version throughout this section.
The predictive risk of a linear predictor ℓ(x) = xTβ is R(β) = E(Y −ℓ(x))2
where (X,Y ) denotes a new observation. Let γ = γ(β) = (−1, β1, . . . , βp)T
and let Γ = E(ZZT ) where Z = (Y,X1, . . . ,Xp). Then, we can write R(β) =
γTΓγ. The lasso estimator can now be written as β̂(Ωn) = argminβ∈B(Ωn) R̂(β)
where R̂(β) = γT Γ̂γ and Γ̂ = n−1
i=1ZiZ
Define
β∗ = argmin
β∈B(Ωn)
R(β),
where
B(Ωn) = {β :‖β‖1 ≤Ωn}.
Thus, ℓ∗(x) = x
Tβ∗ is the best linear predictor in the set B(Ωn). The best
linear predictor is well defined even though E(Y |X) is no longer assumed to
be linear. Greenshtein and Ritov (2004) call an estimator β̂n persistent, or
predictive risk consistent, if
R(β̂n)−R(β∗)
as n→∞.
The assumptions we make in this section are:
(B1) pn ≤ en
, for some 0≤ ξ < 1;
(B2) the elements of Γ̂ satisfy an exponential inequality
P(|Γ̂jk − Γjk|> ε)≤ c3e−nc4ε
for some c3, c4 > 0;
VARIABLE SELECTION 21
(B3) there exists B0 <∞ such that, for all n, maxj,kE(|ZjZk|)≤B0.
Condition (A2) can easily be deduced from more primitive assumptions
as in Greenshtein and Ritov (2004), but for simplicity we take (A2) as an
assumption. Let us review one of the results in Greenshtein and Ritov (2004).
For the moment, replace (A1) with the assumption that pn ≤ nb, for some
b. Under these conditions, it follows that
∆n ≡max
|Γ̂jk − Γjk|=OP
Hence,
β∈B(Ωn)
|R(β)− R̂(β)|= sup
β∈B(Ωn)
|γT (Γ− Γ̂)γ|
≤∆n sup
β∈B(Ωn)
‖γ‖21 =Ω2nOP
The latter term is oP (1) as long as Ωn = o((n/ logn)
1/4). Thus, we have the
following.
Theorem A.1 [Greenshtein and Ritov (2004)]. If Ωn = o((n/ logn)
1/4),
then the lasso estimator is persistent.
For future reference, let us state a slightly different version of their result
that we will need. We omit the proof.
Theorem A.2. Let γ > 0 be such that ξ+γ < 1. Let Ωn =O(n
(1−ξ−γ)/4).
Then, under (B1) and (B2),
β∈B(Ωn)
|R̂(β)−R(β)|> 1
=O(e−cn
)(20)
for some c > 0.
The estimator β̂(Ωn) lies on the boundary of the ball B(Ωn) and is very
sensitive to the exact choice of Ωn. A potential improvement—and something
that reflects actual practice—is to compute the set of lasso estimators β̂(ℓ),
for 0≤ ℓ≤ Ωn and then select from that set based on cross validation. We
now confirm that the resulting estimator preserves persistence. As before,
we split the data into D1 and D2. Construct the lasso estimators {β̂(ℓ) : 0≤
ℓ≤Ωn}. Choose ℓ̂ by cross validation using D2. Let β̂ = β̂(ℓ̂).
22 L. WASSERMAN AND K. ROEDER
Theorem A.3. Let γ > 0 be such that ξ + γ < 1. Under (A1), (A2)
and (A3), if Ωn =O(n
(1−ξ−γ)/4), then the cross-validated lasso estimator β̂
is persistent. Moreover,
R(β̂)− inf
0≤ℓ≤Ωn
R(β̂(ℓ))
P→ 0.(21)
Proof. Let β∗(ℓ) = argminβ∈B(ℓ)R(β). Define h(ℓ) = R(β∗(ℓ)), g(ℓ) =
R(β̂(ℓ)) and c(ℓ) = L̂(β̂(ℓ)). Note that, for any vector b, we can write R(b) =
τ2 + bTΣb− 2bTρ where ρ= (E(Y X1), . . . ,E(Y Xp))T .
Clearly, h is monotone nonincreasing on [0,Ωn]. We claim that |h(ℓ +
δ) − h(ℓ)| ≤ cΩnδ where c depends only on Γ. To see this, let u = β∗(ℓ),
v = β∗(ℓ+ δ) and a= ℓβ∗(ℓ+ δ)/(ℓ+ δ) so that a ∈B(ℓ). Then,
h(ℓ+ δ)≤ h(ℓ)
=R(u)≤R(a)
=R(v) +R(a)−R(v)
= h(ℓ+ δ) +
ρT v− δ(2ℓ+ δ)
(ℓ+ δ)2
≤ h(ℓ+ δ) + 2δC + δ(2Ωn + δ)C,
where C =maxj,k |Γj,k|=O(1).
Next, we claim that g(ℓ) is Lipschitz on [0,Ωn] with probability tending to
1. Let β̂(ℓ) = argminβ∈B(ℓ) R̂(β) denote the lasso estimator and set û= β̂(ℓ)
and v̂ = β̂(ℓ+δ). Let εn = n
−γ/4. From (20), the following chain of equations
hold except on a set of exponentially small probability:
g(ℓ+ δ) =R(v̂)≤ R̂(v̂) + εn ≤ R̂(v) + εn
≤R(v) + 2εn = h(ℓ+ δ) + 2εn
≤ h(ℓ) + cΩnδ +2εn =R(u) + cΩnδ+ 2εn
≤R(û) + cΩnδ+2εn = g(ℓ) + cΩnδ +2εn.
A similar argument can be applied in the other direction. Conclude that
|g(ℓ+ δ)− g(ℓ)| ≤ cΩnδ +2εn
except on a set of small probability.
Now let A= {0, δ,2δ, . . . ,mδ} where m is the smallest integer such that
mδ ≥ Ωn. Thus, m∼ Ωn/δn. Choose δ = δn = n−3(1−ξ−γ)/8. Then, Ωnδn→
0 and Ωn/δn ≤ n3(1−ξ−γ)/4. Using the same argument as in the proof of
Theorem 3.2,
|L̂(β̂(ℓ))−R(β̂(ℓ))|= σn,
VARIABLE SELECTION 23
where σn = oP (1). Then,
R(β∗(Ωn))≤R(β̂)≤ L̂(β̂(ℓ̂)) + σn
≤ L̂(mδn) + σn ≤ g(mδn) + 2σn ≤ g(Ωn) + 2σn + cΩnδn
≤ h(Ωn) + 2σn + εn + cΩnδn
=R(β∗(Ωn)) + 2σn + εn + cΩnδn
and persistence follows. To show the second result, let β̃ = argmin0≤ℓ≤Ωn g(ℓ)
and β = argminℓ∈A g(ℓ). Then,
R(β̃)≤ L̂(β̃) + σn ≤ L̂(β) + σn
≤R(β) + 2σn ≤R(β̃) + 2σn + cδnΩn
and the claim follows. �
Acknowledgments. The authors are grateful for the use of a portion of
the sample from the Osteoporotic Fractures in Men (MrOS) Study to il-
lustrate their methodology. MrOs is supported by the National Institute of
Arthritis and Musculoskeletal and Skin Diseases (NIAMS), the National In-
stitute on Aging (NIA) and the National Cancer Institute (NCI) through
Grants U01 AR45580, U01 AR45614, U01 AR45632, U01 AR45647, U01
AR45654, U01 AR45583, U01 AG18197 and M01 RR000334. Genetic analy-
ses in MrOS were supported by R01-AR051124. We also thank two referees
and an AE for helpful suggestions.
REFERENCES
Barron, A., Cohen, A., Dahmen, W. and DeVore, R. (2008). Approximation and
learning by greedy algorithms. Ann. Statist. 36 64–94. MR2387964
Bühlmann, P. (2006). Boosting for high-dimensional linear models. Ann. Statist. 34 559–
583. MR2281878
Candes, E. and Tao, T. (2007). The Dantzig selector: Statistical estimation when p is
much larger than n. Ann. Statist. 35 2313–2351. MR2382644
Donoho, D. (2006). For most large underdetermined systems of linear equations, the
minimal l1-norm near-solution approximates the sparsest near-solution. Comm. Pure
Appl. Math. 59 797–829. MR2217606
Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004). Least angle regres-
sion. Ann. Statist. 32 407–499. MR2060166
Fan, J. and Lv, J. (2008). Sure independence screening for ultra-high-dimensional feature
space. J. Roy. Statist. Soc. Ser. B 70 849–911.
Greenshtein, E. and Ritov, Y. (2004). Persistence in high-dimensional linear predictor
selection and the virtue of overparametrization. Bernoulli 10 971–988. MR2108039
Meinshausen, N. (2007). Relaxed lasso. Comput. Statist. Data Anal. 52 374–393.
MR2409990
Meinshausen, N. and Bühlmann, P. (2006). High dimensional graphs and variable se-
lection with the lasso. Ann. Statist. 34 1436–1462. MR2278363
http://www.ams.org/mathscinet-getitem?mr=2387964
http://www.ams.org/mathscinet-getitem?mr=2281878
http://www.ams.org/mathscinet-getitem?mr=2382644
http://www.ams.org/mathscinet-getitem?mr=2217606
http://www.ams.org/mathscinet-getitem?mr=2060166
http://www.ams.org/mathscinet-getitem?mr=2108039
http://www.ams.org/mathscinet-getitem?mr=2409990
http://www.ams.org/mathscinet-getitem?mr=2278363
24 L. WASSERMAN AND K. ROEDER
Meinshausen, N. and Yu, B. (2009). Lasso-type recovery of sparse representations of
high-dimensional data. Ann. Statist. To appear.
Orwoll, E., Blank, J. B., Barrett-Connor, E., Cauley, J., Cummings, S., En-
srud, K., Lewis, C., Cawthon, P. M., Marcus, R., Marshall, L. M., McGowan,
J., Phipps, K., Sherman, S., Stefanick, M. L. and Stone, K. (2005). Design and
baseline characteristics of the osteoporotic fractures in men (MrOS) study—a large ob-
servational study of the determinants of fracture in older men. Contemp. Clin. Trials.
26 569–585.
Robins, J., Scheines, R., Spirtes, P. and Wasserman, L. (2003). Uniform consistency
in causal inference. Biometrika 90 491–515. MR2006831
Spirtes, P., Glymour, C. and Scheines, R. (2001). Causation, Prediction, and Search.
MIT Press. MR1815675
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Roy. Statist.
Soc. Ser. B 58 267–288. MR1379242
Tropp, J. A. (2004). Greed is good: Algorithmic results for sparse approximation. IEEE
Trans. Inform. Theory 50 2231–2242. MR2097044
Tropp, J. A. (2006). Just relax: Convex programming methods for identifying sparse
signals in noise. IEEE Trans. Inform. Theory 52 1030–1051. MR2238069
Wainwright, M. (2006). Sharp thresholds for high-dimensional and noisy recovery of
sparsity. Available at arxiv.org/math.ST/0605740.
Wellcome Trust (2007). Genome-wide association study of 14,000 cases of seven com-
mon diseases and 3000 shared controls. Nature 447 661–678.
Zhao, P. and Yu, B. (2006). On model selection consistency of lasso. J. Mach. Learn.
Res. 7 2541–2563. MR2274449
Zhang, C. H. and Huang, J. (2006). Model selection consistency of the lasso in high-
dimensional linear regression. Ann. Statist. 36 1567–1594. MR2435448
Zou, H. (2006). The adaptive lasso and its oracle properties. J. Amer. Statist. Assoc. 101
1418–1429. MR2279469
Department of Statistics
Carnegie Mellon University
Pittsburgh
E-mail: larry@stat.cmu.edu
roeder@stat.cmu.edu
http://www.ams.org/mathscinet-getitem?mr=2006831
http://www.ams.org/mathscinet-getitem?mr=1815675
http://www.ams.org/mathscinet-getitem?mr=1379242
http://www.ams.org/mathscinet-getitem?mr=2097044
http://www.ams.org/mathscinet-getitem?mr=2238069
http://arxiv.org/math.ST/0605740
http://www.ams.org/mathscinet-getitem?mr=2274449
http://www.ams.org/mathscinet-getitem?mr=2435448
http://www.ams.org/mathscinet-getitem?mr=2279469
mailto:larry@stat.cmu.edu
mailto:roeder@stat.cmu.edu
	Introduction
	Error control
	Loss and cross-validation
	Loss
	Cross-validation
	Multi-stage methods
	The lasso
	Stepwise regression
	Marginal regression
	Modifications
	Application
	Simulations
	Proofs
	Discussion
	Appendix: Prediction
	Acknowledgments
	References
	Author's addresses

  This paper explores the following question: what kind of statistical
guarantees can be given when doing variable selection in high-dimensional
models? In particular, we look at the error rates and power of some multi-stage
regression methods. In the first stage we fit a set of candidate models. In the
second stage we select one model by cross-validation. In the third stage we use
hypothesis testing to eliminate some variables. We refer to the first two
stages as "screening" and the last stage as "cleaning." We consider three
screening methods: the lasso, marginal regression, and forward stepwise
regression. Our method gives consistent variable selection under certain
conditions.

Introduction. Several methods have been developed lately for high-
dimensional linear regression such as the lasso [Tibshirani (1996)], Lars
[Efron et al. (2004)] and boosting [Bühlmann (2006)]. There are at least
two different goals when using these methods. The first is to find models
with good prediction error. The second is to estimate the true “sparsity
pattern,” that is, the set of covariates with nonzero regression coefficients.
These goals are quite different and this paper will deal with the second goal.
(Some discussion of prediction is in the Appendix.) Other papers on this
topic include Meinshausen and Bühlmann (2006), Candes and Tao (2007),
Wainwright (2006), Zhao and Yu (2006), Zou (2006), Fan and Lv (2008),
Meinshausen and Yu (2008), Tropp (2004, 2006), Donoho (2006) and Zhang
and Huang (2006). In particular, the current paper builds on ideas in Mein-
shausen and Yu (2008) and Meinshausen (2007).
Let (X1, Y1), . . . , (Xn, Yn) be i.i.d. observations from the regression model
Yi =X
i β + εi,(1)
Received June 2008; revised August 2008.
1Supported by NIH Grant MH057881.
AMS 2000 subject classifications. Primary 62J05; secondary 62J07.
Key words and phrases. Lasso, stepwise regression, sparsity.
This is an electronic reprint of the original article published by the
Institute of Mathematical Statistics in The Annals of Statistics,
2009, Vol. 37, No. 5A, 2178–2201. This reprint differs from the original in
pagination and typographic detail.
http://arxiv.org/abs/0704.1139v2
http://www.imstat.org/aos/
http://dx.doi.org/10.1214/08-AOS646
http://www.imstat.org
http://www.ams.org/msc/
http://www.imstat.org
http://www.imstat.org/aos/
http://dx.doi.org/10.1214/08-AOS646
2 L. WASSERMAN AND K. ROEDER
where ε ∼ N(0, σ2), Xi = (Xi1, . . . ,Xip)T ∈ Rp and p = pn > n. Let X be
the n × p design matrix with jth column X•j = (X1j , . . . ,Xnj)T and let
Y = (Y1, . . . , Yn)
T . Let
D= {j :βj 6= 0}
be the set of covariates with nonzero regression coefficients. Without loss
of generality, assume that D = {1, . . . , s} for some s. A variable selection
procedure D̂n maps the data into subsets of S = {1, . . . , p}.
The main goal of this paper is to derive a procedure D̂n such that
lim sup
P(D̂n ⊂D)≥ 1− α,(2)
that is, the asymptotic type I error is no more than α. Note that throughout
the paper we use ⊂ to denote nonstrict set-inclusion. Moreover, we want
D̂n to have nontrivial power. Meinshausen and Bühlmann (2006) control a
different error measure. Their method guarantees lim supn→∞P(D̂n ∩ V 6=
∅)≤ α where V is the set of variables not connected to Y by any path in
an undirected graph.
Our procedure involves three stages. In stage I we fit a suite of candidate
models, each model depending on a tuning parameter λ,
S = {Ŝn(λ) :λ ∈ Λ}.
In stage II we select one of those models Ŝn using cross-validation to select λ̂.
In stage III we eliminate some variables by hypothesis testing. Schematically,
stage I−→ S stage II−→︸ ︷︷ ︸
screen
stage III−→ D̂n︸ ︷︷ ︸
clean
Genetic epidemiology provides a natural setting for applying screen and
clean. Typically, the number of subjects, n, is in the thousands, while p
ranges from tens of thousands to hundereds of thousands of genetic fea-
tures. The number of genes exhibiting a detectable association with a trait
is extremely small. Indeed, for type I diabetes only ten genes have exhibited
a reproducible signal [Wellcome Trust (2007)]. Hence, it is natural to assume
that the true model is sparse. A common experimental design involves a 2-
stage sampling of data, with stages 1 and 2 corresponding to the screening
and cleaning processes, respectively.
In stage 1 of a genetic association study, n1 subjects are sampled and one
or more traits such as bone mineral density are recorded. Each subject is
also measured at p locations on the chromosomes. These genetic covariates
usually have two forms in the population due to variability at a single nu-
cleotide and hence are called single nucleotide polymorphisms (SNPs). The
distinct forms are called alleles. Each covariate takes on a value (0, 1 or 2)
VARIABLE SELECTION 3
indicating the number of copies of the less common allele observed. For a
well-designed genetic study, individual SNPs are nearly uncorrelated unless
they are physically located in very close proximity. This feature makes it
much easier to draw causal inferences about the relationship between SNPs
and quantitative traits. It is standard in the field to infer that an association
discovered between a SNP and a quantitative trait implies a causal genetic
variant is physically located near the one exhibiting association. In stage 2,
n2 subjects are sampled at a subset of the SNPs assessed in stage 1. SNPs
measured in stage 2 are often those that achieved a test statistic that ex-
ceeded a predetermined threshold of significance in stage 1. In essence, the
two stage design pairs naturally with a screen and clean procedure.
For the screen and clean procedure, it is essential that Ŝn has two prop-
erties as n→∞ as follows:
P(D ⊂ Ŝn)→ 1(3)
|Ŝn|= oP (n),(4)
where |M | denotes the number of elements in a setM . Condition (3) ensures
the validity of the test in stage III while condition (4) ensures that the power
of the test is not too small. Without condition (3), the hypothesis test in
stage III would be biased. We will see that the power goes to 1, so taking
α = αn→ 0 implies consistency: P(D̂n =D)→ 1. For fixed α, the method
also produces a confidence sandwich for D, namely,
lim inf
P(D̂n ⊂D ⊂ Ŝn)≥ 1−α.
To fit the suite of candidate models, we consider three methods. In method
Ŝn(λ) = {j : β̃j(λ) 6= 0},
where β̃j(λ) is the lasso estimator, the value of β that minimizes
(Yi −XTi β)2 + λ
|βj |.
In method 2, take Ŝn(λ) to be the set of variables chosen by forward stepwise
regression after λ steps. In method 3, marginal regression, we take
Ŝn = {j : |µ̂j|> λ},
where µ̂j is the marginal regression coefficient from regressing Y onXj . (This
is equivalent to ordering by the absolute t-statistics since we will assume that
the covariates are standardized.) These three methods are very similar to
basis pursuit, orthogonal matching pursuit and thresholding [see, e.g., Tropp
(2004, 2006) and Donoho (2006)].
4 L. WASSERMAN AND K. ROEDER
Notation. Let ψ =minj∈D |βj |. Define the loss of any estimator β̂ by
L(β̂) =
(β̂ − β)TXTX(β̂ − β) = (β̂ − β)T Σ̂n(β̂ − β),(5)
where Σ̂n = n
−1XTX . For convenience, when β̂ ≡ β̂(λ) depends on λ we
write L(λ) instead of L(β̂(λ)). If M ⊂ S, let XM be the design matrix
with columns (X•j : j ∈M) and let β̂M = (XTMXM )−1XTMY denote the least-
squares estimator, assuming it is well defined. Note that our use ofX•j differs
from standard ANOVA notation. Write Xλ instead of XM whenM = Ŝn(λ).
When convenient, we extend β̂M to length p by setting β̂M (j) = 0, for j /∈M .
We use the norms
v2j , ‖v‖1 =
|vj | and ‖v‖∞ =max
|vj |.
If C is any square matrix, let φ(C) and Φ(C) denote the smallest and
largest eigenvalues of C. Also, if k is an integer, define
φn(k) = min
M : |M |=k
XTMXM
and Φn(k) = max
M : |M |=k
XTMXM
We will write zu for the upper quantile of a standard normal, so that P(Z >
zu) = u where Z ∼N(0,1).
Our method will involve splitting the data randomly into three groups
D1, D2 and D3. For ease of notation, assume the total sample size is 3n and
that the sample size of each group is n.
Summary of assumptions. We will use the following assumptions through-
out except in Section 8:
(A1) Yi =X
i β + εi where εi ∼N(0, σ2), for i= 1, . . . , n.
(A2) The dimension pn of X satisfies pn→∞ and pn ≤ c1en
, for some
c1 > 0 and 0≤ c2 < 1.
(A3) s≡ |{j :βj 6= 0}|=O(1) and ψ =min{|βj | :βj 6= 0}> 0.
(A4) There exist positive constants C0,C1 and κ such that
P(lim supn→∞Φn(n)≤C0) = 1 and P(lim infn→∞ φn(C1 logn)≥ κ) = 1. Also,
P(φn(n)> 0) = 1, for all n.
(A5) The covariates are standardized: E(Xij) = 0 and E(X
ij) = 1. Also,
there exists 0<B <∞ such that P(|Xjk| ≤B) = 1.
For simplicity, we include no intercepts in the regressions. The assump-
tions can be weakened at the expense of more complicated proofs. In par-
ticular, we can let s increase with n and ψ decrease with n. Similarly, the
normality and constant variance assumptions can be relaxed.
VARIABLE SELECTION 5
2. Error control. Define the type I error rate q(D̂n) = P(D̂n ∩ Dc 6=
∅) and the asymptotic error rate lim supn→∞ q(D̂n). We define the power
π(D̂n) = P(D ⊂ D̂n) and the average power
πav =
P(j ∈ D̂n).
It is well known that controlling the error rate is difficult for at least
three reasons: correlation of covariates, high-dimensionality of the covariate
and unfaithfulness (cancellations of correlations due to confounding). Let us
briefly review these issues.
It is easy to construct examples where q(D̂n)≤ α implies that π(D̂n)≈ α.
Consider the following two models for random variables Z = (Y,X1,X2):
Model 1
X1 ∼N(0,1),
Y = ψX1 +N(0,1),
X2 = ρX1 +N(0, τ
Model 2
X2 ∼N(0,1),
Y = ψX2 +N(0,1),
X1 = ρX2 +N(0, τ
Under models 1 and 2, the marginal distribution of Z is P1 =N(0,Σ1)
and P2 =N(0,Σ2), where
ψ2 + 1 ψ ρψ
ψ 1 ρ
ρψ ρ ρ2 + τ2
 , Σ2 =
ψ2 + 1 ρψ ψ
ρψ ρ2 + τ2 ρ
ψ ρ 1
Given any ε > 0, we can choose ρ sufficiently close to 1 and τ sufficiently close
to 0 such that Σ1 and Σ2 are as close as we like, and hence, d(P
1 , P
2 )< ε
where d is total variation distance. It follows that
P2(2 /∈ D̂)≥ P1(2 /∈ D̂)− ε≥ 1− α− ε.
Thus, if q ≤ α, then the power is less than α+ ε.
Dimensionality is less of an issue thanks to recent methods. Most methods,
including those in this paper, allow pn to grow exponentially. But all the
methods require some restrictions on the number s of nonzero βj ’s. In other
words, some sparsity assumption is required. In this paper we take s fixed
and allow pn to grow.
False negatives can occur during screening due to cancellations of cor-
relations. For example, the correlation between Y and X1 can be 0 even
when β1 is huge. This problem is called unfaithfulness in the causality liter-
ature [see Spirtes, Glymour and Scheines (2001) and Robins et al. (2003)].
False negatives during screening can lead to false positives during the second
stage.
6 L. WASSERMAN AND K. ROEDER
Let µ̂j denote the regression coefficient from regressing Y on Xj . Fix j ≤ s
and note that
µj ≡ E(µ̂j) = βj +
k 6=j
1≤k≤s
βkρkj,
where ρkj = corr(Xk,Xj). If
k 6=j
1≤k≤s
βkρkj ≈−βj,
then µj ≈ 0 no matter how large βj is. This problem can occur even when
n is large and p is small.
For example, suppose that β = (10,−10,0,0) and that ρ(Xi,Xj) = 0 ex-
cept that ρ(X1,X2) = ρ(X1,X3) = ρ(X2,X4) = 1− ε, where ε > 0 is small.
Then,
β = (10,−10,0,0), but µ≈ (0,0,10,−10).
Marginal regression is extremely susceptible to unfaithfulness. The lasso
and forward stepwise, less so. However, unobserved covariates can induce
unfaithfulness in all the methods.
3. Loss and cross-validation. Let Xλ = (X•j : j ∈ Ŝn(λ)) denote the de-
sign matrix corresponding to the covariates in Ŝn(λ) and let β̂(λ) be the
least-squares estimator for the regression restricted to Ŝn(λ), assuming the
estimator is well defined. Hence, β̂(λ) = (XTλXλ)
−1XTλ Y . More generally,
β̂M is the least-squares estimator for any subset of variables M . When con-
venient, we extend β̂(λ) to length p by setting β̂j(λ) = 0, for j /∈ Ŝn(λ).
3.1. Loss. Now we record some properties of the loss function. The first
part of the following lemma is essentially Lemma 3 of Meinshausen and Yu
(2008).
Lemma 3.1. Let M+m = {M ⊂ S : |M | ≤m,D⊂M}. Then,
L(β̂M )≤
4m log p
nφn(m)
→ 1.(6)
Let M−m = {M ⊂ S : |M | ≤m,D 6⊂M}. Then,
L(β̂M )≥ ψ2φn(m+ s)
→ 1.(7)
VARIABLE SELECTION 7
3.2. Cross-validation. Recall that the data have been split into groups
D1, D2 and D3 each of size n. Construct β̂(λ) from D1 and let
L̂(λ) =
Xi∈D2
(Yi −XTi β̂(λ))
We would like L̂(λ) to order the models the same way as the true loss
L(λ) [defined after (5)]. This requires that, asymptotically, L̂(λ)−L(λ)≈ δn,
where δn does not involve λ. The following bounds will be useful. Note that
L(λ) and L̂(λ) are both step functions that only change value when a variable
enters or leaves the model.
Theorem 3.2. Suppose that maxλ∈Λn |Ŝn(λ)| ≤ kn. Then, there exists
a sequence of random variables δn =OP (1) that do not depend on λ or X,
such that, with probability tending to 1,
|L(λ)− L̂(λ)− δn|=OP
n1−c2
4. Multi-stage methods. Themulti-stage methods use the following steps.
As mentioned earlier, we randomly split the data into three parts, D1, D2
and D3, which we take to be of equal size:
1. Stage I. Use D1 to find Ŝn(λ), for each λ.
2. Stage II. Use D2 to find λ̂ by cross-validation, and let Ŝn = Ŝn(λ̂).
3. Stage III. Use D3 to find the least-squares estimate β̂ for the model Ŝn.
D̂n = {j ∈ Ŝn : |Tj |> cn},
where Tj is the usual t-statistic, cn = zα/2m and m= |Ŝn|.
4.1. The lasso. The lasso estimator [Tibshirani (1996)] β̃(λ) minimizes
Mλ(λ) =
(Yi −XTi β)2 + λ
|βj |
and let Ŝn(λ) = {j : β̃j(λ) 6= 0}. Recall that β̂(λ) is the least-squares estima-
tor using the covariates in Ŝn(λ).
Let kn =A logn where A> 0 is a positive constant.
Theorem 4.1. Assume that (A1)–(A5) hold. Let Λn = {λ : |Ŝn(λ)| ≤
kn}. Then:
1. The true loss overfits: P(D ⊂ Ŝn(λ∗))→ 1 where λ∗ = argminλ∈Λn L(λ).
8 L. WASSERMAN AND K. ROEDER
2. Cross-validation also overfits: P(D⊂ Ŝn(λ̂))→ 1 where λ̂= argminλ∈Λn L̂(λ).
3. Type I error is controlled: lim supn→∞ P(D
c ∩ D̂n 6=∅) ≤ α.
If we let α= αn→ 0, then D̂n is consistent for variable selection.
Theorem 4.2. Assume that (A1)–(A5) hold. Let αn→ 0 and
∞. Then, the multi-stage lasso is consistent,
P(D̂n =D)→ 1.(10)
The next result follows directly. The proof is thus omitted.
Theorem 4.3. Assume that (A1)–(A5) hold. Let α be fixed. Then,
(D̂n, Ŝn) forms a confidence sandwich
lim inf
P(D̂n ⊂D ⊂ Ŝn)≥ 1−α.(11)
Remark 4.4. This confidence sandwich is expected to be conservative
in the sense that the coverage can be much larger than 1− α.
4.2. Stepwise regression. Let kn = A logn for some A > 0. The version
of stepwise regression we consider is as follows:
1. Initialize: Res = Y , λ= 0, Ŷ = 0 and Ŝn(λ) =∅.
2. Let λ← λ+1. Compute µ̂j = n−1〈Xj ,Res〉 for j = 1, . . . , p.
3. Let J = argmaxj |µ̂j |. Set Ŝn(λ) = {Ŝn(λ−1), J}. Set Ŷ =Xλβ̂(λ) where
β̂λ = (X
−1XTλ Y , and let Res = Y − Ŷ .
4. If λ= kn, stop. Otherwise, go to step 2.
For technical reasons, we assume that the final estimator xT β̂ is truncated
to be no larger than B. Note that λ is discrete and Λn = {0,1, . . . , kn}.
Theorem 4.5. With Ŝn(λ) defined as above, the statements of Theorems
4.1, 4.2 and 4.3 hold.
4.3. Marginal regression. This is probably the oldest, simplest and most
common method. It is quite popular in gene expression analysis. It used to
be regarded with some derision but has enjoyed a revival. A version appears
in a recent paper by Fan and Lv (2008). Let Ŝn(λ) = {j : |µ̂j | ≥ λ} where
µ̂j = n
−1〈Y,X•j〉.
Let µj = E(µ̂j), and let µ(j) denote the value of µ ordered by their absolute
values,
|µ(1)| ≥ |µ(2)| ≥ · · · .
VARIABLE SELECTION 9
Theorem 4.6. Let kn →∞ with kn = o(
n). Let Λn = {λ : |Ŝn(λ)| ≤
kn}. Assume that
|µj|> |µ(kn)|.(12)
Then, the statements of Theorems 4.1, 4.2 and 4.3 hold.
Assumption (12) limits the degree of unfaithfulness (small partial corre-
lations induced by cancellation of parameters). Large values of kn weaken
assumption (12), thus making the method more robust to unfaithfulness, but
at the expense of lower power. Fan and Lv (2008) make similar assumptions.
They assume that there is a C > 0 such that |µj| ≥ C|βj | for all j, which
rules out unfaithfulness. However, they do not explicitly relate the values of
µj for j ∈D to the values outside D as we have done. On the other hand,
they assume that Z =Σ−1/2X has a spherically symmetric distribution. Un-
der this assumption and their faithfulness assumption, they deduce that the
µj ’s outside D cannot strongly dominate the µj ’s within D. We prefer to
simply make this an explicit assumption without placing distributional as-
sumptions on X . At any rate, any method that uses marginal regressions as
a starting point must make some sort of faithfulness assumptions to succeed.
4.4. Modifications. Let us now discuss a few modifications of the basic
method. First, consider splitting the data only into two groups, D1 and D2.
Then do the following steps:
1. Stage I. Find Ŝn(λ) for λ ∈Λn, where |Ŝn(λ)| ≤ kn for each λ ∈ Λn using
2. Stage II. Find λ̂ by cross-validation, and let Ŝn = Ŝn(λ̂) using D2.
3. Stage III. Find the least-squares estimate β̂
using D2. Let D̂n = {j ∈
Ŝn : |Tj |> cn}, where Tj is the usual t-statistic.
Theorem 4.7. Choosing
log logn
2kn log(2pn)
controls asymptotic type I error.
The critical value in (13) is hopelessly large and it does not appear it can
be substantially reduced. We present this mainly to show the value of the
extra data-splitting step. It is tempting to use the same critical value as in
the tri-split case, namely cn = zα/2m where m = |Ŝn|, but we suspect this
will not work in general. However, it may work under extra conditions.
10 L. WASSERMAN AND K. ROEDER
5. Application. As an example, we illustrate an analysis based on part
of the osteoporotic fractures in men study [MrOS, Orwoll et al. (2005)]. A
sample of 860 men were measured at a large number of genes and outcome
measures. We consider only 296 SNPs which span 30 candidate genes for
bone mineral density. An aim of the study was to identify genes associated
with bone mineral density that could help in understanding the genetic basis
of osteoporosis in men. Initial analyses of this subset of the data revealed
no SNPs with a clear pattern of association with the phenotype; however,
three SNPs, numbered (67), (277) and (289), exhibited some association in
the screening of the data. To further explore the effacacy of the lasso screen
and clean procedure, we modified the phenotype to enhance this weak signal
and then reanalyzed the data to see if we could detect this planted signal.
We were interested in testing for main effects and pairwise interactions
in these data; however, including all interactions results in a model with
43,660 additional terms, which is not practical for this sample size. As a
compromise, we selected 2 SNPs per gene to model potential interaction
effects. This resulted in a model with a total of 2066 potential coefficients,
including 296 main effects and 1770 interaction terms. With this model, our
initial screen detected 10 terms, including the 3 enhanced signals, 2 other
main effects and 5 interactions. After cleaning, the final model detected the
3 enhanced signals and no other terms.
6. Simulations. To further explore the screen and clean procedures, we
conducted simulation experiments with four models. For each model Yi =
XTi β + εi where the measurement errors, εi and ε
ij , are i.i.d. normal(0,1)
and the covariates Xij ’s are normal(0,1) (except for model D). Models differ
in how Yi is linked to Xi and the dependence structure of the Xi’s. Models
A, B and C explore scenarios with moderate and large p, while model D
focuses on confounding and unfaithfullness, as follows:
(A) Null model: β = (0, . . . ,0) and the Xij ’s are i.i.d.
(B) Triangle model: βj = δ(10− j), j = 1, . . . ,10, βj = 0, j > 10 and Xij ’s are
i.i.d.
(C) Correlated Triangle model: as B, but withXij(+1) = ρXij+(1−ρ2)1/2ε∗ij ,
for j > 1, and ρ= 0.5.
(D) Unfaithful model: Yi = β1Xi1 + β2Xi2 + εi, for β1 = −β2 = 10, where
the Xij ’s are i.i.d. for j = {1,5,6,7,8,9,10}, but Xi2 = ρXi1 + τε∗i2,
Xi3 = ρXi1 + τε
i10, and Xi4 = ρXi2 + τε
i11, for τ = 0.01 and ρ= 0.95.
We used a maximum model size of kn = n
1/2 which technically goes be-
yond the theory but works well in practice. Prior to analysis, the covariates
are scaled so that each has mean 0 and variance 1. The tests were initially
performed using a third of the data for each of the 3 stages of the proce-
dure (Table 1, top half, 3 splits). For models A, B and C, each approach
VARIABLE SELECTION 11
Table 1
Size and power of screen and clean procedures using lasso, stepwise and marginal
regression for the screening step. For all procedures α= 0.05. For p= 100, δ = 0.5 and
for p= 1000, δ = 1.5. Reported power is πav. The top 8 rows of simulations were
conducted using three stages as described in Section 4, with a third of the data used for
each stage. The bottom 8 rows of simulations were conducted splitting the data in half,
using the first portion with leave-one-out cross validation for stages 1 and 2 and the
second portion for cleaning
Size Power
Splits n p Model Lasso Step Marg Lasso Step Marg
2 100 100 A 0.005 0.001 0.004 0 0 0
2 100 100 B 0.01 0.02 0.03 0.62 0.62 0.31
2 100 100 C 0.001 0.01 0.01 0.77 0.57 0.21
2 100 10 D 0.291 0.283 0.143 0.08 0.08 0.04
2 100 1000 A 0.001 0.002 0.010 0 0 0
2 100 1000 B 0.002 0.020 0.010 0.17 0.09 0.11
2 100 1000 C 0.02 0.14 0.01 0.27 0.15 0.11
2 1000 10 D 0.291 0.283 0.143 0.08 0.08 0.04
3 100 100 A 0.040 0.050 0.030 0 0 0
3 100 100 B 0.02 0.01 0.02 0.91 0.90 0.56
3 100 100 C 0.03 0.04 0.03 0.91 0.88 0.41
3 100 10 D 0.382 0.343 0.183 0.16 0.18 0.09
3 100 1000 A 0.035 0.045 0.040 0 0 0
3 100 1000 B 0.045 0.020 0.035 0.57 0.66 0.29
3 100 1000 C 0.06 0.070 0.020 0.74 0.65 0.19
3 1000 10 D 0.481 0.486 0.187 0.17 0.17 0.13
has type I error less than α, except the stepwise procedure which has trou-
ble with model C when n = p = 100. We also calculated the false positive
rate and found it to be very low (about 10−4 when p= 100 and 10−5 when
p= 1000) indicating that even when a type I error occurs, only a very small
number of terms are included erroneously. The lasso screening procedure
exhibited a slight power advantage over the stepwise procedure. Both meth-
ods dominated the marginal approach. The Markov dependence structure
in model C clearly challenged the marginal approach. For model D, none of
the approaches controlled the type I error.
To determine the sensitivity of the approach to using distinct data for each
stage of the analysis, simulations were conducted screening on the first half
of the data and cleaning on the second half (2 splits). The tuning parameter
was selected using leave-one-out cross validation (Table 1, bottom half). As
expected, this approach lead to a dramatic increase in the power of all the
procedures. More surprising is the fact that the type I error was near α or
12 L. WASSERMAN AND K. ROEDER
below for models A, B and C. Clearly this approach has advantages over
data splitting and merits further investigation.
A natural competitor to screen and clean procedure is a two-stage adap-
tive lasso [Zou (2006)]. In our implementation, we split the data and used
one half for each stage of the analysis. At stage one, leave-one-out cross val-
idation lasso screens the data. In stage two, the adaptive lasso, with weights
wj = |β̂j |−1, cleans the data. The tuning parameter for the lasso was again
chosen using leave-one-out cross validation. Table 2 provides the size, power
and false positive rate (FPR) for this procedure. Naturally, the adaptive
lasso does not control the size of the test, but the FPR is small. The power
of the test is greater than we found for our lasso screen and clean procedure,
but this extra power comes at the cost of a much higher type I error rate.
7. Proofs. Recall that if A is a square matrix, then φ(A) and Φ(A)
denote the smallest and largest eigenvalues of A. Throughout the proofs we
make use of the following fact. If v is a vector and A is a square matrix,
φ(A)‖v‖2 ≤ vTAv ≤Φ(A)‖v‖2.(14)
We use the following standard tail bound: if Z ∼N(0,1), then P(|Z|> t)≤
t−1e−t
2/2. We will also use the following results about the lasso from Mein-
shausen and Yu (2008). Their results are stated and proved for fixed X but,
under the conditions (A1)–(A5), it is easy to see that their conditions hold
with probability tending to one and so their results hold for random X as
well.
Theorem 7.1 [Meinshausen and Yu (2008)]. Let β̃(λ) be the lasso esti-
mator.
Table 2
Size, power and false positive rate (FPR) of two-stage adaptive lasso procedure
n p Model Size Power FPR
100 100 A 0.93 0 0.032
100 100 B 0.84 0.97 0.034
100 100 C 0.81 0.96 0.031
100 10 D 0.67 0.21 0.114
100 1000 A 0.96 0 0.004
100 1000 B 0.89 0.65 0.004
100 1000 C 0.76 0.77 0.002
1000 10 D 0.73 0.24 0.013
VARIABLE SELECTION 13
1. The squared error satisfies
‖β̃(λ)− β‖22 ≤
cm log pn
nφ2n(m)
→ 1,(15)
where m= |Ŝn(λ)| and c > 0 is a constant.
2. The size of Ŝn(λ) satisfies
|Ŝn(λ)| ≤
τ2Cn2
→ 1,(16)
where τ2 = E(Y 2i ).
Proof of Lemma 3.1. Let D ⊂M and φ= φ(n−1XTMXM ). Then,
L(β̂M ) =
εTXM (X
MXM )
−1XTMε≤
‖XTMε‖2 =
Z2j ,
where Zj = n
−1/2XT•jε. Conditional on X , Zi ∼N(0, a2j ) where a2j = n−1 ×∑n
ij . Let A
n =max1≤j≤pn a
j . By Hoeffding’s inequality, (A2) and (A5),
P(En)→ 1 where En = {An ≤
2}. So
1≤j≤pn
|Zj |>
4 log pn
1≤j≤pn
|Zj |>
4 log pn,En
1≤j≤pn
|Zj |>
4 log pn,E
1≤j≤pn
|Zj |>
4 log pn,En
+ P(Ecn)
An max
1≤j≤pn
|Zj |
4 log pn,En
+ o(1)
1≤j≤pn
|Zj |
2 log pn
+ o(1)
1≤j≤pn
|Zj |
2 log pn
)∣∣∣X
+ o(1)
2 log pn
+ o(1) = o(1).
j∈M Z
j ≤mmax1≤j≤pn Z2j and (6) follows.
Now we lower bound L(β̂M ). Let M be such that D 6⊂ M . Let A =
{j : β̂M (j) 6= 0} ∪D. Then, |A| ≤m + s. Therefore, with probability tend-
ing to 1,
L(β̂M ) =
(β̂M − β)TXTX(β̂M − β) =
(β̂M − β)TXTAXA(β̂M − β)
14 L. WASSERMAN AND K. ROEDER
≥ φn(m+ s)‖β̂M − β‖2 = φn(m+ s)
(β̂M (j)− β(j))2
≥ φn(m+ s)
j∈D∩Mc
(0− β(j))2 ≥ φn(m+ s)ψ2.
Proof of Theorem 3.2. Let Ỹ denote the responses, and X̃ the design
matrix, for the second half of the data. Then, Ỹ = X̃β + ε̃. Now
L(λ) =
(β̂(λ)− β)TXTX(β̂(λ)− β) = (β̂(λ)− β)T Σ̂n(β̂(λ)− β)
L̂(λ) = n−1‖Ỹ −X̃β̂(λ)‖2 = (β̂(λ)−β)T Σ̃n(β̂(λ)−β)+δn+
〈ε̃, X̃(β̂(λ)−β)〉,
where δn = ‖ε̃‖2/n, and Σ̂n = n−11 XTX and Σ̃n = n−1X̃T X̃ . By Hoeffding’s
inequality
P(|Σ̂n(j, k)− Σ̃n(j, k)|> ε)≤ e−ncε
for some c > 0, and so
|Σ̂n(j, k)− Σ̃n(j, k)|> ε
≤ p2ne−ncε
Choose εn = 4/(cn
1−c2). It follows that
|Σ̂n(j, k)− Σ̃n(j, k)|>
cn1−c2
≤ e−2nc2 → 0.
Note that
|{j : β̂j(λ) 6= 0} ∪ {j :βj 6= 0}| ≤ kn + s.
Hence, with probability tending to 1
|L(λ)− L̂(λ)− δn| ≤
cn1−c2
‖β̂(λ)− β‖21 + 2ξn(λ)
for all λ ∈Λn, where
ξn(λ) =
ε̃iµi(λ)
and µi(λ) = X̃
i (β̂(λ)−β). Now ‖β̂(λ)−β‖21 =OP ((kn+s)2) since ‖β̂(λ)‖2 =
OP (kn/φ(kn)). Thus, ‖β̂(λ)−β‖1 ≤C(kn+s) with probability tending to 1,
for some C > 0. Also, |µi(λ)| ≤B‖β̂(λ)−β‖1 ≤BC(kn+ s) with probability
tending to 1. Let W ∼N(0,1). Conditional on D1,
|ξn(λ)|
µ2i (λ)|W | ≤
BC(kn + s)|W |,
VARIABLE SELECTION 15
so supλ∈Λn |ξn(λ)|=OP (kn/
n). �
Proof of Theorem 4.1. 1. Let λn = τn
C/kn,M = Ŝn(λn) and m=
|M |. Then, P(m≤ kn)→ 1 due to (16). Hence, P(λn ∈ Λn)→ 1. From (15),
‖β̃(λn)− β‖22 ≤O
kn log pn
= oP (1).
Hence, ‖β̃(λn)− β‖2∞ = oP (1). So, for each j ∈D,
|β̃j(λn)| ≥ |βj | − |β̃j(λn)− βj | ≥ ψ+ oP (1)
and hence, P(minj∈D |β̃j(λn)|> 0)→ 1. Therefore, Γn = {λ ∈Λn :D ⊂ Ŝn(λ)}
is nonempty. By Lemma 3.1,
L(λn)≤ cm log pn/(nφ(m)) =OP (kn log pn/n).(17)
On the other hand, from Lemma 3.1,
λ∈Λn∩Γcn
L(β̂λ)>ψ
2φ(kn)
→ 1.(18)
Now, nφn(kn)/(kn log pn)→∞, and so (17) and (18) imply that
λ∈Λn∩Γcn
L(β̂λ)>L(λn)
Thus, if λ∗ denotes the minimizer of L(λ) over Λn, we conclude that P(λ∗ ∈
Γn)→ 1, and hence, P(D ⊂ Ŝn(λ∗))→ 1.
2. This follows from part 1 and Theorem 3.2.
3. Let A= Ŝn ∩Dc. We want to show that
|Tj|> cn
≤ α+ o(1).
|Tj |> cn
|Tj |> cn,D ⊂ Ŝn
|Tj |> cn,D 6⊂ Ŝn
|Tj |> cn,D ⊂ Ŝn
+ P(D 6⊂ Ŝn)
|Tj |> cn,D ⊂ Ŝn
+ o(1).
Conditional on (D1,D2), β̂A is normally distributed with mean 0 and vari-
ance matrix σ2(XTAXA)
−1 when D⊂ Ŝn. Recall that
Tj(M) =
eTj (X
MXM )
−1XTMY
eTj (X
MXM )
β̂M,j
16 L. WASSERMAN AND K. ROEDER
whereM = Ŝn, s
j = σ̂
2eTj (X
MXM )
−1ej and ej = (0, . . . ,0,1,0, . . . ,0)
T , where
the 1 is in the jth coordinate. When D ⊂ Ŝn, each Tj , for j ∈ A, has a t-
distribution with n−m degrees of freedom wherem= |Ŝn|. Also, cn/tα/2m→
1 where tu denotes the upper tail critical value for the t-distribution. Hence,
|Tj |> cn,D ⊂ Ŝn|D1,D2
|Tj |> tα/2m,D ⊂ Ŝn|D1,D2
≤ α+ an,
where an = o(1), since |A| ≤m. It follows that
|Tj |> cn,D ⊂ Ŝn
≤ α+ o(1).
Proof of Theorem 4.2. From Theorem 4.1, P(D̂n ∩Dc 6= ∅) ≤ αn
and so P(D̂n ∩Dc 6=∅)→ 0. Hence, P(D̂n ⊂D)→ 1. It remains to be shown
P(D⊂ D̂n)→ 1.(19)
The test statistic for testing βj = 0 when Ŝn =M is
Tj(M) =
eTj (X
MXM )
−1XTMY
eTj (X
MXM )
For simplicity in the proof, let us take σ̂ = σ, the extension to unknown σ
being straightforward. Let j ∈D,M= {M : |M | ≤ kn,D⊂M}. Then,
P(j /∈ D̂n) = P(j /∈ D̂n,D ⊂ Ŝn) + P(j /∈ D̂n,D 6⊂ Ŝn)
≤ P(j /∈ D̂n,D ⊂ Ŝn) + P(D 6⊂ Ŝn)
= P(j /∈ D̂n,D ⊂ Ŝn) + o(1)
P(j /∈ D̂n, Ŝn =M) + o(1)
P(|Tj(M)|< cn, Ŝn =M) + o(1)
P(|Tj(M)|< cn) + o(1).
Conditional on D1 ∪D2, for each M ∈M, Tj(M) = (βj/sj) +Z, where Z ∼
N(0,1). Without loss of generality, assume that βj > 0. Hence,
P(|Tj(M)|< cn|D1 ∪D2) = P
−cn −
<Z < cn −
VARIABLE SELECTION 17
Fix a small ε > 0. Note that s2j ≤ σ2/(nκ). It follows that, for all large n,
cn − βj/sj <−ε
n. So,
P(|Tj(M)|< cn|D1 ∪D2)≤ P(Z <−ε
n)≤ e−nε2/2.
The number of models inM is
pn − s
j − s
pn − s
kn − s
(pn − s)e
kn − s
)kn−s
≤ knpknn ,
where we used the inequality
P(|Tj(M)|< cn|D1 ∪D2)≤ knpknn e−nε
2 → 0
by (A2). We have thus shown that P(j /∈ D̂n)→ 0, for each j ∈D. Since |D|
is finite, it follows that P(j /∈ D̂n for some j ∈D)→ 0 and hence (19). �
Proof of Theorem 4.5. A simple modification of Theorem 3.1 of
Barron et al. (2008) shows that
L(kn) =
‖Ŷkn −Xβ‖2 = oP (1).
[The modification is needed because Barron et al. (2008) require Y to be
bounded while we have assumed that Y is normal. By a truncation argument,
we can still derive the bound on L(kn).] So
‖β̂kn − β‖2 ≤
L(kn)
φn(kn + s)
L(kn)
= oP (1).
Hence, for any ε > 0, with probability tending to 1, ‖β̂(kn)−β‖2 < ε so that
|β̂j | > ψ/2 > 0, for all j ∈D. Thus, P(D ⊂ Ŝn(kn))→ 1. The remainder of
the proof of part 1 is the same as in Theorem 4.1. Part 2 follows from the
previous result together with Theorem 3.2. The proof of part 3 is the same
as for Theorem 4.1. �
Proof of Theorem 4.6. Note that µ̂j − µj = n−1
i=1Xijεi. Hence,
µ̂j − µj ∼N(0,1/n). So, for any δ > 0,
|µ̂j − µj |> δ
P(|µ̂j − µj|> δ)
2/2 ≤ c1e
2/2→ 0.
18 L. WASSERMAN AND K. ROEDER
By (12), conclude that D ⊂ Ŝn(λ) when λ = µ̂(kn). The remainder of the
proof is the same as the proof of Theorem 4.5. �
Proof of Theorem 4.7. Let A= Ŝn ∩Dc. We want to show that
|Tj|> cn
≤ α+ o(1).
For fixed A, β̂A is normal with mean 0 but this is not true for random A.
Instead we need to bound Tj . Recall that
Tj(M) =
eTj (X
MXM )
−1XTMY
eTj (X
MXM )
β̂M,j
whereM = Ŝn, s
j = σ̂
2eTj (X
MXM )
−1ej and ej = (0, . . . ,0,1,0, . . . ,0)
T where
the 1 is in the jth coordinate. The probabilities that follow are conditional
on D1 but this is supressed for notational convenience. First, write
|Tj |> cn
|Tj |> cn,D ⊂ Ŝn
|Tj |> cn,D 6⊂ Ŝn
|Tj |> cn,D ⊂ Ŝn
+ P(D 6⊂ Ŝn)
|Tj |> cn,D ⊂ Ŝn
+ o(1).
When D ⊂ Ŝn,
)−1 1
where Q
= ((1/n)XT
)−1, γ
= n−1XT
ε, and β
(j) = 0, for j ∈A.
Now, s2j ≥ σ̂2/(nC) so that
|Tj | ≤
nC|β̂
Ŝn,j
n log logn|β̂
Ŝn,j
for j ∈ Ŝn. Therefore,
|Tj |> cn,D ⊂ Ŝn
Ŝn,j
σ̂cn√
,D ⊂ Ŝn
Let γ = n−1XT ε. Then,
‖β̂A‖2 ≤ γT
knmax1≤j≤pn γ
VARIABLE SELECTION 19
It follows that
Ŝn,j
knmax1≤j≤pn |γj |
kn log logn max
1≤j≤pn
|γj |,
since κ > 0. So,
Ŝn,j
|> σ̂cn√
n log logn
,D ⊂ Ŝn
1≤j≤pn
|γj |>
log logn
Note that γj ∼N(0, σ2/n), and hence,
|γj |
2σ2 log(2pn)
There exists εn → 0 such that P(Bn)→ 1 where Bn = {(1 − εn) ≤ σ̂/σ ≤
(1 + ε)}. So,
1≤j≤pn
|γj|>
log logn
1≤j≤pn
|γj |>
σcn(1− εn)
log logn
σ(1− εn)cn
log logn
|γj |
≤ α+ o(1). �
8. Discussion. The multi-stage method presented in this paper success-
fully controls type I error while giving reasonable power. The lasso and step-
wise have similar performance. Although theoretical results assume indepen-
dent data for each of the three stages, simulations suggest that leave-one-out
cross-validation leads to valid type I error rates and greater power. Screening
the data in one phase of the experiment and cleaning in a followup phase
leads to an efficient experimental design. Certainly this approach deserves
further theoretical investigation. In particular, the question of optimality is
an open question.
The literature on high-dimensional variable selection is growing quickly.
The most important deficiency in much of this work, including this paper,
is the assumption that the model Y = XTβ + ε is correct. In reality, the
model is at best an approximation. It is possible to study linear procedures
when the linear model is not assumed to hold as in Greenshtein and Ritov
(2004). We discuss this point in the Appendix. Nevertheless, it seems useful
to study the problem under the assumption of linearity to gain insight into
these methods. Future work should be directed at exploring the robustness
of the results when the model is wrong.
Other possible extensions include: dropping the normality of the errors,
permitting nonconstant variance, investigating the optimal sample sizes for
each stage and considering other screening methods besides cross-validation.
20 L. WASSERMAN AND K. ROEDER
Finally, let us note that the example involving unfaithfulness, that is,
cancellations of parameters to make the marginal correlation much different
than the regression coefficient, pose a challenge for all the methods and
deserve more attention even in cases of small p.
APPENDIX: PREDICTION
Realistically, there is little reason to believe that the linear model is cor-
rect. Even if we drop the assumption that the linear model is correct, sparse
methods like the lasso can still have good properties as shown in Greenshtein
and Ritov (2004). In particular, they showed that the lasso satisfies a risk
consistency property. In this appendix we show that this property continues
to hold if λ is chosen by cross-validation.
The lasso estimator is the minimizer of
i=1(Yi −XTi β)2 + λ‖β‖1. This
is equivalent to minimizing
i=1(Yi−XTi β)2 subject to ‖β‖1 ≤Ω, for some
Ω. (More precisely, the set of estimators as λ varies is the same as the set of
estimators as Ω varies.) We use this second version throughout this section.
The predictive risk of a linear predictor ℓ(x) = xTβ is R(β) = E(Y −ℓ(x))2
where (X,Y ) denotes a new observation. Let γ = γ(β) = (−1, β1, . . . , βp)T
and let Γ = E(ZZT ) where Z = (Y,X1, . . . ,Xp). Then, we can write R(β) =
γTΓγ. The lasso estimator can now be written as β̂(Ωn) = argminβ∈B(Ωn) R̂(β)
where R̂(β) = γT Γ̂γ and Γ̂ = n−1
i=1ZiZ
Define
β∗ = argmin
β∈B(Ωn)
R(β),
where
B(Ωn) = {β :‖β‖1 ≤Ωn}.
Thus, ℓ∗(x) = x
Tβ∗ is the best linear predictor in the set B(Ωn). The best
linear predictor is well defined even though E(Y |X) is no longer assumed to
be linear. Greenshtein and Ritov (2004) call an estimator β̂n persistent, or
predictive risk consistent, if
R(β̂n)−R(β∗)
as n→∞.
The assumptions we make in this section are:
(B1) pn ≤ en
, for some 0≤ ξ < 1;
(B2) the elements of Γ̂ satisfy an exponential inequality
P(|Γ̂jk − Γjk|> ε)≤ c3e−nc4ε
for some c3, c4 > 0;
VARIABLE SELECTION 21
(B3) there exists B0 <∞ such that, for all n, maxj,kE(|ZjZk|)≤B0.
Condition (A2) can easily be deduced from more primitive assumptions
as in Greenshtein and Ritov (2004), but for simplicity we take (A2) as an
assumption. Let us review one of the results in Greenshtein and Ritov (2004).
For the moment, replace (A1) with the assumption that pn ≤ nb, for some
b. Under these conditions, it follows that
∆n ≡max
|Γ̂jk − Γjk|=OP
Hence,
β∈B(Ωn)
|R(β)− R̂(β)|= sup
β∈B(Ωn)
|γT (Γ− Γ̂)γ|
≤∆n sup
β∈B(Ωn)
‖γ‖21 =Ω2nOP
The latter term is oP (1) as long as Ωn = o((n/ logn)
1/4). Thus, we have the
following.
Theorem A.1 [Greenshtein and Ritov (2004)]. If Ωn = o((n/ logn)
1/4),
then the lasso estimator is persistent.
For future reference, let us state a slightly different version of their result
that we will need. We omit the proof.
Theorem A.2. Let γ > 0 be such that ξ+γ < 1. Let Ωn =O(n
(1−ξ−γ)/4).
Then, under (B1) and (B2),
β∈B(Ωn)
|R̂(β)−R(β)|> 1
=O(e−cn
)(20)
for some c > 0.
The estimator β̂(Ωn) lies on the boundary of the ball B(Ωn) and is very
sensitive to the exact choice of Ωn. A potential improvement—and something
that reflects actual practice—is to compute the set of lasso estimators β̂(ℓ),
for 0≤ ℓ≤ Ωn and then select from that set based on cross validation. We
now confirm that the resulting estimator preserves persistence. As before,
we split the data into D1 and D2. Construct the lasso estimators {β̂(ℓ) : 0≤
ℓ≤Ωn}. Choose ℓ̂ by cross validation using D2. Let β̂ = β̂(ℓ̂).
22 L. WASSERMAN AND K. ROEDER
Theorem A.3. Let γ > 0 be such that ξ + γ < 1. Under (A1), (A2)
and (A3), if Ωn =O(n
(1−ξ−γ)/4), then the cross-validated lasso estimator β̂
is persistent. Moreover,
R(β̂)− inf
0≤ℓ≤Ωn
R(β̂(ℓ))
P→ 0.(21)
Proof. Let β∗(ℓ) = argminβ∈B(ℓ)R(β). Define h(ℓ) = R(β∗(ℓ)), g(ℓ) =
R(β̂(ℓ)) and c(ℓ) = L̂(β̂(ℓ)). Note that, for any vector b, we can write R(b) =
τ2 + bTΣb− 2bTρ where ρ= (E(Y X1), . . . ,E(Y Xp))T .
Clearly, h is monotone nonincreasing on [0,Ωn]. We claim that |h(ℓ +
δ) − h(ℓ)| ≤ cΩnδ where c depends only on Γ. To see this, let u = β∗(ℓ),
v = β∗(ℓ+ δ) and a= ℓβ∗(ℓ+ δ)/(ℓ+ δ) so that a ∈B(ℓ). Then,
h(ℓ+ δ)≤ h(ℓ)
=R(u)≤R(a)
=R(v) +R(a)−R(v)
= h(ℓ+ δ) +
ρT v− δ(2ℓ+ δ)
(ℓ+ δ)2
≤ h(ℓ+ δ) + 2δC + δ(2Ωn + δ)C,
where C =maxj,k |Γj,k|=O(1).
Next, we claim that g(ℓ) is Lipschitz on [0,Ωn] with probability tending to
1. Let β̂(ℓ) = argminβ∈B(ℓ) R̂(β) denote the lasso estimator and set û= β̂(ℓ)
and v̂ = β̂(ℓ+δ). Let εn = n
−γ/4. From (20), the following chain of equations
hold except on a set of exponentially small probability:
g(ℓ+ δ) =R(v̂)≤ R̂(v̂) + εn ≤ R̂(v) + εn
≤R(v) + 2εn = h(ℓ+ δ) + 2εn
≤ h(ℓ) + cΩnδ +2εn =R(u) + cΩnδ+ 2εn
≤R(û) + cΩnδ+2εn = g(ℓ) + cΩnδ +2εn.
A similar argument can be applied in the other direction. Conclude that
|g(ℓ+ δ)− g(ℓ)| ≤ cΩnδ +2εn
except on a set of small probability.
Now let A= {0, δ,2δ, . . . ,mδ} where m is the smallest integer such that
mδ ≥ Ωn. Thus, m∼ Ωn/δn. Choose δ = δn = n−3(1−ξ−γ)/8. Then, Ωnδn→
0 and Ωn/δn ≤ n3(1−ξ−γ)/4. Using the same argument as in the proof of
Theorem 3.2,
|L̂(β̂(ℓ))−R(β̂(ℓ))|= σn,
VARIABLE SELECTION 23
where σn = oP (1). Then,
R(β∗(Ωn))≤R(β̂)≤ L̂(β̂(ℓ̂)) + σn
≤ L̂(mδn) + σn ≤ g(mδn) + 2σn ≤ g(Ωn) + 2σn + cΩnδn
≤ h(Ωn) + 2σn + εn + cΩnδn
=R(β∗(Ωn)) + 2σn + εn + cΩnδn
and persistence follows. To show the second result, let β̃ = argmin0≤ℓ≤Ωn g(ℓ)
and β = argminℓ∈A g(ℓ). Then,
R(β̃)≤ L̂(β̃) + σn ≤ L̂(β) + σn
≤R(β) + 2σn ≤R(β̃) + 2σn + cδnΩn
and the claim follows. �
Acknowledgments. The authors are grateful for the use of a portion of
the sample from the Osteoporotic Fractures in Men (MrOS) Study to il-
lustrate their methodology. MrOs is supported by the National Institute of
Arthritis and Musculoskeletal and Skin Diseases (NIAMS), the National In-
stitute on Aging (NIA) and the National Cancer Institute (NCI) through
Grants U01 AR45580, U01 AR45614, U01 AR45632, U01 AR45647, U01
AR45654, U01 AR45583, U01 AG18197 and M01 RR000334. Genetic analy-
ses in MrOS were supported by R01-AR051124. We also thank two referees
and an AE for helpful suggestions.
REFERENCES
Barron, A., Cohen, A., Dahmen, W. and DeVore, R. (2008). Approximation and
learning by greedy algorithms. Ann. Statist. 36 64–94. MR2387964
Bühlmann, P. (2006). Boosting for high-dimensional linear models. Ann. Statist. 34 559–
583. MR2281878
Candes, E. and Tao, T. (2007). The Dantzig selector: Statistical estimation when p is
much larger than n. Ann. Statist. 35 2313–2351. MR2382644
Donoho, D. (2006). For most large underdetermined systems of linear equations, the
minimal l1-norm near-solution approximates the sparsest near-solution. Comm. Pure
Appl. Math. 59 797–829. MR2217606
Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004). Least angle regres-
sion. Ann. Statist. 32 407–499. MR2060166
Fan, J. and Lv, J. (2008). Sure independence screening for ultra-high-dimensional feature
space. J. Roy. Statist. Soc. Ser. B 70 849–911.
Greenshtein, E. and Ritov, Y. (2004). Persistence in high-dimensional linear predictor
selection and the virtue of overparametrization. Bernoulli 10 971–988. MR2108039
Meinshausen, N. (2007). Relaxed lasso. Comput. Statist. Data Anal. 52 374–393.
MR2409990
Meinshausen, N. and Bühlmann, P. (2006). High dimensional graphs and variable se-
lection with the lasso. Ann. Statist. 34 1436–1462. MR2278363
http://www.ams.org/mathscinet-getitem?mr=2387964
http://www.ams.org/mathscinet-getitem?mr=2281878
http://www.ams.org/mathscinet-getitem?mr=2382644
http://www.ams.org/mathscinet-getitem?mr=2217606
http://www.ams.org/mathscinet-getitem?mr=2060166
http://www.ams.org/mathscinet-getitem?mr=2108039
http://www.ams.org/mathscinet-getitem?mr=2409990
http://www.ams.org/mathscinet-getitem?mr=2278363
24 L. WASSERMAN AND K. ROEDER
Meinshausen, N. and Yu, B. (2009). Lasso-type recovery of sparse representations of
high-dimensional data. Ann. Statist. To appear.
Orwoll, E., Blank, J. B., Barrett-Connor, E., Cauley, J., Cummings, S., En-
srud, K., Lewis, C., Cawthon, P. M., Marcus, R., Marshall, L. M., McGowan,
J., Phipps, K., Sherman, S., Stefanick, M. L. and Stone, K. (2005). Design and
baseline characteristics of the osteoporotic fractures in men (MrOS) study—a large ob-
servational study of the determinants of fracture in older men. Contemp. Clin. Trials.
26 569–585.
Robins, J., Scheines, R., Spirtes, P. and Wasserman, L. (2003). Uniform consistency
in causal inference. Biometrika 90 491–515. MR2006831
Spirtes, P., Glymour, C. and Scheines, R. (2001). Causation, Prediction, and Search.
MIT Press. MR1815675
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Roy. Statist.
Soc. Ser. B 58 267–288. MR1379242
Tropp, J. A. (2004). Greed is good: Algorithmic results for sparse approximation. IEEE
Trans. Inform. Theory 50 2231–2242. MR2097044
Tropp, J. A. (2006). Just relax: Convex programming methods for identifying sparse
signals in noise. IEEE Trans. Inform. Theory 52 1030–1051. MR2238069
Wainwright, M. (2006). Sharp thresholds for high-dimensional and noisy recovery of
sparsity. Available at arxiv.org/math.ST/0605740.
Wellcome Trust (2007). Genome-wide association study of 14,000 cases of seven com-
mon diseases and 3000 shared controls. Nature 447 661–678.
Zhao, P. and Yu, B. (2006). On model selection consistency of lasso. J. Mach. Learn.
Res. 7 2541–2563. MR2274449
Zhang, C. H. and Huang, J. (2006). Model selection consistency of the lasso in high-
dimensional linear regression. Ann. Statist. 36 1567–1594. MR2435448
Zou, H. (2006). The adaptive lasso and its oracle properties. J. Amer. Statist. Assoc. 101
1418–1429. MR2279469
Department of Statistics
Carnegie Mellon University
Pittsburgh
E-mail: larry@stat.cmu.edu
roeder@stat.cmu.edu
http://www.ams.org/mathscinet-getitem?mr=2006831
http://www.ams.org/mathscinet-getitem?mr=1815675
http://www.ams.org/mathscinet-getitem?mr=1379242
http://www.ams.org/mathscinet-getitem?mr=2097044
http://www.ams.org/mathscinet-getitem?mr=2238069
http://arxiv.org/math.ST/0605740
http://www.ams.org/mathscinet-getitem?mr=2274449
http://www.ams.org/mathscinet-getitem?mr=2435448
http://www.ams.org/mathscinet-getitem?mr=2279469
mailto:larry@stat.cmu.edu
mailto:roeder@stat.cmu.edu
	Introduction
	Error control
	Loss and cross-validation
	Loss
	Cross-validation
	Multi-stage methods
	The lasso
	Stepwise regression
	Marginal regression
	Modifications
	Application
	Simulations
	Proofs
	Discussion
	Appendix: Prediction
	Acknowledgments
	References
	Author's addresses
