Marginal false discovery rates for penalized regression models
Patrick Breheny
Department of Biostatistics
University of Iowa
April 10, 2017

arXiv:1607.05636v2 [math.ST] 7 Apr 2017

Abstract

1

Penalized regression methods are an attractive tool for high-dimensional data
analysis, but their widespread adoption has been hampered by the difficulty
of applying inferential tools. In particular, the question “How reliable is the
selection of those features?” has proved difficult to address. In part, this difficulty arises from defining false discoveries in the classical, fully conditional sense,
which is possible in low dimensions but does not scale well to high-dimensional
settings. Here, we consider the analysis of marginal false discovery rates for penalized regression methods. Restricting attention to the marginal FDR permits
straightforward estimation of the number of selections that would likely have
occurred by chance alone, and therefore provides a useful summary of selection
reliability. Theoretical analysis and simulation studies demonstrate that this
approach is quite accurate when the correlation among predictors is mild, and
only slightly conservative when the correlation is stronger. Finally, the practical utility of the proposed method and its considerable advantages over other
approaches are illustrated using gene expression data from The Cancer Genome
Atlas and GWAS data from the Myocardial Applied Genomics Network.

Introduction

Penalized regression is an attractive methodology for dealing with high-dimensional data where classical likelihood
approaches to modeling break down. However, its widespread adoption has been hindered by a lack of inferential
tools. In particular, penalized regression is very useful for variable selection, but the question “How reliable are
those selections?” has proved difficult to address. As I will argue in this paper, this difficulty is partially due to
the complexity of defining a “false discovery” in the penalized regression setting.
In this paper, I will focus mainly on the lasso for linear regression models (Tibshirani, 1996), although the idea
is very general and can be extended to a variety of other regression models and penalty functions. In particular,
suppose we use the lasso to select 10 variables from a pool of potentially important features. This paper addresses
the question, “How many of those selections would likely have occurred by chance alone?”
There has been a fair amount of recent work on the idea of hypothesis testing in the high-dimensional penalized
regression setting. A comprehensive review is beyond the scope of this paper, but it is worth introducing two
approaches to which I will compare the method proposed in this paper. One approach is to split the sample
into two parts, using the first part for variable selection and the second part for hypothesis testing. This idea
was first introduced in Wasserman and Roeder (2009), who studied the problem using a single split of the data
set. Meinshausen et al. (2009) extended this approach by considering multiple random splits and combining the
results. Dezeure et al. (2015) provides a comprehensive review of this approach and the details involved, along
with procedures for limiting the overall false discovery rate through this form of testing.
An alternative approach is to test the significance of adding a variable along the solution path as the degree
of penalization is relaxed, conditional on the other variables already included in the model. Several tests fall into
1

this category, including the covariance test (Lockhart et al., 2014), the spacing test, and an exact test (Tibshirani
et al., 2016, which also introduces the spacing test). The details of the tests differ, but all involve modifying
classical tests for the significance of an added variable by conditioning on the fact that the added variable was
not prespecified, but rather selected from a pool of potential variables. G’Sell et al. (2016) provide an important
extension to this work by deriving a rule, forwardStop, for selecting the stopping point along this sequence of
sequential tests so as to preserve a specified false discovery rate.
False discoveries are straightforward to define in single-variable hypothesis testing: a false discovery is one
that is independent of the outcome. In regression models, however, this idea is complicated by the consideration
of various kinds of conditional independence. Typically, in regression a feature Xj is considered a false discovery
if it is independent of the outcome Y given all other features; symbolically, if Xj ⊥⊥ Y |{Xk }k6=j . This is the
perspective adopted by most work in this area, including sample splitting; here, we will refer to this as the
fully conditional perspective. The pathwise approaches use a different definition, which we will call the pathwise
conditional perspective. Letting M(λ) denote the set of variables with nonzero coefficients in the model at the
point in the path where feature j is selected, in these approaches a feature j is considered a false discovery if
Xj ⊥⊥ Y |{Xk : k ∈ Mj− }. Here, λ indexes the the various models along the lasso solution path; see (2.1) for an
exact definition.
This paper introduces a weaker definition of false discovery than those considered in the existing literature.
Rather than the conditional definitions given in the previous paragraph, here we consider a marginal perspective
in which a selected feature j is false if it is marginally independent of the outcome, the definition used in singlefeature testing: Xj ⊥
⊥ Y . Adopting a simpler definition makes it possible to estimate the expected number of
false discoveries as well as their rate, which we call the marginal false discovery rate, or mFDR. Our goal here
is not to argue that the marginal FDR is always superior to either of the conditional FDRs; clearly, there are
times where a conditional FDR is an important quantity of interest. However, the mFDR is also an interesting
and useful summary of feature selection accuracy. This is especially true in high dimensions where, as we will see,
attempting to control a conditional false discovery rate is often hopeless.

2

Marginal false discovery rates

Consider the linear model with normally distributed errors:
y = Xβ + ε
iid

εi ∼ N(0, σ 2 ),
where y denotes the response and X the n × p design matrix, with n denoting the number of independent
observations and p the number of features. Without loss of generality, we assume throughout that the response
and covariates
are centered so that the intercept term can be ignored, and that the features are standardized so
P
b defined as the quantity that
that n1 i x2ij = 1 for all j. We are interested in studying the lasso estimator β,
minimizes
1
ky − Xβk22 + λ kβk1 ,
2n

(2.1)

P
where kβk1 = j |βj |. Here, the least-squares loss is used to measure the fit of the model and the L1 norm is used
to penalize large values of β, with λ controlling the tradeoff between the two.
An appealing aspect of using the lasso to estimate β is that the resulting estimates are sparse: some coefficients
will be nonzero, but for many, βbj = 0. In this paper, we say that a feature for which βbj 6= 0 has been “selected.”
b changes with λ although we will suppress this in the notation; or a large enough value of λ, βbj = 0
Note that β
for all j, but as we lower λ, more variables will be included. To address the expected number of features included
in a lasso model by chance alone, we begin by considering the orthonormal case, then turn our attention to the
general case.

2

2.1

Orthonormal case

b denote the residuals. The lasso (2.1) is a
For a given value of the regularization parameter λ, let r = y − Xβ
convex optimization problem, so the Karush-Kuhn-Tucker (KKT) conditions are both necessary and sufficient for
b
any solution β:
1 0
x r = λ sign(βbj )
n j

for all βbj =
6 0

1 0
x r ≤λ
n j

for all βbj = 0.

Letting X−j and β −j denote the portions of the design matrix and coefficient vector that remain after removing
b denote the partial residuals with respect to feature j. The KKT conditions
the jth feature, let rj = y − X−j β
−j
thus imply that
1 0
x rj > λ
n j
1 0
x rj ≤ λ
n j

for all βbj 6= 0
(2.2)
for all βbj = 0

and therefore the probability that variable j is selected is


1 0
b
P(βj 6= 0) = P
x rj > λ .
n j
This indicates that if we are able to characterize the distribution of n1 x0j rj under the null, we can estimate
the number of false discoveries in the model. Indeed, this is straightforward in the case of orthonormal design
( n1 X0 X = I):
1 0
1
b )
x rj = x0j (y − Xβ + X−j β −j + xj βj − X−j β
−j
n j
n
1
= x0j ε + βj
n
∼ N (βj , σ 2 /n).

(2.3)

Thus, if βj = 0, we have

P


√
1 0
xj rj > λ = 2Φ(−λ n/σ).
n

These results are related to the expected number of false discoveries in the following theorem, the proof of
which follows directly from the above by summing P(βbj 6= 0) over the set of null variables.
Theorem 1. Suppose

1 0
nX X

= I. Then for any value of λ,
√
E |S ∩ N | = 2 |N | Φ(−λ n/σ),

where S = {j : βbj 6= 0} is the set of selected variables and N = {j : βj = 0} is the set of null variables.
To use this as an estimate, the unknown quantities |N | and σ 2 must be estimated. First, |N | can be replaced by
p, using the total number of variables as an upper bound for the null variables. The variance σ 2 can be estimated
by
σ̂ 2 =

r0 r
.
n − |S|

3

Dividing the residual sum of squares by the degrees of freedom of the lasso (Zou et al., 2007) is a simple approach
for estimating the residual variance, but other possibilities exist (e.g., Fan et al., 2012). This implies the following
estimate for the expected number of false discoveries:
√
c = 2pΦ(− nλ/σ̂)
FD
(2.4)
and, as an estimate of the false discovery rate:
c
[ = FD .
FDR
|S|

(2.5)

This straightforward estimate of the false discovery rate is made possible by the orthonormality condition
assumed in (2.3). When the features that comprise X are correlated, however, the distribution of n1 x0j rj is
considerably more complex, as is the proper definition of what is meant by a “false discovery.” Nevertheless, as
the rest of this paper will show, with a suitable definition of false discovery, estimator (2.5) can still be used to
estimate false discovery rates for variable selection in non-orthonormal settings.

2.2

Definition in the general case

Consider the causal diagram presented in Figure 1. In this situation, variable A could never be considered a false
discovery: it has a direct causal relationship with the outcome Y . Likewise, if variable N were selected, this would
obviously count as a false discovery – N has no relationship, direct or indirect, to the outcome.

B
A

N
Y

Figure 1: Causal diagram depicting three types of features and their relationship to the outcome.
Variable B, however, occupies a gray area as far as false discoveries are concerned. From a marginal perspective,
B would not be considered a false discovery, since it is not independent of Y . However, B and Y are conditionally
independent given A, so from a fully conditional regression perspective, B is a false discovery. Lastly, from a
pathwise conditional perspective, B may or may not be a false discovery, depending on whether A is already
included in the model or not at the point at which B enters.
The central argument of this paper is that estimating the number of false selections arising from variables like
B is inherently complicated and requires complex approaches (either in the mathematical or computational sense);
however, simple approaches like the one derived in Section 2.1 may still be used to estimate the number of false
selections arising from variables like N .
Here, I define a noise feature to be a variable like N , that has no causal path (direct or indirect) between it and
the outcome, and the marginal false discovery rate as the proportion of selected features that are noise variables.
Again, this definition is consistent with how false discoveries are defined in univariate testing, but differs from the
existing literature in regression modeling.
It is worth addressing a technical point here: it is possible for two variables to be conditionally dependent despite
being marginally independent. Such a relationship would appear in a directed acyclic graph as the following:

X1

Y
X2

4

As the figure indicates, such a relationship would require a feature (or features) to be caused by the outcome
Y . This is contrary to the usual framework in which regression models are applied, where Y is (or is assumed
to be) a consequence of the features used for prediction (e.g., a genetic variant can lead to heart disease, but
heart disease will not change an individual’s genetic sequence). Thus, although we are actually making a stronger
assumption than mere marginal independence here, in the typical regression setting the two are equivalent. At any
rate, although the term may be slightly imprecise from a technical perspective, we feel that the name “marginal
false discovery” is easily understood and best distinguishes this quantity from false discoveries under the various
conditional perspectives.
The proposed definition has several advantages. First, when two variables (like A and B in Figure 1) are
correlated, it is often difficult to distinguish between which of them is driving changes in Y and which is merely
correlated with Y. As we will see, this causes methods using a conditional perspective to be conservative, especially
in high dimensions.
Second, in many scientific applications, discovering variables like B in Figure 1 is not problematic. For example,
two genetic variants in close proximity to each other on a chromosome will be highly correlated. Although it is
obviously desirable to identify which of the two is the causal variant, locating a nearby variant is also an important
scientific achievement, as it narrows the search to a small region of the genome for future follow-up studies.
The final advantage is clarity of interpretation. Whether or not a feature would quality as a marginal false
discovery depends only on the relationship between it and the outcome, not on whether another feature has been
included in the model. In contrast, interpreting the results from pathway-based tests such as those in Lockhart
et al. (2014) and Tibshirani et al. (2016) can be challenging, as the definition of the null hypothesis is constantly
changing with λ.

3

Independent noise features

The orthonormal conditions of Section 2.1 clearly do not hold in most settings for which penalized regression
is typically used. Thankfully, they can be relaxed in two important ways that make the results more widely
applicable. First, the predictors do not have to be strictly orthogonal in order for the estimator to work; they
can simply be uncorrelated. Second, this condition of being uncorrelated applies only to the noise features – i.e.,
the variables like N in Figure 1; variables like A and B can have any correlation structure. These statements are
justified theoretically in Section 3.1 and via simulation in 3.2.

3.1

Theory

To make these statements concrete, let A, N partition {1, 2, . . . , p} such that βj = 0 for all j ∈ N and the following
condition holds:


1 0
ΣA 0
lim X X =
.
n→∞ n
0 ΣN
Under this definition, the opening remarks of this section can be stated precisely in the following theorem.
√
Theorem 2. Suppose ΣN = I. Then for any j ∈ N and for λn such that the sequence nλn is bounded,
1
d
√ x0j rj −→ N (0, σ 2 ).
n
Theorem 2 shows that if the noise features are uncorrelated, then in the limit n1 x0j rj behaves just as it did in
Section 2.1 (equation 2.3). Thus, estimators (2.4) and (2.5) still provide approximate estimators for the marginal
false discovery rate in this setting, at an accuracy that improves with the sample size. The technical condition
√
b will be √n-consistent. Without it (i.e., for
requiring nλn to be bounded is only necessary so that the estimate β
large values of λ), n1 x0j rj will converge to a random variable with a variance larger than σ 2 due to underfitting.
Note that Theorem 2 treats p as fixed; extending this result to allow p → ∞ would be of interest as future work,
but lies outside the scope of this paper.
5

3.2

Simulation

To illustrate the consequences of Theorem 2, let us carry out the following simulation study, with both a “lowdimensional” (n > p) and “high-dimensional” (n < p) component. Motivated by Figure 1, three types of features
will be included:
• Causative: Six variables with βj = 1
• Correlated: Each causative feature is correlated (ρ = 0.5) with m other features; m = 2 for the lowdimensional case and 9 for the high-dimensional case
• Noise: Independent noise features are added to bring the total number of variables up to 60 in the lowdimensional case and 600 in the high-dimensional case
The causative, correlated, and noise features correspond to variables A, B, and N, respectively, in Figure 1. In each
setting, the sample size was n = 100, while the total number of causative/correlated/noise features was 6/12/42
for the low-dimensional setting and 6/54/540 for the high-dimensional setting. The results of the simulation are
shown in Figure 2.
Type

Actual

Estimated

p: 60

p: 600

Noise features selected

80

60

40

20

0
0

−1

−2

−3

−4

log(λ)

0

−1

p: 60

−2

−3

−4

−3

−4

p: 600

mFDR

0.75

0.50

0.25

0.00
0

−1

−2

−3

−4

log(λ)

0

−1

−2

Figure 2: Accuracy of estimators (2.4) and (2.5) in the case of independent noise features.
As Theorem 2 implies, estimators (2.4) and (2.5) are quite accurate, on average, when the noise features are
independent. The estimated number of noise features selected and the marginal false discovery rate (mFDR) are
6

Group

Causative (A)

Correlated (B)

Noise (N)

Lasso (mFDR)

p: 60

Method

Univariate
LassoCV
Selective inference
Sample splitting
0

5

10

15

20

7.5

10.0

Avg

Univariate

p: 600

Method

Lasso (mFDR)

Selective inference
Sample splitting
0.0

2.5

5.0

Avg

Figure 3: Average number of each type of feature selected by various methods for the simulation setup in
Section 3.2. Cross-validation is omitted from the high-dimensional (p = 600) plot due to the very large number of
noise features it selects, which dominates the plot when included.
both somewhat conservative, as we would expect from using p as an upper bound for the number of noise features
(e.g., in the high-dimensional case, p = 600 but |N | = 540). However, the effect is slight in this setting. For
example, in the high-dimensional case at λ = 0.55, the actual inclusion rate for noise variables was 5%, while the
estimated rate was 6.5%.
Being able to estimate the mFDR means we can use it to select the regularization parameter λ. For example,
\
we could choose λ to be the smallest value of λ such that mFDR(λ)
< 0.1. Figure 3 compares this approach
with several other methods for selecting λ in terms of the number of each type of feature the method selects
on average. For Lasso (mFDR), univariate testing (i.e., marginal regression), sample splitting (using the hdi
package, Dezeure et al., 2015), and the exact conditional path-based test (using larInf and forwardStop from
the selectiveInference package, Tibshirani et al., 2016; G’Sell et al., 2016), the nominal false discovery rates
were set to 10%. For cross-validation, the value of λ minimizing the cross-validation error was selected.
It is worth noting that both Lasso (mFDR) and univariate testing limit the fraction of selections due to noise
features (to 5% and 7%, respectively, in the p = 60 simulation) to the nominal rate, but claim nothing about
the fraction arising from correlated features. This is obvious from the figure for univariate testing; for Lasso
(mFDR), the fraction of selected features coming from either the Correlated or Noise groups (i.e., all features with
βj = 0) was 17% in the low-dimensional setting and 23% in the high-dimensional setting. Nevertheless, compared
to univariate testing, the Lasso (mFDR) approach has two distinct advantages. Figure 3 shows that using a
penalized regression approach both diminishes the number of merely correlated features selected and improves
power to detect the truly causative features.
Among the penalized regression approaches, cross-validation is a considerable outlier, with no protection against
the selection of large numbers of noise features. In the high-dimensional case, cross-validation selected over 30
noise features on average, and is omitted from the plot so as not to obscure the performance of the other methods.
The sample splitting and selective inference approaches, on the other hand, greatly limit the number of both
noise features and correlated features selected by the lasso model. The more restrictive definition of a false
7

discovery used by these approaches, however, causes them to behave quite conservatively, and have considerably
less power to detect the causative features than any of the other methods considered.

4

Correlated noise features

The simulation results of Section 3.2 are something of a best case scenario for the proposed method, since the
variables in N were independent and Theorem 2 establishes that the estimator is consistent in this case. As we will
see, when noise features are correlated, estimator 2.5 becomes conservative. The case of correlated noise features
is significantly less mathematically tractable, making it harder to construct a rigorous proof of this claim, but here
we offer some insight into why this is true and investigate the performance of the proposed method via simulation.
To illustrate why the proposed mFDR estimator becomes conservative when noise features are correlated, let
us consider the p = 2 case. Letting zj = n1 x0j y, we begin by noting that z follows a multivariate normal distribution
with mean 0 and


σ2 1 ρ
,
Var(z) =
n ρ 1
where ρ denotes the correlation between x1 and x2 , and
√
2 |N | Φ(−λ n/σ) = P(|z1 | > λ) + P(|z2 | > λ).

(4.1)

A visual representation of this quantity is given in Figure 4.
The shading represents whether 0 (white), 1 (light gray), or 2 (dark gray) features will be selected under this
orthogonal approximation; the quantity in (4.1) can be found by integrating these constant regions with respect
to the joint density of z.
ρ = 0.5

ρ = −0.5

A
A

z2 0

0

0

λ

B

−λ

B

−λ

0
z1

λ

−λ

0
z1

λ

Figure 4: Illustration of how correlation affects selection in the bivariate case. Features x1 and x2 are positively
correlated (ρ = 0.5) on the left and negatively correlated (ρ = −0.5) on the right.
In the case where x1 and x2 are non-orthogonal, the white region (where no features are selected) remains
the same, but the conditions under which both features are selected (the light and dark gray regions) change.
By working through the KKT conditions for the bivariate case under various conditions, we can determine these
selection boundaries. For example, the boundary for βb2 > 0 given that βb1 > 0 (i.e., given that z1 > λ) is
1
n

x02 r2 > λ =⇒ z2 − ρβb1 > λ
=⇒ z2 > ρz1 + λ(1 − ρ).

8

These boundaries are drawn in black on Figure 4; the preceding equation is the line just above the letter “A”. As the
figure illustrates, when ρ > 0 the region over which β1 , β2 6= 0 narrows in the upper right and lower left quadrants,
and widens in the other two quadrants. Considering the difference E |S ∩ N | − {P(|z1 | > λ) + P(|z2 | > λ)}, we
see that in many regions, the terms cancel out by integrating the same quantity over the same distribution. The
differences lie in four pairs of triangular regions; one such pair is labeled “A” and “B” in the figure. The region
A is where two variables are selected under orthogonality, but only one in the correlated case, while region B is
where one variable is selected under orthogonality, but two are selected in the correlated case. However, because
z1 and z2 are positively correlated in this scenario, the density at each point in A is higher than its corresponding
point in B, and therefore estimator (2.4) is conservative in the sense that it overestimates the expected number of
false discoveries. The opposite scenario happens for ρ < 0, where region B now has the higher probability density,
which again leads to a larger integral under orthogonality, with equality between the two quantities occurring only
at ρ = 0.
In terms of Theorem 2, this argument implies that when ΣN 6= I, the quantity √1n x0j rj converges to a distribution with thinner tails than N(0, σ 2 ). Intuitively, this makes sense: if two noise features are correlated, a
regression-based method such as the lasso will tend to select a single feature rather than both. Thus, the uncorrelated case is not just mathematically convenient, it also represents a worst case scenario with respect to the
number of noise features that we can expect to be selected by chance alone.
To investigate the robustness of the proposed mFDR estimator in the presence of moderate correlation, let
us carry out the following simulation. The generating model contains 6 independent causative features and 494
correlated noise features, with a 1:1 signal-to-noise ratio (n = 100, p = 500, R2 = 0.5). The noise features are
given an autoregressive correlation structure with Cor(xj , xk ) = 0.8|j−k| . The results of the simulation are shown
in Figure 5.
Type

True

Type

Estimated

Estimated

1.00

50
40

0.75

30

mFDR

Noise features selected

True

20

0.50

0.25

10
0

0.00
−0.5

−1.0

−1.5

log(λ)

−2.0

−2.5

−0.5

−1.0

−1.5

log(λ)

−2.0

−2.5

Figure 5: Accuracy of estimators (2.4) and (2.5) in the case of (autoregressive) correlated noise features.
Compared to Figure 2, the mFDR estimates are somewhat more conservative in this case, although still quite
accurate – certainly accurate enough to be useful in practice. For example, at λ = 0.43, the true inclusion rate for
noise variables was 14%, while the estimated rate according to (2.4) was 20%.
This simulation illustrates that although its derivation is based on independent noise features, the proposed
mFDR estimator is reasonably robust to the presence of correlation. Furthermore, to the extent that it is inaccurate, it provides a conservative estimate of the marginal false discovery rate, suggesting (as the conjecture
indicates) that the approach provides control over the mFDR, at least on average.

9

5

Permutation approach

As correlation between noise features becomes more extreme, the proposed estimator becomes increasingly conservative. As an extreme example, suppose that the correlation structure from Section 4 was exchangeable instead of
autoregressive: Cor(xj , xk ) = 0.8 for all j, k. The results of this modification to the simulation (all other aspects
remaining the same) are shown in Figure 6.
As one might expect, the estimates are far more conservative in this case. For example, at λ = 0.43, the
estimated mFDR is 23% even though the true mFDR is only 1%. This is a consequence of the independence
approximation, where (2.5) operates under the simplifying assumption that the selection of one noise feature does
not affect the probability of other noise features being selected. This is reasonably accurate in most cases, but
not in this situation, where all noise features are highly correlated with each other. As a result, the lasso tends
to select only a single noise feature from this highly correlated set, while the independence-based estimate (2.4)
indicates that it has likely selected, say, 7 or 8.
This phenomenon is not unique to penalized regression; substantial correlation among features causes problems
with conventional false discovery rates as well (Efron, 2007). One widely used approach for controlling error rates
while preserving correlation structures is to use a permutation approach (Westfall and Young, 1993), and a similar
strategy may be applied in order to calculate marginal false discovery rates for penalized regression as well. The
primary advantage of this approach is that it eliminates the conservatism of the analytic approach developed in
Section 2, while the primary disadvantage is a greatly increased computational burden. In this section, I describe
two permutation-based methods and apply them to the simulation shown in Figure 6.
1.00

0.75

mFDR

Type
True
Analytic

0.50

PermY
PermR
0.25

0.00
−0.5

−1.0

−1.5

log(λ)

−2.0

−2.5

Figure 6: Accuracy of the analytic (2.5), permute-the-outcome (“PermY”, (5.1)), and permute-the-residuals
(“PermR”) estimators in the case of (exchangeable) correlated noise features.

5.1

Permuting the outcome

The simplest approach is simply to randomly permute the outcome y, creating new outcomes ỹ(b) for b =
e (b) (λ; X, ỹ(b) ), estimate the average num1, 2, . . . , B. Then, for each permutation b, solve for the lasso path β
ber of noise features included in the model for a given value of λ
P
e(b)
b #{βj (λ) 6= 0}
c
FD(λ) =
,
B
and the marginal false discovery rate using
c
\ Y (λ) = FD(λ) ,
mFDR
S(λ)
10

(5.1)

where S(λ) is the number of variables selected by the lasso using the original (i.e., not permuted) data.
This method is applied to the case of heavy exchangeable correlation in Figure 6. The method is considerably
less conservative than the analytic approach, although still conservative for reasons that will be discussed in
Section 5.2. For example, at λ = 0.43, the average estimated mFDR is 4%, where the true value was 1% and the
analytic approach yielded 23%.
By permuting Y , this approach constructs realizations of the data in which all features belong to N , the set of
noise features. Like (2.5), it limits the number of noise features selected but cannot be used to control the number
of false discoveries in the sense of limiting selections from all variables with βj = 0. In the hypothesis testing
literature, this is referred to as “weak control” over the error rate.

5.2

Permuting the residuals

As seen in Figure 6, permuting the outcomes is still conservative in its estimation of mFDR. The primary reason
for this is that, provided that at least some features are not noise, the variance of Y can be separated into signal
and noise; ideally we would permute only the noise, but by permuting the outcome we permute the signal as
well. This has the effect, in a sense, of overestimating the noise present in the model; this is essentially the same
√
phenomenon that occurs in Theorem 2 if the sequence nλn is not bounded.
b
One alternative is, rather than permuting the outcome, to permute the residuals, r(λ) = y − Xβ(λ),
of the
original lasso fit. The method is otherwise identical to that described in Section 5.1, with the notable exception
e (b) (λ; X, r̃(b) ) must be
that the residuals, unlike the outcomes, depend on λ and thus, B separate lasso solutions β
calculated at each value of λ (in Section 5.1, the same solutions could be used for all values of λ), substantially
increasing the computational burden.
Nevertheless, this increased computational cost does offer a benefit, as seen in Figure 6. Unlike the analytic
and permute-the-outcome approaches, permuting the residuals provides an essentially unbiased estimate of the
true false inclusion rate except at very small λ. For example, at λ = 0.29, the average estimated mFDR is 8%, as
was the true mFDR, whereas the permute-the-outcome and analytic methods estimated mFDRs of 16% and 90%,
respectively.
The purpose of this manuscript is primarily to investigate the analytic approach to calculating false inclusion
rates outlined in Section 2, with relatively less emphasis on the permutation approaches. In the opinion of the
author, the analytic approach is almost always useful, since it can be calculated instantly. Whether one wishes to
take the time to carry out the permutation approach, however, depends on context: how correlated the features
are, how problematic a somewhat conservative estimate of mFDR is, and how far along in the analytic process one
is (initial exploration or final publication results). For a more detailed comparison of the analytic and permutation
approaches in the context of genetic association studies, see Yi et al. (2015).

6
6.1

Case studies
Breast cancer gene expression study

As a case study in applying the proposed method to real data, we will analyze data on gene expression in breast
cancer patients from The Cancer Genome Atlas (TCGA) project, available at http://cancergenome.nih.gov.
In this dataset, expression measurements of 17,814 genes, including BRCA1, from 536 patients are recorded on
the log scale.
BRCA1 is a well-studied tumor suppressor gene with a strong relationship to breast cancer risk. Because
BRCA1 is likely to interact with many other genes, including tumor suppressors and regulators of the cell division
cycle, it is of interest to find genes with expression levels related to that of BRCA1. These genes may be functionally
related to BRCA1 and are useful candidates for further studies.
For this analysis, 491 genes with missing data were excluded, resulting in a design matrix with p = 17, 322
predictors. The resulting mFDR estimates for the lasso solution path are presented in Figure 7. Here, the mFDR
estimates indicate that many genes are predictive of BRCA1 expression – we can safely select 55 variables before

11

Selected

Variables selected

0

8

24

58

1.0
# of variables

250

0.8
mFDR

E(False)

134 255

0.6
0.4
0.2
0.0

200
150
100
50
0

−1.0

−2.0

−3.0

−4.0

−1.0

log(λ)

−2.0

−3.0

−4.0

log(λ)

Figure 7: False inclusion rate estimates for a lasso model applied to the breast cancer TCGA data.
the false inclusion rate exceeds 10%. This makes sense scientifically, as a large number of genes are known to affect
BRCA1 expression through a variety of mechanisms, and the sample size here is sufficient that we should be able
to identify many of them. For the sake of comparison, the univariate approach of separately testing each feature
identifies 7,903 genes that have a significant correlation with BRCA1 at a false discovery rate of 10%.
In contrast, the selective inference and sample-splitting procedures each select just a single feature. Although
the mFDR estimates are slightly conservative as discussed in Section 4, this matters far less in practice than
whether one adopts a marginal or conditional perspective with respect to false discoveries. This example clearly
illustrates how conservative the conditional perspective is in practice, especially with high-dimensional data.
The mFDR approach is also vastly more convenient than sampling splitting or covariance testing from the
standpoint of computational burden. Calculating the mFDR is essentially instantaneous after fitting the lasso
path; the entire analysis requires just 1.3 seconds. Meanwhile, sample splitting required 9 minutes, and applying
the pathwise conditional test with the forwardStop rule using the selectiveInference package required 1.5 hours.
Interestingly, although sample splitting and pathwise conditional testing each identified a single feature, they
did not select the same feature. Sample splitting selected the gene NBR2, which is unquestionably associated
with BRCA1 expression – the two genes are adjacent to one another on chromosome 17 and share a promoter (Di
et al., 2010). However, NBR2 is the third gene to enter the lasso model. The forwardStop rule based on pathwise
conditional testing recommends stopping after the first variable is added to the model, and thus fails to identify
NBR2.
Finally, it is worth comparing these results to the selection of λ by cross-validation (CV). For the TCGA data,
λ = 0.0436 minimizes the CV error. The estimated mFDR at this value, however, is 77%, indicating that although
this value of λ may be attractive from a prediction perspective, we cannot be confident that the variables selected
by the model are truly related to the outcome. Indeed, it has long been recognized from a theoretical perspective
that while the lasso has attractive variable selection and prediction properties, it cannot achieve both those aims
simultaneously (Fan and Li, 2001). The mFDR estimates illustrate this concretely: λ = 0.0436 produces accurate
predictions, but a larger value, λ = 0.0627, is required in order to have confidence that no more than 10% of the
variables selected are false discoveries.

6.2

Minimax concave penalty

The mFDR estimator follows directly from the KKT conditions for a given penalty, and it is straightforward to
extend to other penalties. In fact, the KKT conditions for many penalties lead to the same expression (2.2) and
therefore the same mFDR estimator (although as a technical note, for nonconvex optimization problems these
conditions are not usually referred to as KKT conditions; they are still necessary but no longer sufficient).
One such penalty is the minimax concave penalty, or MCP (Zhang, 2010). Briefly, the MCP produces sparse
estimates like the lasso, but modifies the penalty such that the selected variables are estimated with less shrinkage
towards zero (see the original paper for details). The main consequence of this is that MCP solutions tend to be
more sparse than those of the lasso, in that they allow features to have larger regression coefficients, requiring

12

fewer features to achieve the same predictive ability. To put it differently, for the lasso to estimate large coefficients
accurately, it must lower λ, thereby increasing the number of noise features in the model and increasing the mFDR.
MCP, on the other hand, can estimate large coefficients accurately while still screening out the noise features, at
least asymptotically (this is known as the “oracle property”).
To see how this works in practice, we can fit an MCP model to the TCGA data and compare its results to those
of the lasso. For the value of λ minimizing CV error, both methods had similar predictive accuracy (cross-validated
R2 = 0.61 for the lasso and 0.58 for MCP), but the MCP model used far fewer features (38 compared to 96 for the
lasso). Consequently, for these values of λ the estimated mFDR for MCP was much lower than that of the lasso:
5% compared to the lasso’s 77%. If we restrict the lasso mFDR to 5% by increasing λ, its predictive accuracy falls
to R2 = 0.57.
This is representative of the relationship between lasso and MCP: the lasso can typically achieve slightly better
prediction accuracy, but in doing so selects a large number of noise features. This pattern has been observed in
simulation studies (e.g., Breheny and Huang, 2011), but mFDR estimates offer a way to observe and assess this
tradeoff in the analysis of real data.
This example also illustrates the ease with which the proposed estimator can be extended to other penalties.
In addition to MCP, the SCAD (Fan and Li, 2001) and elastic net (Zou and Hastie, 2005) also have similar KKT
conditions leading to (2.2) and thus require only trivial modifications to the mFDR estimator.

6.3

Genetic association study of cardiac fibrosis

As an additional example, let us also look at a genetic association study of cardiac fibrosis. The data come
from the Myocardial Applied Genomics Network (MAGNet), which collected tissue and gene expression data on
313 human hearts along with genetic data using the Affymetrix Genome-Wide Human SNP Array 6.0. We are
interested in the the ratio of cardiomyocytes to fibroblasts in the heart tissue. An abundance of fibroblasts is
indicative of cardiac fibrosis, which leads to heart failure. In this analysis, the goal is to discover single nucleotide
polymorphisms (SNPs) associated with increased fibrosis. Here, we use the log of the cardiomyoctye:fibroblast
ratio as the response, and SNPs with < 10% missing data as the predictors, resulting in p = 660, 496 features.
mimimac2 was used to impute the missing genotypes that remained (Fuchsberger et al., 2015). Neither sample
splitting or covariance testing was attempted here due to the size of the data set.
In contrast to the TCGA gata, no features can be selected with any degree of confidence in the MAGNet data.
For all values of λ in which features are selected, the mFDR estimate is 100%. This is consistent with a traditional
genome-wide association analysis, which also fails to identify any SNPs that are significant following a Bonferroni
correction for multiple testing.
This negative example illustrates the usefulness of mFDR estimates in terms of guarding against false positives.
The efficiency with which methods like the lasso extend to high dimensions make analyses like this (with n = 313
and p = 660, 496) feasible from a computational standpoint. However, there are important statistical considerations
here that need to be accounted for; namely, with p so large, the probability of selecting noise features is so great
that one cannot place any trust in the variables that the lasso selects. These considerations are clearly illustrated
by the estimation of mFDR, but are lost if one simply uses the lasso to identify the top |S| SNPs (as in, e.g., Wu
et al., 2009).

7

Discussion

The false inclusion rate estimator presented in this article is a straightforward, useful way of quantifying the
reliability of feature selection for penalized regression models such as the lasso. Compared with other proposals
such as sample splitting and the pathwise covariance and selective inference tests, mFDR uses a more relaxed
definition of a false discovery, in which only variables marginally independent of the outcome are considered
to be false discoveries. The relaxed definition pursued here has several advantages in practice, particularly in
high dimensions. These advantages include greater power to detect features, no added computation cost, and a

13

simple rate with a straightforward interpretation. Although the estimate can be conservative for highly correlated
features, in most realistic settings this conservatism is mild.
The proposed estimator is easy to implement and is available (along with the permutation methods of Section 5)
in the R package ncvreg (Breheny and Huang, 2011), which was used to fit all of the models presented in this paper
(except where the hdi and selectiveInference packages were used). The following bit of code demonstrates its
use:
fit <- ncvreg(X, y, penalty="lasso")
obj <- fir(fit)
plot(obj)

This will fit a lasso model, estimate the mFDR at each value of λ and assign it to an object, then plot the results,
as in Figure 7.
The simplicity of the method makes it available at no added computational cost and very easy to generalize to
new methods: the fir function in ncvreg also works for elastic net, SCAD, MCP, and Mnet (Huang et al., 2016)
penalties. As with any statistic, there are limitations and one needs to be careful with interpretation, but mFDR
is a convenient summary measure that is helpful in the selection of λ for penalized regression models as well as a
useful estimate for the reliability of feature selection in such models.

Acknowledgments
I would like to thank Jian Huang for fruitful discussions of this idea when it was still in its early stages, Ina
Hoeschele for discussions of the permutation approach, and Ryan Boudreau for the MAGNet data of Section 6.3.

Appendix
All of the code and data to reproduce the analyses in this manuscript, with the exception of the human genetic
data considered in Section 6.3, are available at https://github.com/pbreheny/lassoFDRpaper.
b , for any j ∈ N we have
Proof of Theorem 2. Starting with the observation that rj = Xβ + ε − X−j β
−j
√
1
1
b )}
√ x0j rj = √ x0j ε + ( n1 x0j X−j ){ n(β −j − β
−j
n
n
The first term, √1n x0j ε, follows a N (0, σ 2 ) distribution by the independence of xj and ε. The second term, n1 x0j X−j ,
converges to zero by the assumption that ΣN = I (xj is also independent of variables in A since j ∈ N ). Finally,
√
the third term is bounded in probability provided that nλn is bounded (Fan and Li, 2001). Therefore, by
Slutsky’s Theorem, the entire quantity converges in distribution to N (0, σ 2 ).

References
Breheny, P. and Huang, J. (2011). Coordinate descent algorithms for nonconvex penalized regression, with
applications to biological feature selection. Annals of Applied Statistics, 5 232–253.
Dezeure, R., Bühlmann, P., Meier, L. and Meinshausen, N. (2015). High-dimensional inference: Confidence
intervals, p-values and R-software hdi. Statist. Sci., 30 533–558.
Di, L.-J., Fernandez, A. G., De Siervi, A., Longo, D. L. and Gardner, K. (2010). Transcriptional
regulation of brca1 expression by a metabolic switch. Nat Struct Mol Biol, 17 1406–1413.
Efron, B. (2007). Correlation and large-scale simultaneous significance testing. Journal of the American Statistical Association, 102 93–103.
14

Fan, J., Guo, S. and Hao, N. (2012). Variance estimation using refitted cross-validation in ultrahigh dimensional
regression. Journal of the Royal Statistical Society Series B, 74 37–65.
Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle properties. Journal
of the American Statistical Association, 96 1348–1360.
Fuchsberger, C., Abecasis, G. R. and Hinds, D. A. (2015). minimac2: faster genotype imputation. Bioinformatics, 31 782–784.
G’Sell, M. G., Wager, S., Chouldechova, A. and Tibshirani, R. (2016). Sequential selection procedures
and false discovery rate control. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 78
423–444.
Huang, J., Breheny, P., Lee, S., Ma, S. and Zhang, C.-H. (2016). The Mnet method for variable selection.
Statistica Sinica, 26 903–923.
Lockhart, R., Taylor, J., Tibshirani, R. J. and Tibshirani, R. (2014). A significance test for the lasso.
The Annals of Statistics, 42 413–468.
Meinshausen, N., Meier, L. and Bühlmann, P. (2009). p-values for high-dimensional regression. Journal of
the American Statistical Association, 104 1671–1681.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society
Series B, 58 267–288.
Tibshirani, R. J., Taylor, J., Lockhart, R. and Tibshirani, R. (2016). Exact post-selection inference for
sequential regression procedures. Journal of the American Statistical Association, 111 600–620.
Wasserman, L. and Roeder, K. (2009). High dimensional variable selection. Annals of Statistics, 37 2178.
Westfall, P. H. and Young, S. S. (1993). Resampling-based multiple testing: Examples and methods for
p-value adjustment. Wiley, New York.
Wu, T., Chen, Y., Hastie, T., Sobel, E. and Lange, K. (2009). Genome-wide association analysis by lasso
penalized logistic regression. Bioinformatics, 25 714–721.
Yi, H., Breheny, P., Imam, N., Liu, Y. and Hoeschele, I. (2015). Penalized multimarker vs. single-marker
regression methods for genome-wide association studies of quantitative traits. Genetics, 199 205–222.
Zhang, C. (2010). Nearly unbiased variable selection under minimax concave penalty. Annals of Statistics, 38
894–942.
Zou, H. and Hastie, T. (2005). Regularization and variable selection via the elastic net. Journal of the Royal
Statistical Society Series B, 67 301–320.
Zou, H., Hastie, T. and Tibshirani, R. (2007). On the “degrees of freedom” of the lasso. Annals of Statistics,
35 2173–2192.

15

