1

A knockoff filter for high-dimensional selective inference
Rina Foygel Barber∗ and Emmanuel J. Candès†

arXiv:1602.03574v2 [stat.ME] 18 Sep 2017

February 2016; Revised September, 2017

Abstract
This paper develops a framework for testing for associations in a possibly high-dimensional linear model
where the number of features/variables may far exceed the number of observational units. In this framework,
the observations are split into two groups, where the first group is used to screen for a set of potentially relevant
variables, whereas the second is used for inference over this reduced set of variables; we also develop strategies
for leveraging information from the first part of the data at the inference step for greater power. In our work,
the inferential step is carried out by applying the recently introduced knockoff filter, which creates a knockoff
copy—a fake variable serving as a control—for each screened variable. We prove that this procedure controls
the directional false discovery rate (FDR) in the reduced model controlling for all screened variables; this says
that our high-dimensional knockoff procedure ‘discovers’ important variables as well as the directions (signs)
of their effects, in such a way that the expected proportion of wrongly chosen signs is below the user-specified
level (thereby controlling a notion of Type S error averaged over the selected set). This result is non-asymptotic,
and holds for any distribution of the original features and any values of the unknown regression coefficients, so
that inference is not calibrated under hypothesized values of the effect sizes. We demonstrate the performance
of our general and flexible approach through numerical studies. Finally, we apply our method to a genome-wide
association study to find locations on the genome that are possibly associated with a continuous phenotype.
Keywords: knockoffs, inference in high-dimensional regression models, errors of type S, Lasso, variable selection, (directional) false discovery rate, martingales.

1

Introduction

Many modern studies in the sciences take the following form: collect a large amount of data, and then ask which
of the potentially many measured variables are possibly associated with a response of interest. A typical example
would be a genome-wide association study (GWAS), where one wishes to ‘discover’ which genetic variants are
possibly associated with a trait. In such studies, it is common to read genetic variants using single-nucleotide polymorphism (SNP) arrays and then look for associations between the trait and the hundreds of thousands of SNPs.
In these types of studies and others, we typically have many more variables or features (SNPs) than observations
(individuals in the study).
Finding statistically significant associations between a response and a large set of potentially explanatory variables
requires the specification of a statistical model; in this work, we consider the classical linear regression model,
which takes the form
p
X
Y =
βj Xj + ,
(1)
j=1

where Y is the response variable of interest, the Xj ’s are the explanatory variables, and  is a (stochastic) error
term. In the GWAS example, Y may be the level of HDL cholesterol and Xj the number of recessive alleles at a
given location on the genome. In such applications it is important to keep in mind that the linear relationship (1) is
just an approximation of a more complicated and unknown model.
∗ Department

of Statistics, University of Chicago
of Statistics and Mathematics, Stanford University

† Departments

2

02/2016; Revised 09/2017

In classical statistics, when we speak of finding associations, we think of testing each of the p hypotheses Hj :
βj = 0, and one can see that there are a few complications with this viewpoint. The first is that for fixed designs,
the linear model (1) with p parameters is not identifiable whenever the number n of samples is smaller than this
number p of features. Formally, testing an unidentifiable model may seem like a problematic proposition. A second
complication is that in our GWAS example, it may very well be the case that causal mutations affecting levels of
HDL cholesterol, for instance, have not been typed (i.e. the array of SNPs is a subsample of the genome, and
does not include the causal mutation). In such studies one is, therefore, led to search for indirect effects. In other
words, while the causal factors have not been measured, we nevertheless hope to have measured variables that are
sufficiently close or correlated to be picked up by association. In these settings, it is understood that even if there
is a true sparse ‘causal’ model, we may end up working with an approximate and non-sparse statistical model of
the form (1) in which few, if any, of the features have precisely zero effect (i.e. for most features j, βj = 0 is not
exactly true); this phenomenon occurs because the causal effects have been folded onto the observed variables.
There are two consequences of this: (1) sparsity assumptions, which are traditionally assumed in the literature to
relieve ourselves from the identifiability problem, are not applicable1 and (2) when few or none of the regression
coefficients take on the value zero, controlling a classical Type I error typically makes little sense since nearly any
selected feature is technically a true positive, and we need a different measure of performance to capture the idea
that we would like to correctly identify the meaningful large effects while screening out features with near-zero
effects. To summarize our discussion, we have at least two problems:
• Lack of identifiability (at least when p > n).
• Possibly correlated features and lack of model sparsity so that few, if any, of the null hypotheses actually
hold.
These complications should not occlude the relatively simple goal of a statistical selection procedure, which in
the GWAS context can be stated as follows: understanding that SNPs may be proxies for causal mutations, we
would like to find as many proxies or locations on the genome without too many false positives—a false positive
being a reported location (SNP) with no causal mutation sitting nearby. Being guaranteed a sufficiently small
fraction of false positives would assure the scientist that most of her discoveries are indeed true and, in some
sense, replicable. We are thus interested in reliable methods for finding the large effects, those SNPs with large
coefficient magnitudes. This paper develops a framework and methods for getting close to this goal. At a high
level, we propose a two-stage procedure: first, using a portion of the data set to screen for a large set of potentially
relevant features (which will likely include a high proportion of false positives), and second, applying the recently
developed knockoff filter to the second portion of the data in order to perform inference within the submodel
selected at the screening stage.

1.1

Errors of Type S and directional FDR

In a series of papers, Gelman and his collaborators introduce a useful measure they call error of Type S (S stands
for sign). In their work [10], they argue that “in classical statistics, the significance of comparisons (e. g. , θ1 − θ2 )
is calibrated using Type I error rate, relying on the assumption that the true difference is zero, which makes no
sense in many applications.” They then propose “a more relevant framework in which a true comparison can be
positive or negative, and, based on the data, you can state “θ1 > θ2 with confidence”, “θ2 > θ1 with confidence”,
or “no claim with confidence”.” In this framework, a Type S error (sometimes called a Type III error in other
works) occurs when we claim with confidence that the comparison goes one way when, in fact, it goes the other
way. This echoes an earlier point made by Tukey and we quote from [31]: “All we know about the world teaches
us that the effects of A and B are always different—in some decimal place—for any A and B. Thus asking ‘Are
the effects different’ is foolish. What we should be answering first is ‘Can we tell the direction in which the effects
of A differ from the effects of B?’.” We refer to [31] as well as [18] and [5] for additional discussion.
This point of view is extremely relevant to our (possibly high-dimensional) regression problem and our goal of
finding the important effects. When we choose to report a variable, we should be able to state with confidence
1 The

literature also typically assumes low correlations between variables, which is also not applicable in the GWAS example.

3

02/2016; Revised 09/2017

its direction of effect, i.e. whether βj > 0 or βj < 0. To press this point, we would surely have little faith in a
procedure that would control errors of Type I but would not be able to tell the directions of the effects. Another
useful aspect of Type S errors is that small effects βj ≈ 0 are generally those for which we cannot have much
certainty about their direction, since the size of the effect is comparable to or below the noise level of our estimate.
Therefore, if we were concerned by signed errors, we would probably not report these effects as discoveries (even
if we believe that no effect is exactly zero)—arguably a wise thing to do.
To measure our success at selecting only those effects that are large enough to be meaningfully distinguished from
noise, we might like to control the mixed directional false discovery rate (FDRdir ) [25, 6] defined as follows: letting
d j ∈ {±1} of the direction of effect,
Sb ⊂ {1, . . . , p} be the set of selected features together with estimates sign
FDRdir = E [FDPdir ] ,

where the mixed directional false discovery proportion (FDPdir ) is given by
n
o
d j 6= sign(βj )
j ∈ Sb : sign
FDPdir =
b ∨1
|S|

(2)

with the convention that sign(0) = 0. This definition subsumes Type I and Type S errors: a false discovery or error
occurs either when a zero effect is selected, or when a nonzero effect is selected but with the incorrect sign. (In
b ∨ 1 denotes max{|S|,
b 1}; this maximum is taken so that if we make zero discoveries, then the
the denominator, |S|
FDP is calculated as zero.) The term “mixed” comes from this combination of two error types; for our purposes,
we consider declaring that βj > 0 to be a sign error whether the truth is that βj = 0 or βj < 0, and will not
distinguish between these two types of errors, and so we will drop the term “mixed” and refer to it simply as the
“directional FDR”. The directional FDP is then the total number of errors of this (combined) type averaged over
the selected set. In regression settings in which many of the βj ’s are approximately zero but perhaps not exactly
zero, controlling the directional FDR may be more appropriate than the “classical” FDR (which counts Type I
errors only):
n
o
j ∈ Sb and βj = 0
FDR = E [FDP] , FDP =
.
b ∨1
|S|
Note that controlling the directional FDR has a different flavor than classical FDR, since paraphrasing Gelman
et al., the significance is not calibrated using Type I errors which implicitly assumes that some of the βj ’s are
zero. In contrast, we are after a form of inference holding no matter the values of the βj ’s. In classical settings
ind

with independent statistics, where yi ∼ N (µi , 1), say, and we wish to test for the means, [6] establishes that the
Benjamini-Hochberg procedure [4] achieves directional FDR control.
In earlier work [1] we introduced the knockoff filter, a new variable selection procedure, which rigorously controls
the classical FDR, in the low-dimensional setting where p ≤ n. It should be clear that the equivalent notion (2)
associated with Type S errors is in principle more difficult to control since by definition we have
FDRdir ≥ FDR,
due to the fact that the classical FDR does not record an error when the jth feature is correctly included in the
model but with the incorrect sign, while the directional FDR does. A first contribution of this paper is to show
that the knockoff filter offers more than advertised in our earlier work [1]: indeed, we will find that this filter,
with no modifications, actually controls the directional FDR as well as the FDR. In particular, this implies that the
knockoff filter is applicable in settings where we do not expect any sparsity but nonetheless wish to test whether
our conclusions regarding effect signs are reliable. This addresses one of the complications discussed earlier.

1.2

High dimensionality

To address the other complication, namely, the high-dimensionality issue (p > n), we propose the following
general strategy:

4

02/2016; Revised 09/2017

• Screening step. In a first step, screen all the features {Xj }, j = 1, . . . , p, as to identify a set Sb0 ⊂ {1, . . . , p}
of potentially relevant features with |Sb0 | < n. In this paper, we shall regard this step as being quite liberal in
the sense that it typically produces a long list, hopefully containing most of the important features (with large
effect sizes) but also possibly many features with zero effects βj = 0 or nearly vanishing effects βj ≈ 0.
• Inference/selection step. The screening step yields a reduced model
X partial
Xj + ,
βj
Y =

(3)

b0
j∈S

where we use the notation β partial to indicate that both the definition and the meaning of the regression
coefficients has changed (see below). Then test for associations in this reduced model by controlling the
directional FDR.
In our GWAS example, this means that we would first screen for promising SNPs, and then extract from this set a
final and smaller list of SNPs we are confident are associated with the phenotype under study (recall that reporting
a SNP only means that we are confident a mutation lies nearby). This approach—a screening step, followed by
inference on the screened submodel—has been studied by many researchers; see [35, 33] for a non-exhaustive list
of important works among those lines.
We would like to emphasize the importance of considering the directional FDR rather than the classical (unsigned)
FDR in the framework of a screening step followed by an inference step. Imagine that our sampled observations
y ∈ Rn have mean Xβ and uncorrelated errors. In the reduced model, we wish to provide inference about the
b
partial regression coefficients, β partial ∈ R|S0 | , when we regress the response y onto the features XSb0 only. Even
if the original coefficient vector β for the full regression model is exactly sparse with many entries exactly equal
to zero, the vector of partial regression coefficients β partial is nonetheless likely to be dense unless there is special
structure within X (e.g. if columns of X are orthogonal to each other).
Now, consider a feature Xj not appearing in the original full model, i.e. βj = 0. In the partial regression, we may
find one of two scenarios:
• The coefficient in the partial regression may be large, that is, βjpartial 6≈ 0. This typically occurs if Xj is highly
correlated with some strong signal Xk that was missed by the screening step. (For instance, we may have
missed a causal SNP but included a close neighbor, instead, in the screened set.) Here, we would generally
prefer to include Xj in the selected model in any case, since it is a good proxy for the missed relevant feature
Xk .
• Alternately, the coefficient might remain near zero, βjpartial ≈ 0. This is likely whenever Xj is not a proxy for
any missed feature. Here, the sign of βjpartial cannot be estimated with much certainty, and so with directional
FDR control, we are likely to exclude Xj from our final model.
These considerations hopefully make clear that careful screening followed by directional FDR control might yield
a valuable selection procedure, which in our GWAS example would likely be able to report SNPs (locations on the
genome) associated with a trait while controlling a very natural notion of false positive rate. On the other hand,
a setting where this framework does not apply is that of causal inference, in which one would like to report true
causal effects and would prefer not to select variables that, due to correlation, act as proxies to the true causal
features.

1.3

Comparison with other works

There is of course a vast literature on inference in the regression setting. When n > p, we have available the
sampling distribution of least-squares estimates, which makes the construction of p-values and confidence intervals
possible, although selective inference procedures controlling the FDR and similar criteria are more tricky because
of the correlations between such test statistics. For an overview of some of the work when n > p, we refer the
reader to the companion article [1].

5

02/2016; Revised 09/2017

There has also been much recent work on the high-dimensional setting p > n, which is the focus of our work here.
One type of approach to high-dimensional inference is to avoid the identifiability problem by making assumptions
on the data distribution such as (1) a highly sparse coefficient vector β, and/or (2) a “beta-min” condition which
requires nonzero entries of β to not be too close to zero, and/or (3) assumptions on the design matrix such as lowpairwise correlations or placing a random distribution on the features. These types of conditions ensure that, with
high probability (asymptotically), a sure screening property will hold: that is, at any place in the method requiring a
model selection step, with high probability all the relevant features (i.e. all Xj with βj 6= 0) will be selected. Under
these conditions, it is possible to produce an asymptotically normal estimator for each coefficient βj which can
then yield a p-value or confidence interval which is asymptotically valid. Methods in this category include the work
of Voorman et al. [32], Belloni et al. [3], and Huang et al. [13]. A different type of result, relying on rapid decay
of the error in estimating β in an asymptotic setting instead of requiring a sure screening property, is the debiasing
work of Zhang and Zhang [37] and Javanmard and Montanari [17]; see also Voorman et al. [32] for a different but
related approach. Lockhart et al. [21] propose a sequential testing procedure, which moves along the Lasso path as
variables are added into the model, under the assumption that the signals appear before the false positives along the
Lasso path (which holds with high probability asymptotically under similar assumptions, i.e.. highly sparse β, the
beta-min condition, and a well-behaved design matrix). To contrast these lines of work with our method presented
here, empirically we observe reliable performance in terms of FDR control at even relatively small sample sizes,
and our theoretical framework avoids strong assumptions on signal strength or sparsity level or the design matrix.
In contrast with the asymptotic setting described above, other approaches focus on finite-sample p-values and
confidence intervals which do not rely on “high probability” selection events but instead are calculated by conditioning on the screened or selected model. Lee et al. [19] develops exact post-selection inference for the Lasso
when applied with a fixed penalty parameter, with methodology that extends to other model selection procedures
which can be expressed as linear constraints on y; Tian et al. [28] extend this methodology to the setting of the
square-root-Lasso, where the noise level σ is no longer assumed to be known. In contrast, Tibshirani et al. [30]
propose methods for sequential selection procedures, such as the path of selected models obtained by tuning the
Lasso penalty parameter; in this setting, FDR control can be obtained with sequential hypothesis testing procedures e.g. [12], and the tests can be viewed as asking whether the most recent added feature has a true effect in the
partial model fitted so far. Fithian et al. [9] develops general theory for this type of inference after conditioning on
(sequential or nonsequential) selection events, termed “data carving”. In contrast to these methods, which perform
inference on the outcome of specific model selection methods (such as the Lasso), the knockoff framework gives
the flexibility of choosing (nearly) any test statistics for model selection.
In earlier work, Leeb and Pötscher [20] studied the impossibility of characterizing the (unconditional) distribution
of post-selection estimators, and the PoSI (“post-selection inference”) framework of Berk et al. [7] proposed handling the range of possible selected models by taking a maximum over the set of possibilities. In practice, this type
of method may have high power when considering extremely small submodels, but lose power when the selected
model being considered grows larger.
Finally, we have called attention on Type S errors as a useful error metric for controlled variable selection. (Methods such as Lee et al. [19]’s post-selection inference for the Lasso can provide p-values testing a one-sided
(i.e. signed) hypothesis, and can therefore be used to control directional FDR, but since this literature focuses
on inference for individual hypotheses rather than FDR control, this distinction is not discussed explicitly.) We
believe this is a correct framework for the model selection problem in the high-dimensional setting, where strict
sparsity assumptions may be violated whenever a reduced model is considered.

2

Knockoffs

Since this paper builds upon and extends the knockoff methodology introduced in [1], we first review the essentials
of this method. Imagine we have a response y ∼ N (Xβ, σ 2 I) and have reasons to believe that the coefficient
vector β ∈ Rp is sparse (so that, for now, we restrict our attention to controlling Type I errors only, and do not yet
consider sign errors). To then estimate β, the literature offers many tools for estimating a sparse model in a highdimensional setting. For instance, we may use a forward stepwise selection method such as orthogonal matching
pursuit (OMP) [22]. We initialize with the empty selected model, S0 = ∅, and a zero estimate, b(0) = 0. Then at

6

02/2016; Revised 09/2017

each time k ≥ 1, we choose feature jk = arg maxj |X>
j (y − Xb(k − 1))|, add this feature to the set by setting
Sk = Sk−1 ∪ {jk }, and refit the coefficients via least squares,
b(k) = arg minb∈Rp ky − Xbk2

subject to bj = 0 for all j ∈
/ Sk .

Fixing some desired number of steps K, or some stopping criterion for choosing K adaptively, we then take our
estimated coefficient vector to be β̂ = β̂OMP (K) := b(K). Alternately, we might choose a convex optimization
method such as the Lasso [29], where

β̂Lasso (λ) = arg minb∈Rp 21 ky − Xbk22 + λkbk1
(4)

A critical question is this: can we estimate the number of false discoveries in the coefficient vector produced by
OMP or by the Lasso, i.e. in β̂ = β̂OMP or β̂ = β̂Lasso ? That is to say, can we guess the number of coordinates
with β̂j 6= 0 but βj = 0? The knockoff filter provides answers to questions of this kind by constructing a knockoff
companion for each variable, which essentially acts as a control.

2.1

Knockoff variables

e j is constructed such that the augmented n × (2p) design matrix of features
For each feature Xj , a knockoff copy X
together with knockoffs satisfies
"
# 

>
>e

> 

X
X
X
X
Σ
Σ − diag{s}
e
e
X X
X X = e>
,G,
(5)
e >X
e = Σ − diag{s}
Σ
X X X

where the vector s ≥ 0 of course must be small enough that the resulting Gram matrix is positive semidefinite.2
By definition, we see that the knockoffs exhibit the same correlation structure as the original variables, and that
the cross-correlation is also preserved in the sense that for all pairs of distinct variables j 6= k, we have X>
j Xk =
>
e
Xj Xk . This property has the important consequence that if j is null, i.e. βj = 0, then
d e>
X>
j y = Xj y,

d

(6)

where = denotes equality in distribution (to see this, note that both random variables are Gaussian with the same
mean and variance). The distributional equality (6) also holds more generally: for any subset S of nulls (i.e. a set
of features j with βj = 0), we have



d 
e >
e > y,
X X
y= X X
(7)
swap(S)





e
e by swapping the columns Xj and X
e j for each j ∈ S. In
where the matrix X X
is obtained from X X
swap(S)
fact, one can say that knockoffs are constructed according to (5) precisely for (7) to hold. To see the value of the
exchangeability property (7), suppose we were using marginal correlations to test for nonzero effects. Since the
knockoffs were man-made with a total disregard for the response, they are known to be null and hence we could
use them as controls. For instance, to declare that βj 6= 0 with some confidence we would certainly want to see
e>
e>
|X>
j y| > |Xj y|. In other words, the values of Xj y provide meaningful comparisons since they have the same
distributions as the corresponding null statistics.
We can now return to our problem of estimating how many of the discoveries in our sparse coefficient vector
estimate β̂ are, in fact, null. To do this, we run the OMP, the Lasso, or any desired sparse estimation method on
the original design augmented with knockoffs. For instance, for OMP, at each time k we would choose a feature




e > (y − X X
e · b(k − 1)) ,
jk = arg max X X
j
j=1,...,2p

where b(k − 1) ∈ R2p is the solution from the previous time step, while for the Lasso we would optimize


e bk2 + λkbk1 .
minb∈R2p 12 ky − X X
2
2 If

e can be constructed to satisfy this matrix equality as long as n ≥ 2p.
s is chosen so that the Gram matrix is positive semidefinite, then X

7

02/2016; Revised 09/2017

In either case, we can denote the solution by [β̂ β̃] (the first p components are the coefficients of the original

 



e > X X
e and X X
e > y are invariant
variables and the last p are for the knockoffs). Due to (5) and (7), X X
and distributionally invariant, respectively, to any swap of original null features and knockoffs. As a consequence,
for any subset of nulls,
d
[β̂ β̃]swap(S) = [β̂ β̃].
Hence, β̂j and β̃j have identical distributions whenever j is null and, in particular,
d

|{null j : β̂j 6= 0}| = |{null j : β̃j 6= 0}| ≤ |{j : β̃j 6= 0}|.
Therefore, by looking at the number of knockoff variables that are in our estimated coefficient vector, we have an
upward biased estimate of the number of null variables selected by this estimator.

2.2

The knockoff filter

After constructing knockoff copies of the features, we can then use them to create a feature selection rule that will
achieve exact FDR control. (One can also control other types of error such as the family-wise error rate (FWER)
or the k-FWER as in [14]). For each variable j ∈ {1, . . . , p}, we work with a general statistic Wj obeying a
sufficiency and an antisymmetry property. The sufficiency property says that Wj depends upon the design matrix
only through its covariance, and upon the response y through marginal correlations: namely,

 
 

e > X X
e , X X
e > y)
Wj = wj ( X X

(8)

e j has the effect of changing the
for some function wj . The antisymmetry property says that swapping Xj and X
sign of Wj ,
(


e , y),


wj ( X X
j 6∈ S,
e


wj ( X X
,
y)
=
(9)
swap(S)
e , y), j ∈ S.
−wj ( X X

In [1], a much discussed example of such a statistic concerns knots at which variables enter the Lasso path: for
e j , record the first time that this feature or its knockoff enters the Lasso path,
each feature Xj and each knockoff X
i.e. the largest penalty parameter value λ such that β̂j 6= 0 or β̃j 6= 0; then set
(


e j,
+1, if Xj enters before X
e
(10)
Wj = largest λ such that Xj or Xj enters Lasso path ·
e
−1, if Xj enters after Xj .

e j is selected by the OMP method or another
We can analogously record the first time that a variable Xj or X
e
forward stepwise algorithm (setting |Wj | to be large if Xj or Xj enters early). Many other statistics are of course
possible, for instance
Wj = |β̂j | − |β̃j |
(11)
where [β̂ β̃] ∈ R2p is the estimated coefficient vector from any penalized method of the form {min 12 ky −


e bk2 + P (b)} (as long as P (b) treats entries j and j + p of the coefficient vector b exchangeably).
X X
2

By construction, for a null feature, the original variable is equally likely to enter before or after its knockoff. In
general, this says that the sign of a null statistic Wj is a fair coin flip, and it also holds that the signs of the null
Wj ’s are independent of each other as well as of the magnitudes |W|.
Now that we have test statistics for each variable, we need a selection rule. Choose a threshold T > 0 by setting3


#{j : Wj ≤ −t}
T = min t > 0 :
≤q ,
(12)
#{j : Wj ≥ t}

3 In practice, we will only choose a value t from the set of magnitudes of W. Abusing notation, we write min{t > 0 : . . . } to mean
min{t ∈ W+ : . . . }, where W+ = {|Wj | : |Wj | > 0, j = 1, . . . , p}.

8

02/2016; Revised 09/2017

where q is the target FDR level. The output of the procedure is the selected model
Sb = {j : Wj ≥ T } .

(13)

This selection makes sense since the numerator in (12) is an (over)estimate for the number of false discoveries
among all reported variables with Wj ≥ t, that is,
d

#{j : Wj ≤ −t} ≥ #{null j : Wj ≤ −t} = #{null j : Wj ≥ t}.
Hence, the ratio appearing in the right-hand side of (12) is an estimate of the false discovery proportion (FDP)
if we were to use the threshold t. Hence, we can view our selection rule as a step-up rule stopping the first time
our estimate falls below our target level. A slightly more conservative procedure, the knockoff+ filter, is given by
incrementing the number of negatives by one, replacing the threshold in (12) with the choice


1 + #{j : Wj ≤ −t}
≤q ,
(14)
T+ = min t > 0 :
#{j : Wj ≥ t}
and setting

Sb = {j : Wj ≥ T+ }.

(15)

The key results in [1] state that knockoff+ controls the FDR, and the knockoff controls a slight modification of the
FDR, at the level q. These results hold without any assumptions on the design or the value of the unknown regression coefficients. Crucially, we emphasize that neither the knockoff procedure nor its FDR controlling property
assume any knowledge of the noise level σ.

3

Controlling Errors of Type S

Although they are more difficult to control, we have argued that Type S errors are more meaningful for regression
problems since we would like to be able to tell the direction of effect with confidence (and should avoid reporting
those variables whose direction we cannot reliably determine). Consider the setting of Section 2 where n ≥ 2p
and y ∼ N (Xβ, σ 2 I). Then to estimate the direction of effect—the sign of βj —note that
⊥
⊥
e j )> y ∼
(Xj − X
N (sj βj , 2sj σ 2 ) for j = 1, . . . , p,

where s ≥ 0 is from the knockoff construction (5). Hence, a natural estimate for the sign of βj is


e j )> y .
d j = sign (Xj − X
sign

(16)

With this, the knockoff filter with no adjustment whatsoever controls the directional FDR as well as the FDR. To
understand intuitively why this error measure, which in principle is harder to control than FDR, is still bounded
by the same method, note that for a nonnegative effect (βj ≥ 0), the probability of choosing a negative sign,
d j = −1, is in fact highest when βj = 0; this error is less likely if βj > 0 since in that case we would have
i.e. sign
e j )> y be normally distributed with a positive mean.
(Xj − X

Theorem 1. Assume y ∼ N (Xβ, σ 2 I) and fix any desired FDR level q. With the estimated direction effects (16),
knockoff+ controls the directional FDR (2); that is, FDRdir ≤ q for the selected set Sb defined in (15). If we use
knockoff instead, then mFDRdir ≤ q, where mFDRdir is the slightly modified version defined as
o 
 n
d j 6= sign(βj )
j ∈ Sb : sign
 ,
mFDRdir = E 
(17)
b + q −1
|S|

where the selected set Sb is defined as in (13).

9

02/2016; Revised 09/2017

This result (proved in Appendix A.1) is more powerful than Theorems 1 and 2 from [1], which established that
knockoff+ controls the FDR at level q, i.e. FDR ≤ q, and knockoff controls the modified FDR, i.e. the expected
b + q −1 (this denominator is
ratio between the number of (unsigned) false discoveries |{j ∈ Sb and βj = 0}| and |S|
as in (17)).
Our new Theorem 1 is also more delicate to prove. The reason is that when we work with Type I errors as in the
original result, the probability of claiming that βj 6= 0 is calculated under the assumption that βj = 0. When
working with Type S errors, however, there is no such calibration since there is no null, and each feature Xj might
have a nonzero coefficient βj , although it may be close to zero. In that case, the test statistic Wj will have a
different probability of appearing with a positive or negative sign, and we are no longer dealing with i.i.d. and
unbiased signs for the null Wj ’s.
Thus in the proof of Theorem 1, we will see that the result is ultimately a consequence of the following martingale
inequality regarding a sequence of Bernoulli variables which have varying probabilities of success:
Lemma 1. Suppose that B1 , . . . , Bn are independent variables, with Bi ∼ Bernoulli(ρi ) for each i, where
mini ρi ≥ ρ > 0. Let J be a stopping time in reverse time with respect to the filtration {Fj }, where
Fj = {B1 + · · · + Bj , Bj+1 , . . . , Bn } .
Then
E



1+J
1 + B1 + · · · + BJ



≤ ρ−1 .

The proof of this lemma is given in Appendix A.4. In the simpler setting where ρ1 = · · · = ρn , this result is
proved (in a slightly different form) in [1], with the proof relying heavily on the exchangeability of the Bj ’s (when
the Bj ’s are determined by the signs of the null Wj ’s, each has distribution Bernoulli(0.5)). The result stated here
is more subtle due to the lack of exchangeability.
To see the connection between Lemma 1 and (directional) FDR control for the knockoff+ method, observe that
FDP =

#{j null : Wj ≥ T+ }
1 + #{j null : Wj ≤ −T+ }
#{j null : Wj ≥ T+ }
=
·
1 ∨ #{j : Wj ≥ T+ }
1 + #{j null : Wj ≤ −T+ }
1 ∨ #{j : Wj ≥ T+ }
#{j null : Wj ≥ T+ }
,
≤q·
1 + #{j null : Wj ≤ −T+ }

where the inequality follows from the definition of T+ (14). Now reorder the indices of the null Wj ’s so that
|W(1) | ≥ |W(2) | ≥ . . . , where |W(1) | is the null with the largest magnitude, |W(2) | the second largest, and so on,
and let Bj = 1W(j) <0 . Then we have
#{j null : Wj ≥ T+ }
(1 − B1 ) + · · · + (1 − BJ )
1+J
=
=
−1,
1 + #{j null : Wj ≤ −T+ }
1 + B1 + · · · + BJ
1 + B1 + · · · + BJ
where J is the index such that |W(1) | ≥ · · · ≥ |W(J) | ≥ T+ > |W(J+1) | ≥ . . . ; this last expression is the
quantity in Lemma 1. If each Wj , for a null feature, is equally likely to be positive or negative, then the Bj ’s
are i.i.d. Bernoulli(0.5) variables, and we obtain FDR control at level q by applying Lemma 1 with ρi = 0.5.
However, when considering directional FDR control, the βj ’s may be nonzerzo and so the Wj ’s will typically not
be symmetric; therefore we may need to apply Lemma 1 with varying ρi ’s. We refer to Appendix A.2 for details.

4

Knockoffs in High Dimensions

In high dimensions, where p > n, the knockoff construction is no longer possible—in fact, the Gram matrix
e j = Xj for each feature j = 1, . . . , p, so that the
condition (5) would be possible only if s = 0, that is, if X
knockoff procedure would have zero power. In this setting, one straightforward approach would be to use part
of the data to reduce the number of features, and a disjoint part of the data to run the knockoff filter. In the next
section, we develop this data splitting approach, and then find that we can gain substantial power by using a subtle
way of ‘recycling’ the data.

10

02/2016; Revised 09/2017

4.1

Feature screening and reusing data

Consider splitting the n observations into two disjoint groups of size n0 and n1 = n − n0 , used for the purpose of
first screening for a smaller set of potentially relevant features, then running a model selection procedure over this
reduced list of features, respectively. We denote the two disjoint portions of the data as (X(0) , y(0) ) ∈ Rn0 ×p ×Rn0
and (X(1) , y(1) ) ∈ Rn1 ×p × Rn1 , and then follow the next two steps:
• Screening step: using (X(0) , y(0) ), we identify a subset Sb0 ⊂ [p] of potentially relevant features, such that
|Sb0 | < n1 .

• Selection step (splitting): ignoring any features that were discarded in the screening step, we run the knock(1)
off procedure on the remaining data, that is, on (XSb , y(1) ).
0

This straightforward data splitting approach is a natural extension of the low-dimensional knockoff filter, and it is
clear that the approach will control the directional FDR in the final model selection step as long as the screening
step correctly captures all the relevant features—those with non-vanishing regression coefficients—a property often
referred to as sure screening in the literature [8]. (False positives in the screening step, i.e. including too many
null features, do not pose a problem). If the sure screening property fails, we can still obtain inference guarantees
relative to the new, screened submodel—see Section 4.3 below.
There is an inherent loss of power due to the split of the data, since the model selection step uses n1 rather than n
observations. However, using disjoint data sets for the screening and selection steps is critical, since the distribution
of y(0) cannot be treated as a Gaussian linear model after the screening step has taken place. The screening step
is a function of the random variable y(0) , and so the following modification of the selection step would not control
the FDR: ignoring any features that were discarded in the screening step, we run the knockoff procedure on the full
data set, i.e. on (XSb0 , y). The loss of FDR control is not merely theoretical: a null feature Xj that is chosen by
the screening step is generally more likely to appear as a false positive when running the knockoff filter, leading to
a much higher FDR.
In light of this, we propose two mechanisms for increasing power relative to the data splitting procedure described
so far:
1. data recycling, where the first portion of the split data can be reused to some extent without losing any of
the guaranteed FDR control,
2. and signed statistics, where we can decide adaptively to focus our search on a positive effect only or a
negative effect only.
4.1.1

Increasing power with data recycling

Surprisingly, data recycling retains the FDR control properties of the split-data procedure, while raising power
substantially to approach the sensitivity of the full-data procedure.
e (1) for
• Selection step (with recycling): we begin by constructing a knockoff matrix on the remaining data X
b
S
0

(1)

XSb just as we would in the data splitting version. The difference now is that we concatenate the original
0
design matrix on the first n0 observations with the knockoff matrix on the next n1 observations,
" (0) #
XSb
e
0
XSb0 =
.
(18)
e (1)
X
b
S
0

(Note that on the first part of the split data, knockoffs are exact copies.) We then run the knockoff filter on the
e b , and the original
entire data set with all n samples, using the design matrix XSb0 , the knockoff matrix X
S0
response y.

11

02/2016; Revised 09/2017

(0)

The term “recycling” refers to the way that we incorporate the data (XSb , y(0) ), which was already used in the
0
screening step, for the selection step. Here, we think of the mechanism for choosing Sb0 as being completely
arbitrary as opposed to having some kind of ‘pre-registered’ screening method. For this reason, we need to treat
y(0) as a fixed vector (or rather, to condition on its value) in the construction and analysis of the knockoff filter for
the selection step (see Section 4.1.3 for additional discussion about this point).
Of course, knockoffs with data splitting are actually a special case of knockoffs with data recycling. Any knockoff
(1) e (1)
(1)
statistics Wj that are a function of XSb , X
(i.e. data splitting) can trivially be expressed as a function
b0 , and y
S
0
e b , and y by simply ignoring the first n0 many data points.
of X b , X
S0

S0

Increased power Why do we expect the knockoff method to exhibit higher power when using data recycling
rather than data splitting, when we see in (18) that the knockoff matrix differs from the original features only on
the second part of the split data? Indeed, later on in Section 5 we will see empirically that the power gain can be
quite substantial.

To understand this phenomenon, we return to the mechanism of the knockoff filter itself. Recall that we construct
e j ), where |Wj | is large if either Xj or X
e j appear highly
a statistic Wj for each feature/knockoff pair (Xj , X
significant, and sign(Wj ) indicates which of the two appears most significant. For instance, if we are using a
e j was one of the
forward stepwise selection method as described in Section 2, |Wj | will be large if either Xj or X
e
e j.
first variables to enter the model; we will have Wj > 0 if Xj enters before Xj , and Wj < 0 if Xj enters after X

If the early portion of the knockoff path (where we order the features in order of magnitudes |Wj |) has many
positive Wj ’s and few negative Wj ’s, the stopping rule (12) will allow us to choose a threshold T that is not too
large, and we will make many rejections. If instead the high-magnitude |Wj |’s have many negative signs, though,
we will be forced to stop early, choosing a high threshold T and making few rejections. For any null feature j
(i.e. βj = 0), the sign of Wj is equally likely to be negative as positive, so in order to obtain high power, we need
two things:
1. Good separation in the feature ordering—most of the features that appear early in the path, i.e. most j with
|Wj | large, are non-nulls;

2. Good power in the non-null signs—for most non-nulls, the statistics have a positive sign, Wj > 0.

Next, we ask how data splitting and data recycling compare in light of these two considerations. For the second,
both methods suffer from power loss relative to an unscreened method—specifically, obtaining Wj > 0 rather than
Wj < 0 depends on the model selection method being able to distinguish between feature j and its knockoff copy.
For data splitting, by the sufficiency (8) and antisymmetry (9) properties, we see that if
(1) > (1)

Xj

y

e (1) > y(1)
and X
j

(19)

are equal in distribution, then sign(Wj ) is equally likely to be +1 or −1; in other words, it is only if these two inner
products have substantially different distributions, that we can hope to see sign(Wj ) = +1 with high probability.
For data recycling, the same statement holds, except instead we compare the inner products
(0) > (0)

X>
j y = Xj

y

(1) > (1)

+ Xj

y

and

e > y = X(0) > y(0) + X
e (1) > y(1) .
X
j
j
j

(20)

Examining (19) versus (20), we see that these two comparisons are equivalent, as (20) simply adds the same
(0)
constant (i.e. Xj > y(0) ) to each term of (19). In other words, these two tasks are equally difficult.
However, for the first component of power—that is, good separation in the feature ordering—we can expect that
knockoffs with data recycling will have the advantage, since the first n0 data points carry substantial information
separating signals from nulls. Of course, if some non-null feature j has strong correlation with the response,
e j towards the top of the path—since
then the information from the first n0 data points will push both Xj and X
(0)
(0)
e
e (1) 6= X(1) , will hopefully enable Xj
X
= Xj by construction—but the second portion of the data, where X
j
j
j
e j (i.e. Wj > 0).
to enter the path before X
To summarize, the knockoff filter’s power relies on the non-null Wj ’s ability to appear early and with positive sign
along the path of statistics; data recycling helps the non-null Wj ’s appear early, thus increasing power.

12

4.1.2

02/2016; Revised 09/2017

Increasing power with signed statistics

Since the knockoff filter with data splitting or with data recycling treats the first portion of the data, (X(0) , y(0) ), as
fixed once the screening step has been performed, we are free to use this data in any way we wish, to try to gather
information about the true signals. While the screening step identifies the indices of the variables which are likely
(0)
to contain signals, we can also use this step to identify probable signs—that is, if feature Xj is selected in the
screening step, which is run on the data set (X(0) , y(0) ), we can at the same time specify the sign of the estimated
(0)
d j . We then use this information to look specifically for an effect consistent with
effect of this feature, denoted sign
this sign in the second phase of our procedure. This will in general increase power, since we are extracting more
information from the first part of the data before running the selection procedure on the second part, (X(1) , y(1) ).
Example: penalized least-squares To illustrate the use of the sign information in practice, suppose that we are
using a penalized least-squares regression method for computing the statistics Wj —for instance, as described in
e j entered the Lasso path (or the path of
Section 2, the Wj ’s may indicate the time at which features Xj and X
some other sequential selection procedure), or may be given by Wj = |β̂j | − |β̃j | where [β̂ β̃] is the solution to a
penalized least squares optimization problem, with some penalty function P (b) (with the Lasso as a special case if
P (b) = λkbk1 ). In either case, if we are not using sign information, we might begin by computing the penalized
least-squares solution [β̂ β̃] via either

for data splitting, or

minb∈R2|Sb0 |
minb∈R2|Sb0 |

1
(1)
2 ky
1
2 ky

(1) e (1)
2
− [XSb X
b ]bk2 + P (b)
S
0

0

e b ]bk2 + P (b)
− [XSb0 X
2
S0

for data recycling. Instead, to make use of the sign information gathered at the screening phase, we can consider a
sign-restricted version of this penalized least-squares problem: in the data splitting version we would be interested
in the solution to
(1) e (1)
1
2
(1)
minb∈R2|Sb0 |
− [XSb X
b ]bk2 + P (b)
2 ky
S
s.t.

(0)

0

0

bj · signj ≥ 0
(0)
bj+|Sb0 | · signj ≥ 0

and similarly for the recycling version. In other words, we are running the same penalized least-squares optimization, but with the added restriction that we will only select the jth feature or jth knockoff feature if its estimated
effect direction agrees with the sign information gathered at the screening stage.
Increased power Why should we expect the knockoff method to exhibit higher power when leveraging the sign
(0)
information, signj for each selected feature j ∈ Sb0 ? Intuitively, we think of the screening step as choosing from
p possible hypotheses, where hypothesis Hj is described by the question, “Does Xj appear in the true model?”.
However, we can instead frame this step as considering twice as many hypotheses: hypothesis Hj+ is the question,
“Does Xj appear in the true model with a positive effect?”, and hypothesis Hj− is the question “Does Xj appear in
the true model with a negative effect?”. With this in mind, if the screening step fits a model to the first part of the
data (X(0) , y(0) ), and this fitted model includes the jth feature with a positive coefficient, then this indicates that
hypothesis Hj+ is more likely to contain a signal. Therefore, it makes sense that in our selection step we should
focus our attention on hypothesis Hj+ . Of course, it is possible that the true effect is in fact negative, but since this
is less likely, we would generally be increasing the noise if we give equal attention to hypotheses Hj+ and Hj− .
4.1.3

Relation to existing work

We pause here to compare our data recycling technique to the data carving methods of Fithian et al. [9]. In that
work, as in ours, some part of the data (X(0) , y(0) ) is used to perform an initial screening step, and is then reused
for inference along with the remaining data (X(1) , y(1) ). However, the mechanism behind reusing the first part of

13

02/2016; Revised 09/2017

the data is highly different. In Fithian et al. [9]’s work, the first part of the response y(0) is reused by leveraging
the “remaining randomness” after the screening step: y(0) is regarded as a random vector, with distribution given
by conditioning on the outcome of the screening step, i.e. (y(0) | Sb0 ). In contrast, in our “recycling” procedure,
y(0) can be used in an arbitrary way for the screening step, and so there is no “remaining randomness” in y(0) .
Instead, we reuse y(0) by treating it as a fixed vector; y(1) is the only variable treated as random for the purpose of
our FDR control results.
More broadly, the knockoff filter with either data splitting or data recycling, can be considered as similar in flavor
to the “screen and clean” methodology of Wasserman and Roeder [33], where an initial screening step (performed
via a high-dimensional Lasso, with a tuning parameter chosen by validation), is followed by an inference step on
an independent portion of the data, with inference performed via a least-squares regression. This methodology
has been applied to genome-wide association studies (GWAS) by Wu et al. [35]. At a high level, our method
follows this same overall framework of screening for a low-dimensional submodel, then using new data for lowdimensional inference; however, we will gain power by “recycling” the first part of the data, and more significantly,
by using the knockoff filter rather than least-squares for the inference step, which in the low-dimensional setting
gives substantial gains in power as it leverages the sparse structure of the true model [1].

4.2

Sure screening and directional FDR control

In order to demonstrate the correctness of the selection step with recycling, we first assume that the screening
step is highly likely to find all true signals—that is, Sb0 ⊇ support(β) with high probability. Note that we do
not assume that the screening step exactly selects the true support—while this type of sparsistency result holds for
the Lasso in the asymptotic setting (see e.g. [38]), in practice, even with ideal simulated data, it is generally the
case that we cannot perfectly separate true signals from false positives using the limited data available, unless the
signal is extremely strong and extremely sparse [26]. Instead, this assumption merely requires that there is some
sufficiently liberal screening procedure that is likely to capture all the true signals, along with many false positives.
From a theoretical standpoint, this property is known to hold under conditions far weaker than those needed for
exact recovery of the true support [36, 8].
Suppose we fix some chosen method for the screening step, which is constrained only in that it must be a function
of the first portion of the data, (X(0) , y(0) ). Define the sure screening event as
E = {Sb0 ⊇ support(β) and |Sb0 | ≤ n1 /2} ,

and note that 1E is a function of y(0) when we treat X and β as fixed. Conditioning on this event, directional FDR
control of the selection step holds by the following theorem (proved in Appendix A.1):
Theorem 2. Suppose y ∼ N (Xβ, σ 2 I). Then the knockoff procedure (using either data splitting or data recyd j as in (16), controls the modified directional FDR at the level
cling), with estimated signs sign
E [mFDRdir | E] ≤ q ,

while if knockoff+ is used, then the directional FDR is controlled as
E [FDRdir | E] ≤ q .
In particular, if E occurs with probability near 1, then the various forms of the directional FDR are controlled even
without conditioning on E, by reformulating the results as FDRdir ≤ q + P {E c } in the case of knockoff+ and
mFDRdir ≤ q + P {E c } for knockoff.

4.3

Directional FDR control in the reduced model

If the screening step misses variables from the model, a bias is introduced which creates a fundamental difficulty.
To discuss this, recall that knockoffs are constructed to guarantee the crucial pairwise exchangeability for the nulls
(7). Among other things, we have seen that exchangeability has the consequence that null variables are equally

14

02/2016; Revised 09/2017

likely to be selected as their knockoff companions. When nonzero effects have not been screened, this may no
longer be the case. Concretely, imagine that Xk has a nonzero effect but was not identified in the screening step.
e>
Then since the knockoff construction used only the screened set Sb0 , we cannot guarantee that X>
j Xk = Xj Xk
for j ∈ Sb0 , which in turn may lead to
d

e>
X>
j y 6= Xj y

even when j is a null feature, i.e. βj = 0. For example, if Xj were strongly correlated with a true signal Xk ,
e j . In general, and in
then we could easily imagine that Xj would be more likely to be selected than its knockoff X
contrast to the lower dimensional setting, the bias implies that when we count the number of knockoffs selected at
some threshold level, we cannot be sure that this number is a good (over)estimate of the number of selected nulls.
(As explained before, however, it may be desirable to select Xj in this case since it is a proxy for Xk .)
We now reframe this issue in terms of the partial regression coefficients—that is, we will perform inference on
the coefficients in the reduced model defined by the subset of features Sb0 selected in the screening step, a point of
view in line with much of the recent literature discussed earlier [19, 7, 11, 9]. Thus we are interested in the partial
regression coefficients—the coefficients of y regressed onto only the selected features. Throughout this section,
we will consider a random design model, and will define the population-level partial regression coefficients as
h
i
(1)
(1)
(1)
(1)
β partial = (XSb > XSb )−1 XSb > E y(1) XSb ; X(0) , y(0) ,
(21)
0

0

0

0

(1)

that is, the expected coefficients when y(1) regressed on the selected features XSb . In other words, this is the
0
quantity for which the least-squares solution in the reduced model is an unbiased estimate. Here, expectation is
taken conditional on the first part of the data (i.e. (X(0) , y(0) ), the data used for the screening step), and on the
(1)
selected features XSb .4
0

FDR control with a Gaussian design In this paper, we give results for a random feature model X, which
assumes that the rows X[i] of X are drawn i.i.d. from a multivariate Gaussian distribution,

iid
X[i] ∼ N ν, Ψ ;

(22)

here, the parameters (i.e. the mean vector ν ∈ Rp and the covariance matrix Ψ ∈ Rp×p ) are arbitrary and
completely unknown to us. The following theorem (proved in Appendix A.1) guarantees directed FDR control in
this setting.
Theorem 3. Assume that the rows of X follow (22), where ν ∈ Rp and Ψ ∈ Rp×p are arbitrary and unknown.
Define the expected partial regression coefficients β partial as in (21). Then the knockoff-with-recycling procedure
controls the modified directional FDR at the level
mFDRdir ≤ q ,
while if knockoff+ is used, then the directional FDR is controlled as
FDRdir ≤ q ;
d j are defined as in (16), and a false discovery is any j ∈ Sb such that sign
d j 6=
here, the estimated signs sign
partial
sign(βj
).

4 This definition is slightly different from the notion of partial regression coefficients found in much of the selective inference literature
−1 X> E [y | X]. We could alternately
discussed here. i.e. [19, 7, 11, 9], which treats the design X as fixed and studies β partial = (X>
b0 )
b XS
b0
S
h
i−1 h S0 i
> y . The difference is between these various
consider a population-level version of these coefficients, β partial = E X>
X
E
X
b
b
b
S
S0

0

S0

definitions is typically not substantial, and indeed under many random models for the design X, can be proved to be vanishingly small.

15

02/2016; Revised 09/2017

While the Gaussian design assumption is strong, we believe that FDR control would be maintained across a far
broader range of models.5 We emphasize that we do not assume knowledge of the mean and covariance parameters
of the Gaussian distribution of X, and do not need to estimate these parameters in the method; in fact the Gaussian
assumption is used only for the technical details of the proof, to show that we can model the response y(1) as
y(1) = XSb βSb0 + XSbc βSbc + (1) = XSb βSpartial
+ 0(1) ,
b
(1)

(1)

0

0

(1)

0

0

0

where 
is i.i.d. Gaussian noise (with a larger variance than the original noise , due to the noise added from
the true signals missed by the screening step). This is where the fact that knockoffs do not need the value of the
noise level (in low dimensions) is immensely useful: the knockoff method will provide valid inference no matter
the value of the ‘new noise level’.
0(1)

5

Simulations

We now examine the performance of our method on a simulated high-dimensional data set, and compare to several
other techniques for inference in low-dimensional and high-dimensional regression.

5.1

Data

The size of our simulated data set is given by dimensions n = 2000 and p = 2500. We generate a true vector
of coefficients β ∈ Rp by choosing k = 50 locations at random for the true signals, with coefficients βj chosen
randomly from {±4.5}, and the other p − k entries set to zero. We then generate a design matrix X ∈ Rn×p whose
rows are i.i.d. draws from a N (0, Σ) distribution, where Σ is defined by a tapered structure, Σij = ρ|i−j| , for a
correlation parameter ρ ∈ {0, 0.25, 0.5, 0.75}; we then normalize the columns of X to have unit norm. For each ρ,
we generate a single matrix X and then generate 100 i.i.d. draws of the error  ∼ N (0, In ), and set y = Xβ + .
All procedures are run with target FDR level q = 0.2.

5.2

Methods

Our simulated experiments compare the following methods:
• The Benjamini-Hochberg (BH) [4] procedure (applied to the least-squares regression coefficients, after
screening for a low-dimensional submodel).
• Exact selective inference for the square-root Lasso, developed by Tian et al. [28]—this method builds on the
earlier work of Lee et al. [19] developing selective inference for the Lasso.
• Knockoff filter with data splitting or data recycling (applied after screening for a low-dimensional submodel).
We now give details for the implementation of each method.
Benjamini-Hochberg We use n0 = 750 data points for an initial screening step, which is carried out by solving
the Lasso across the entire range of λ values, from λ = ∞ to λ = 0. We then choose a screened set Sb0 consisting
of the first 450 features entering the path, and record the sign sj ∈ {±1} for feature j ∈ Sb0 (i.e. the sign of its
estimated coefficient, at the time when it first enters the Lasso path). We next use the remaining n1 = 1250 data
(1)
(1)
(1)
points for inference. We calculate the least-squares regression coefficients β̂LS = (XSb > XSb )−1 XSb > y(1) , and
0
0
0
compute t-scores
(β̂LS )j
q
tj =
(1)
(1)
σ̂ · (XSb > XSb )−1
jj
0

5 An

0

additional result appears in an earlier draft of this work, available at https://arxiv.org/abs/1602.03574v1 (see Theorem
3), treating the case where the design matrix X is fixed and sure screening is not assumed, but instead assuming that the noise level σ is known.

16

02/2016; Revised 09/2017

for each j ∈ Sb0 . Above, σ̂ 2 is the classical estimate of variance,
(1)

2

σ̂ =

ky(1) − XSb β̂LS k22
0

n1 − |Sb0 |

.

Under the sure screening property (where Sb0 contains all true signals, as discussed in Section 4.2), for each null
feature j, tj would follow a t-distribution with n1 − Sb0 degrees of freedom. We thus convert each t-score to a
one-sided p-value Pj (with the side chosen based on the sign sj from the screening step). Finally, we run the
Benjamini-Hochberg procedure [4] at the level q to select a model.
Selective inference for the square-root Lasso The square-root Lasso [2] is given by
 >

kX gk∞
β̂ = arg min {ky − Xbk2 + λ · kbk1 } where λ = κ · Eg∼N (0,In )
kgk2
b∈Rp

(23)

for a constant κ typically chosen in [0.5, 1]. The solution β̂ is in fact exactly equal to the solution to the Lasso (4)
at some different penalty parameter value; the benefit of the square-root Lasso is that the penalty parameter λ can
be chosen without reference to the unknown noise level σ 2 .
Tian et al. [28] derive exact post-selection confidence intervals and p-values for the coefficients βj , after conditioning on the selection event—that is, after observing the signed support of β̂, the solution to the square-root
Lasso—based on the truncated t distribution, without knowledge of the noise level σ. (This work builds on the
polytope method developed in Lee et al. [19], Tibshirani et al. [30] for the Lasso and other procedures.) An implementation of this method, in Python, is available in [27], where the inference is performed approximately by
estimating the noise σ and then using truncated normal distributions, as in Lee et al. [19]’s inference for the Lasso.
We use this available code for our experiments, in order to calculate a one-sided selective p-value Pj for each
feature selected by the square-root Lasso (based on the sign sj with which feature j was selected). We then apply
the Benjamini-Hochberg procedure at level q to the resulting set of p-values in order to select a final model. The
only tuning parameter in this method is the constant κ appearing in (23). We test κ ∈ {0.5, 0.6, . . . , 1} and find
the highest power at κ = 0.6; only this value of κ is displayed in our results.6
Knockoff filter We use n0 = 750 data points for the screening step, carried out exactly as for the BH procedure,
resulting in a submodel Sb0 and signs sj for each j ∈ Sb0 .7 We next use the remaining n1 = 1250 data points for
the knockoff filter, using sign information as described in Section 4.1.2, and either data splitting or data recycling.
Our statistics Wj are formed using the square-root Lasso. Specifically, for data splitting, we solve a sign-restricted
square-root Lasso,
n
o
 (1) (1) 
e
[β̂ β̃] =
arg min
ky(1) − XSb X
bk
+
λ
·
kbk
,
(24)
2
1
b
S
b

b∈R2|S0 |
bj ·sj ≥0 for all j = 1, . . . , p
bj+|Sb | ·sj ≥0 for all j = 1, . . . , p

0

0

0

 (1) (1) 
e
where the penalty parameter λ is defined as in (23) (but with XSb X
b0 in place of X). For data recycling, we
0


 S(1)

e
e (1) and y(1) . In each case, we then
optimize the same problem but with X b X b and y in place of X
X
S0

b0
S

S0

define Wj = |β̂j | − |β̃j | for each j ∈ Sb0 , and then apply the knockoff filter.

b0
S

We run each method with κ ∈ {0.1, 0.2, . . . , 1} and find the highest power at κ = 0.5 for data splitting and κ = 0.7
for data recycling; only these values of κ are displayed in our results.
6 The

kX> gk

software [27] sets λ as κ times the 95th percentile of kgk ∞ in lieu of what is diplayed in (23).
2
7 We note that the two-stage knockoff method of course requires us to choose the data split proportions, i.e. n and n = n − n , as well
0
1
0
b0 |; in our simulated data experiments, changing these values across a reasonable range (e.g. n0 = n/2
as the size of the screened submodel, |S
versus n0 = n/4) leads to mild changes in performance.

17

02/2016; Revised 09/2017

5.3

Results

The knockoff filter and the BH procedure both work within the reduced model defined by the subset of features
Sb0 chosen at the initial screening step; the inference guarantees of these methods are with reference to this large
submodel. On the other hand, the p-values computed by Tian et al. [28]’s method are calculated with reference to
a relatively small submodel, given by features selected by the square-root Lasso. Therefore, in comparing against
Tian et al. [28]’s method, it does not make sense to compare their FDRs over the partial regression coefficients,
since these partial regression coefficients are defined within vastly different submodels. Instead, for each correlation setting and for each method tested, we compute the following measures of performance (averaged over 100
iterations):
• Directional false discovery rate relative to the true full model: the average proportion of the selected features
that were selected with an incorrect sign; and
• Power relative to the true full model: the average proportion of the true support, support(β), which is
selected by the method.
We remark that, in this experiment, the sure screening event discussed in Section 4.2 typically does not occur; even
with uncorrelated features (ρ = 0), we usually screen around 80% of the true signals. Nonetheless, in our results,
the directional FDR and the ordinary, non-directional FDR coincided for all methods—that is, these methods made
only Type I, not Type S, errors. This is an artifact of the clean simulated data, where there are no weak nonzero
signals—any Type S (rather than Type I) error would require that, for some feature j with βj = +4.5 we incorrectly
conclude that βj < 0 (or vice versa), which of course is highly unlikely.
For all methods, the full-model directional FDR and power, over the varying correlation parameter ρ, are displayed
in Table 1 and in Figure 1. We see that at ρ = 0 all four methods control the FDR successfully, and indeed some
are conservative, but all the methods increase in FDR as the correlation ρ grows, with least-squares + BH showing
substantially higher FDR than the target level q. The knockoff methods show slightly higher power than the other
two methods, with a slight gain when using data recycling rather than data splitting.
Next, we consider assessing our false discoveries in the reduced model, that is, Type I and Type S errors for the
partial regression coefficients β partial as defined in (21). When ρ = 0 and the features are, therefore, independent
from each other, the population-based definition of β partial in (21) means that βjpartial = βj ; that is, testing the partial
regression coefficients is equivalent to testing the coefficients in the full model. This is why in this setting there
are many features with exactly zero effect as before, and there are no weak nonzero signals; we again expect to
see directional FDR equal to FDR. On the other hand, for ρ ∈ {0.25, 0.5, 0.75}, we will see βjpartial 6= βj for some
features j, because the positive correlation between nearby features (at the population level) means that, if a true
signal k is missed while a feature j near k is selected at the screening stage, then Xj will act as a proxy for Xk
in the reduced model. This then means that βjpartial can be small but nonzero in the reduced model, and so we may
plausibly make some errors of Type S. (In fact, the reader may verify that in our setup, βjpartial will generally not
vanish—even if βj = 0 in the full model—whenever j is the element from Sb0 closest to k from either the left or
the right.)
In Table 2, we see the FDR and directional FDR for the screened methods—that is, for least-squares + BH, and for
knockoffs with data splitting and with data recycling. As expected, the FDR and directional FDR are equal when
ρ = 0, while when ρ > 0, the FDR is generally slightly smaller than the directional FDR. In all settings and for all
methods, both measures of FDR are below the target level q = 20%.

6

Real Data Experiments

We next implement our high-dimensional knockoff method on a genome wide association study (GWAS) data set.
The data comes from the Northern Finland Birth Cohort 1966 (NFBC1966) [24, 16], made available through the
dbGaP database (accession number phs000276.v2.p1). Here, 5402 SNP arrays from subjects born in Northern
Finland in 1966 are recorded, as well as a number of phenotype variables measured when the subjects were 31
years old.

18

02/2016; Revised 09/2017

Least-squares + BH
(after screening)
Selective inference for
sqrt Lasso [28]
Knockoffs (with
data splitting)
Knockoffs (with
data recycling)

ρ
Dir. FDR
Power
Dir. FDR
Power
Dir. FDR
Power
Dir. FDR
Power

0
18.5 (0.8)
56.2 (0.9)
7.8 (0.8)
63.2 (1.9)
13.0 (0.9)
56.7 (1.2)
13.3 (0.9)
64.4 (1.0)

0.25
19.0 (0.8)
54.3 (0.9)
8.0 (0.6)
59.5 (1.9)
13.2 (0.9)
55.5 (1.0)
12.3 (0.9)
61.6 (1.0)

0.5
23.0 (0.8)
45.4 (1.0)
9.4 (0.8)
50.3 (1.9)
16.0 (1.0)
47.2 (1.5)
16.3 (1.0)
53.8 (1.3)

0.75
30.3 (1.3)
25.9 (0.9)
13.6 (1.1)
31.3 (1.6)
24.6 (1.6)
25.8(1.5)
22.9 (1.4)
29.0 (1.7)

Table 1: Full directional FDR and power (in percentage points) from the simulated data experiment, for various levels of the
feature correlation ρ. (Target FDR level is q = 20%.) Estimated Monte Carlo standard errors are in parentheses.

Least-squares + BH
(after screening)
Knockoffs (with
data splitting)
Knockoffs (with
data recycling)

ρ
FDR
Dir. FDR
FDR
Dir. FDR
FDR
Dir. FDR

0
18.5 (0.8)
18.5 (0.8)
13.0 (0.9)
13.0 (0.9)
13.3 (0.9)
13.3 (0.9)

0.25
17.5 (0.7)
18.0 (0.7)
12.4 (0.8)
12.6 (0.8)
11.6 (0.8)
11.7 (0.8)

0.5
18.5 (0.8)
18.7 (0.8)
11.8 (0.9)
11.9 (0.9)
12.0 (0.9)
12.0 (0.9)

0.75
15.7 (1.0)
15.8 (1.0)
10.2 (1.0)
10.3 (1.0)
7.8 (0.8)
7.8 (0.8)

100

●

80
60

●

●

Least−squares + BH
Selective inf. for sqrt Lasso
Knockoffs + data splitting
Knockoffs + data recycling

●
●

10

15

●

Power (%)

20

●

20

●

0

0

5

Directional FDR (%)

25

●

40

30

Table 2: Reduced FDR and directional FDR (in percentage points) from the simulated data experiment, for various levels of
the feature correlation ρ. (Target FDR level is q = 20%.) Estimated Monte Carlo standard errors are in parentheses.

0.00

0.25

0.50

Feature correlation

0.75

0.00

0.25

0.50

0.75

Feature correlation

Figure 1: Results from the simulated data experiment, plotted against feature correlation ρ. (Target FDR level is q = 20%.)

19

02/2016; Revised 09/2017

6.1

Methods

Data pre-processing. Before applying knockoffs, the data needs to be prepared and our early pre-processing
follows Sabatti et al. [24] and Janson et al. [15]. Genotype features from the original data set were removed if they
met any of the following conditions:
• Not a SNP (some features were, e.g., copy number variations)
• Greater than 5% of values were missing
• All nonmissing values belonged to the same nucleotide
• SNP location could not be aligned to the genome
• A χ2 test rejected Hardy-Weinberg equilibrium at the 0.01% level
• On chromosome 23 (sex chromosome)
The remaining missing values were assumed to take the population frequency of the major allele. In the end, we
have a total of p = 328, 934 SNP features.
For each phenotype, we performed further processing on the response variables. Triglycerides were log-transformed.
C-reactive protein (CRP) was also log-transformed after adding 0.002 mg/l (half the detection limit) to any values
that were recorded as zero. Subjects were excluded from the triglycerides, HDL cholesterol, and LDL cholesterol
if they were on diabetic medication or had not fasted before blood collection (or if either value was missing), or
if they were pregnant, or if their respective phenotype measurement was more than three standard deviations from
the mean, after correcting for sex, oral contraceptive use, and pregnancy. For the four phenotypes of interest to us,
namely, CRP, HDL, LDL and triglycerides, the sample sizes are respectively 5290, 4700, 4682 and 4644. Due to
space limitation, we report our findings on HDL and LDL only.
Finally, in order to correct for population stratification, we regressed the response and the features (SNPs) on the
top five principal components of the design matrix in exactly the same fashion as suggested in Price et al. [23].
Data processing with knockoffs. Now that we have a response/design pair (y, X), we discuss the implementation
of our full GWAS processing pipeline.
1. Data splitting with random rotations. In Section 4.1, we discuss splitting the data (X, y) into two parts,
(X(0) , y(0) ) and (X(1) , y(1) ), with the first part used for screening and the second for inference; in that
section, our discussion centers on splitting by partitioning the n observations into two groups. In practice,
for settings such as GWAS where the design matrix X is naturally somewhat sparse, this may be problematic
as splitting the set of observations can dramatically increase conditioning problems caused by sparsity. Here
we take a different approach to splitting. Under the Gaussian linear model, since the i.i.d. Gaussian noise
term  ∈ Rn has a distribution that is unchanged by rotation, we can rotate the data with a randomly chosen
orthonormal matrix U ∈ Rn×n , to obtain (X0 , y0 ) = (U X, U y). The original linear model is preserved;
as long as the rotation was chosen independently of the data, our model is now
y 0 = X0 β +  0 ,
where 0 = U  ∼ N (0, σ 2 In ). (Note that β is unchanged from the original model.) After rotation, the
new design matrix X0 is now dense, and we are able to split the rotated data into a set of n0 data points (for
screening) and n1 = n − n0 data points (for inference) without matrix conditioning issues. From this point
on, the partitioned data sets (X(0) , y(0) ) and (X(1) , y(1) ) are assumed to be taken from the rotated data,
without further mention. We take n0 = 1900.
2. Feature pre-screening with correlations. The original number of features (SNPs) is 328,934, which is a
bit large to easily run the Lasso for feature screening. As suggested in [35] (where the screen-and-clean
methodology is applied to GWAS), we begin with a pre-screening step that looks only at marginal correlations: using the first part of the data (X(0) , y(0) ), we select the set Sbpre ⊂ [p] of 26,300 features with the
(0) > (0)

largest magnitude correlations Xj

y

.

# times SNP in this region is selected

20

02/2016; Revised 09/2017

LDL
10

Boundary between chromosomes
Appears in meta-analysis
Does not appear in meta-analysis
Appears in Sabatti et al, but not in meta-analysis
Multiple SNPs discovered in this region

*
8
*

*
*

6
4

*
*

2

*

0
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16 17 18 19 20 2122

Location on genome (labeled by chromosome number)

Figure 2: Results of GWAS experiment for the LDL phenotype. The heights of the lines show the number of trials (out of 10
trials) for which at least one SNP in this region was selected; the line type shows whether this same region was identified in
[34] and/or [24]. Regions where more than one SNP was selected (across the 10 trials) are marked with a ‘*’. (See Section 6.2
for details.)
(0)

3. Feature screening with Lasso. We then run the Lasso on (XSb , y(0) ) using the 26,300 pre-screened
pre

features as covariates, and define the screened set of features Sb0 to be the first n1 /4 features, which enter the
Lasso path on this reduced data set.

4. Knockoff filter. Next we bring back the held out part of the data, (X(1) , y(1) ). We apply the knockoff filter
with recycling as defined in Section 4.1, with statistics determined by the sign-restricted Lasso introduced in
Section 4.1.2. Specifically„ we fit the sign-restricted Lasso path
minb∈R2|Sb0 |
s.t.

e b ]bk2 + λkbk1
− [XSb0 X
2
S0
(0)
bj · signj ≥ 0
(0)
bj+|Sb0 | · signj ≥ 0
1
2 ky

e j enters the Lasso path, and set
and set Wj as in (10). (That is, we record the largest λ such that Xj or X
Wj = +λ if Xj enters first, and Wj = −λ otherwise.) Throughout, we use target FDR level q = 0.2.

5. Repeat. We then repeat all these steps, with a new random rotation matrix U ∈ Rn×n , for a total of 10
repetitions. In particular, we are interested in examining whether any selected SNPs appear consistently
across these 10 repetitions of our method.

6.2

Results

Our results for the LDL phenotype are reported in Table 3 and displayed in Figure 2. In total, over the 10 trials (the
random splits of the data), 44 different SNPs are selected. However, as nearby SNPs are extremely correlated, we
consider SNPs whose positions are within 106 base pairs of each other, to be in the same “region”; at this level, in
total over the 10 trials, there are 29 distinct regions discovered. In Table 3 we show, for each region, the selection
frequency of any SNP in that region, i.e. the number of trials out of 10 for which at least one SNP in this region
was selected (we also display the selection frequency for the individual SNPs). We also show any SNPs in the
same region (defined again as a distance of ≤ 106 base pairs) which were identified as associated with LDL in
the meta-analysis of Willer et al. [34]; this meta-analysis works with an extremely large sample size, and therefore
can be viewed as providing a form of “ground truth”. Comparing against the findings of this meta-analysis, we
can try to estimate the FDR of our method as follows: for each trial, we count each selected SNP as a discovery,
and label it as a false discovery if there is no SNP in the meta-analysis from the same region. We then average
this false discovery proportion over the 10 trials, and find an estimated FDR of 32.61%. We also compare against
the findings of Sabatti et al. [24], which analyzes the same data set that we use; this comparison therefore should

21

02/2016; Revised 09/2017

Selection
frequency
(by region)
9

7
7
7
7

6

4
4
3
3
3
3
2
2
2
2
2
2
1
1
1
1
1
1
1
1
1
1
1

Selection
frequency
(by SNP)
7
4
4
2
1
7
7
7
5
2
2
1
1
2
2
1
1

SNP name

Chromosome

rs693
rs3923037
rs754524
rs6754295
rs1429974
rs646776
rs4844614
rs157580
rs11668477
rs11878377
rs1541596
rs10409243
rs688
rs174556
rs174450
rs579383
rs1535

2

4
4
3
3
3
2
1
1
2
2
2
2
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

rs2728487
rs9696070
rs1342165
rs557435
rs4906908
rs579163
rs518181
rs11216267
rs2802955
rs10953541
rs905502
rs1897318
rs10062361
rs6896136
rs12427378
rs2139930
rs7574918
rs6749903
rs3733262
rs981862
rs995124
rs2351643
rs945559
rs1955105
rs2456930
rs4433842
rs4799847

1
1
19
19

11

7
9
1
1
15
11
1
7
8
8
5
12
2
2
4
5
7
8
10
12
15
17
18

Base pair
position
21232195
21148274
21311541
21206183
21300770
109818530
207875175
45395266
11195030
10965270
10987013
10332988
11227602
61580635
61641542
61536583
61597972
47209007
89230779
34800325
55520864
27040082
116713630
116772787
116952392
235015199
107244545
3134810
124864404
74565153
74787310
51074199
51089287
165939179
223653986
129439663
93615466
41764992
63076996
89823147
115344782
62687339
64322683
33686203

Nearby SNPs
found in
Willer et al. [34]
rs1367117

Nearby SNPs
found in
Sabatti et al. [24]
rs693

rs629301

rs646776
rs4844614
rs157580
rs11668477

rs4420638
rs6511720

rs174546

rs174537
rs102275
rs174546
rs174556
rs1535

rs2479409
rs964184
rs514230

rs12916

rs1801689

Table 3: Results of GWAS experiment for the LDL phenotype (see Section 6.2 for details and interpretation of the results).
Physical reference positions for SNPs were drawn from Human Genome Build 37/HG19.

22

02/2016; Revised 09/2017

Selection
frequency
(by region)
9

3

Selection
frequency
(by SNP)
9
9
8
2
2
2
2
3

1
1
1
1
1

1
1
1
1
1

8
4
4

SNP name

Chromosome

rs3764261
rs7499892
rs1532085
rs2575875
rs2740486
rs255049
rs255052
rs7120118

16

11

Base pair
position
56993324
57006590
58683366
107662494
107666513
68013471
68024995
47286290

rs12139970
rs6728178
rs173738
rs2817056
rs2375016

1
2
5
6
12

230406460
21193946
16725880
35734051
109227914

15
9
16

Nearby SNPs
found in
Willer et al. [34]
rs3764261

Nearby SNPs
found in
Sabatti et al. [24]
rs3764261

rs1532085
rs1883025

rs1532085

rs16942887

rs255049

rs3136441

rs2167079
rs7120118

rs4846914

rs7134594

# times SNP in this region is selected

Table 4: Results of GWAS experiment for the HDL phenotype (see Section 6.2 for details and interpretation of the results).
Physical reference positions for SNPs were drawn from Human Genome Build 37/HG19.

HDL
10

Boundary between chromosomes
Appears in meta-analysis
Does not appear in meta-analysis
Appears in Sabatti et al, but not in meta-analysis
Multiple SNPs discovered in this region

8
*

*

6
*

4

*

2
0
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16 17 18 19 20 2122

Location on genome (labeled by chromosome number)

Figure 3: Results of GWAS experiment for the HDL phenotype (same interpretation as in Figure 2).

23

02/2016; Revised 09/2017

Knockoff (data split)
Knockoff (data recy.)
Knockoff (data split)
Knockoff (data recy.)

Target q

dir. FDR

Power

0.1
0.1
0.2
0.2

0.090
0.085
0.180
0.175

0.403
0.470
0.503
0.533

Power (rel. to
screening)
0.632
0.737
0.785
0.833

Table 5: Performance of the knockoff filter averaged over 50 trials. The power is the fraction of true variables selected among
all true positives while the power (rel. to screening) is the fraction of true variables selected among all true positives that have
been screened in.

not be viewed as independent validation, but rather as checking that our analysis agrees with existing work on the
same data set. In Figure 2, we display these findings visually.
For the HDL phenotype, our results are shown in Table 4 and in Figure 3. There are 13 SNPs discovered at
least once in the 10 trials; once grouped into regions, there are 10 distinct regions. Again estimating the FDR
by comparing against Willer et al. [34]’s meta-analysis, we estimate an FDR of 4.29%. This low number of
discoveries, and perhaps conservative FDR, can probably be attributed to fact that this data set has a relatively
small sample size, and fairly weak signal; we would not expect a large number of discoveries in this setting.

6.3

Additional simulations

We conclude this section with a short simulation study assessing directional FDR control and power in a setting
which employs a real genotype matrix as a design but a computer generated response.
• We extract a subset of SNPs with fairly low correlations from the data set processed after performing the
five tasks described at the top of Section 6.1. This is done by applying the PLINK clump8 command using
a significance threshold for index SNPs of 0.2, a physical distance window of 1 kb, a linkage disequilibrium
threshold for clumping of r2 = 0.1, and HDL as a response. This yields a matrix with 26,315 SNPs
with a maximum pairwise correlation equal to 0.316. We proceed by regressing out the top five principal
components of each SNP feature as in Section 6.1, and then normalize the columns. The resulting 4907 ×
26315 design X is fixed throughout the rest of the study.
• For each run, we generate a response y from the model y = Xβ + , where  ∼ N (0, In ) and β has k = 40
nonzero amplitudes equal to ±4.48, each with equal probability, at locations selected uniformly at random.
• We then apply the knockoff procedure, applying either data splitting or data recycling, with random rotations
as earlier. Throughout, we work with n1 = 3000 samples for the inference step, n0 = 1907 samples for
screening step, and a target size for the screened set equal to n1 /4 = 750 SNPs.
The results, averaged over 50 trials, are reported in Table 5. These results confirm our earlier findings: both
knockoff methods control the FDR and the directional FDR, and as expected, we see a gain in power as we move
from data splitting to data recycling.

7

Summary

We have discussed a ‘screen + knockoff’ approach to controlled variable selection, which should appear very
natural to practitioners; in addition, we develop the “recycling” technique to reuse data rather than discarding the
data used for the screening step, thus improving power. A natural question is thus what sort of inference properties
would such a procedure offer? By focusing on Type I + Type S errors, we have shown that screen + knockoff is
guaranteed to have a form of reproducibility in that the directional FDR is rigorously under control. This result
8 For

details about this procedure, see http://pngu.mgh.harvard.edu/~purcell/plink/clump.shtml.

24

02/2016; Revised 09/2017

makes no assumption on the design matrix nor on the values of the regression coefficients, and holds in finite
samples; our contribution is, therefore, very different in spirit from most of the results published in the literature
on high-dimensional inference (the work of Taylor and colleagues referenced earlier being a notable exception),
and employs mathematical arguments centered in martingale theory which are naturally also very different from
the techniques used in the related literature.

Acknowledgements
E. C. was partially supported by the NSF via grant CCF-0963835 and by the Math + X Award from the Simons Foundation.
R. F. B. was partially supported by the NSF via grant DMS-1654076 and by an Alfred P. Sloan fellowship. E. C. would like
to thank Lucas Janson, Chiara Sabatti and Jonathan Taylor for many inspiring discussions about topics ranging from highdimensional inference to GWAS, and useful comments about an early version of the manuscript. We are indebted to Chiara
Sabatti and Christine Peterson for all their help with the NFBC data set and its pre-processing, and to Pragya Sur for her editorial
comments and assistance. E. C. would like to thank Mert Pilanci for his help in running simulations on the Sherlock cluster.

References
[1] Rina Foygel Barber and Emmanuel J. Candès. Controlling the false discovery rate via knockoffs. Ann.
Statist., 43(5):2055–2085, 2015.
[2] Alexandre Belloni, Victor Chernozhukov, and Lie Wang. Square-root lasso: pivotal recovery of sparse signals
via conic programming. Biometrika, 98(4):791–806, 2011.
[3] Alexandre Belloni, Victor Chernozhukov, and Christian Hansen. Inference on treatment effects after selection
among high-dimensional controls. The Review of Economic Studies, 81(2):608–650, 2014.
[4] Y. Benjamini and Y. Hochberg. Controlling the false discovery rate: a practical and powerful approach to
multiple testing. J. Roy. Statist. Soc. Ser. B, 57(1):289–300, 1995.
[5] Yoav Benjamini and Henry Braun. John tukey’s contributions to multiple comparisons. ETS Research Report
Series, 2002(2), 2002.
[6] Yoav Benjamini and Daniel Yekutieli. False discovery rate–adjusted multiple confidence intervals for selected
parameters. Journal of the American Statistical Association, 100(469):71–81, 2005.
[7] Richard Berk, Lawrence Brown, Andreas Buja, Kai Zhang, and Linda Zhao. Valid post-selection inference.
The Annals of Statistics, 41(2):802–837, 2013.
[8] Jianqing Fan and Jinchi Lv. Sure independence screening for ultrahigh dimensional feature space. Journal
of the Royal Statistical Society: Series B (Statistical Methodology), 70(5):849–911, 2008.
[9] William Fithian, Dennis Sun, and Jonathan Taylor. Optimal inference after model selection. arXiv preprint
arXiv:1410.2597, 2014.
[10] Andrew Gelman and Francis Tuerlinckx. Type S error rates for classical and Bayesian single and multiple
comparison procedures. Computational Statistics, 15(3):373–390, 2000.
[11] Max Grazier G’Sell, Trevor Hastie, and Robert Tibshirani. False variable selection rates in regression. arXiv
preprint arXiv:1302.2303, 2013.
[12] Max Grazier G’Sell, Stefan Wager, Alexandra Chouldechova, and Robert Tibshirani. Sequential selection
procedures and false discovery rate control. arXiv preprint arXiv:1309.5352, 2013.
[13] Jian Huang, Shuangge Ma, Cun-Hui Zhang, and Yong Zhou. Semi-penalized inference with direct false
discovery rate control in high-dimensions. arXiv preprint arXiv:1311.7455, 2013.

25

02/2016; Revised 09/2017

[14] Lucas Janson and Weijie Su. Familywise error rate control via knockoffs. arXiv preprint arXiv:1505.06549,
2015.
[15] Lucas Janson, Rina Foygel Barber, and Emmanuel Candès. Eigenprism: Inference for high-dimensional
signal to noise ratios. arXiv preprint arXiv:1505.02097, 2015.
[16] Marjo-Riitta Järvelin, Ulla Sovio, Vanessa King, Liisa Lauren, Baizhuang Xu, Mark I. McCarthy, AnnaLiisa Hartikainen, Jaana Laitinen, Paavo Zitting, Paula Rantakallio, and Paul Elliott. Early life factors and
blood pressure at age 31 years in the 1966 northern finland birth cohort. Hypertension, 44(6):838–846, 2004.
doi: 10.1161/01.HYP.0000148304.33869.ee. URL http://hyper.ahajournals.org/content/
44/6/838.abstract.
[17] Adel Javanmard and Andrea Montanari. Confidence intervals and hypothesis testing for high-dimensional
regression. The Journal of Machine Learning Research, 15(1):2869–2909, 2014.
[18] Lyle V Jones and John W Tukey. A sensible formulation of the significance test. Psychological methods, 5
(4):411, 2000.
[19] Jason D Lee, Dennis L Sun, Yuekai Sun, and Jonathan E Taylor. Exact post-selection inference with the
lasso. arXiv preprint arXiv:1311.6238, 2013.
[20] Hannes Leeb and Benedikt M Pötscher. Can one estimate the conditional distribution of post-model-selection
estimators? The Annals of Statistics, pages 2554–2591, 2006.
[21] Richard Lockhart, Jonathan Taylor, Ryan J. Tibshirani, and Robert Tibshirani. A significance test for the
Lasso. The Annals of Statistics, 42(2):413–468, 2014.
[22] Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad. Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition. In Proceedings of 27th Asilomar Conference on Signals,
Systems and Computers, pages 40–44. IEEE, 1993.
[23] Alkes L Price, Nick J Patterson, Robert M Plenge, Michael E Weinblatt, Nancy A Shadick, and David Reich.
Principal components analysis corrects for stratification in genome-wide association studies. Nature genetics,
38(8):904–909, 2006.
[24] Chiara Sabatti, Susan K Service, Anna-Liisa Hartikainen, Anneli Pouta, Samuli Ripatti, Jae Brodsky, Chris G
Jones, Noah A Zaitlen, Teppo Varilo, Marika Kaakinen, Ulla Sovio, Aimo Ruokonen, Jaana Laitinen, Eveliina Jakkula, Lachlan Coin, Clive Hoggart, Andrew Collins, Hannu Turunen, Stacey Gabriel, Paul Elliot,
Mark I McCarthy, Mark J Daly, Marjo-Riitta Järvelin, Nelson B Freimer, and Leena Peltonen. Genome-wide
association analysis of metabolic traits in a birth cohort from a founder population. Nature genetics, 41(1):
35–46, 2009.
[25] Juliet Popper Shaffer. Multiplicity, directional (type iii) errors, and the null hypothesis. Psychological methods, 7(3):356, 2002.
[26] Weijie Su, Malgorzata Bogdan, and Emmanuel Candès. False discoveries occur early on the lasso path. arXiv
preprint arXiv:1511.01957, 2015. Technical Report, Department of Statistics, Stanford University.
[27] J.
Taylor.
selective-inference.
selective-inference, 2017.

https://github.com/jonathan-taylor/

[28] Xiaoying Tian, Joshua R Loftus, and Jonathan E Taylor. Selective inference with unknown variance via the
square-root lasso. arXiv preprint arXiv:1504.08031, 2015.
[29] R. Tibshirani. Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society.
Series B (Methodological), pages 267–288, 1996.

26

02/2016; Revised 09/2017

[30] Ryan J Tibshirani, Jonathan Taylor, Richard Lockhart, and Robert Tibshirani. Exact post-selection inference
for sequential regression procedures. Journal of the American Statistical Association, 111(514):600–620,
2016.
[31] John W Tukey. The philosophy of multiple comparisons. Statistical science, pages 100–116, 1991.
[32] Arend Voorman, Ali Shojaie, and Daniela Witten. Inference in high dimensions with the penalized score test.
arXiv preprint arXiv:1401.2678, 2014.
[33] Larry Wasserman and Kathryn Roeder. High dimensional variable selection. Annals of statistics, 37(5A):
2178, 2009.
[34] Cristen J. Willer, Ellen M. Schmidt, Sebanti Sengupta, et al. Discovery and refinement of loci associated with
lipid levels. Nature genetics, 45(11):1274–1283, 2013.
[35] Jing Wu, Bernie Devlin, Steven Ringquist, Massimo Trucco, and Kathryn Roeder. Screen and clean: a tool
for identifying interactions in genome-wide association studies. Genetic epidemiology, 34(3):275–285, 2010.
[36] Cun-Hui Zhang and Jian Huang. The sparsity and bias of the lasso selection in high-dimensional linear
regression. The Annals of Statistics, pages 1567–1594, 2008.
[37] Cun-Hui Zhang and Stephanie S. Zhang. Confidence intervals for low dimensional parameters in high dimensional linear models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):
217–242, 2014.
[38] Peng Zhao and Bin Yu. On model selection consistency of lasso. The Journal of Machine Learning Research,
7:2541–2563, 2006.

A
A.1

Proofs
A general form of directional FDR control

d j is
Our results, namely, Theorems 1–3, are all special cases of the general statement below. Throughout, sign


e j )> y . Also, we shall say that a positive semidefinite matrix
d j = sign (Xj − X
defined as in (16), i.e. sign

M ∈ R2p×2p satisfies the pairwise exchangeability condition if Mj,k = Mj,k+p = Mj+p,k = Mj+p,k+p for all
j 6= k ∈ {1, . . . , p} and Mj,j = Mj+p,j+p for all j ∈ {1, . . . , p}.

Theorem 4. Suppose y ∼ N (µ, Θ) for some mean vector µ ∈ Rn and some covariance Θ ∈ Rn×n . Let




e ∈ Rn×p be any fixed matrices such that X X
e > ·Θ· X X
e obeys the pairwise exchangeability condition.
X, X
Let W = (W
the sufficiency and antisymmetry properties (8) and (9). Finally,
n 1 , . . . , Wp ) be any statistic satisfying
o
>
e
d
let Vj = 1 signj 6= sign((Xj − Xj ) µ be the indicator of a sign error on feature j. Then the knockoff+ method
controls a directional FDR defined as

FDRdir

o 
 n
j ∈ Sb and Vj = 1
≤q.
= E
b
1 ∨ |S|

If the knockoff method is used instead, then a modified directional FDR is controlled:
o 
 n
j ∈ Sb and Vj = 1
≤q.
mFDRdir = E 
b + q −1
|S|

We prove this theorem in Appendix A.3. The rest of this section shows how all of our results follow as special
cases.

27

02/2016; Revised 09/2017

Proof of Theorem 1 This is straightforward since the result is obtained by taking µ = Xβ and Θ = σ 2 In in
e j )> µ) = sign(βj ) for all j such that Xj 6= X
e j . (We ignore
the distribution of y, and noting that sign((Xj − X
e j , for in this case we would get Wj = 0 and thus feature j could
the trivial case where sj = 0 implying Xj = X
never be selected.)
The more general form of Theorem 4 also allows us to consider other settings, such as those involving variable
screening in high dimensions. For instance, our next proof concerns the testing of regression coefficients in the full
model under the sure screening property.
Proof of Theorem 2 We prove directional FDR control by conditioning on y(0) . For any y(0) such that the event
E holds, we will show that for the knockoff procedure,
i
h
(25)
E mFDRdir y(0) · 1E ≤ q ,
and for the knockoff+,

i
h
E FDRdir y(0) · 1E ≤ q .

(26)

From this point on we treat y(0) as fixed. The conditional distribution of the response y is given by
 (0)  


y
0
0
(0)
y|y
∼N
,
.
0 σ 2 I n1
X(1) β
If y(0) is such that the event E occurs, then we can simplify this to
#
"


y(0)
0
(0)
and Θ =
y|y
∼ N (µ, Θ), where µ =
(1)
0
XSb βSb0

0
σ 2 I n1

0

(1)

(1)



(0)

(1)

(27)

(28)

,
(0)

e
since βj = 0 for all j 6∈ Sb0 when E occurs. Define ΣSb = XSb > XSb . Recall that X
b0 by (18), and that
b0 = XS
S
0
0
0
(1)
e satisfies
X
b0
S
#
"
(1)
(1)
ΣSb
ΣSb − diag{s}
 (1) (1) >  (1) (1) 
e
e
0
0
(29)
XSb X
XSb X
(1)
(1)
b0 =
b0
S
S
0
0
ΣSb
ΣSb − diag{s}
0

0

for some vector s ≥ 0, by the knockoff construction. We then see that




 



e (1) > X(1) X
e (1) ,
e b >Θ X b X
e b = σ 2 X(1) X
XSb0 X
b
b
b
b
S0
S0
S0
S
S
S
S
0

0

0

0

is a pairwise exchangeable matrix. At this point, we have satisfied the conditions of Theorem 4. Now what is a
sign error in Theorem 4? We calculate that for any j ∈ Sb0 ,
X (1)
e j )> µ = (X(1) − X
e (1) )> X(1) β b =
e (1) )> X(1) · βk = sj · βj ,
(Xj − X
(Xj − X
j
j
j
k
b
S0
S
0

b0
k∈S

where the first equality holds on the event E (i.e. βk = 0 for k 6∈ Sb0 ), and the third equality uses the Gram matrix
e j , an error sign
e j )> µ) is the same as
d j 6= sign((Xj − X
condition (29). Hence, since sj > 0 whenever Xj 6= X
d j 6= sign(βj ), which completes the proof.
sign
We now move on to our result concerning regression coefficients from the reduced model.

28

02/2016; Revised 09/2017

Proof of Theorem 3 This last proof is more subtle and for this, we treat X as random in addition to y, but will
(1)
e (1) as well as on y(0) . By the assumption
condition on X(0) and on the screened features XSb and knockoffs X
b0
S
0
that the rows of X are i.i.d. draws from N (0, Ψ), we can write
(1)

(1)

0

0

XSbc = XSb · Γ + G · Λ /2 ,
1

where Γ and Λ are the fixed unknown matrices given by

Γ = (ΨSb0 ,Sb0 )−1 · ΨSb0 ,Sbc
0

and the Schur complement

Λ = ΨSbc ,Sbc − ΨSbc ,Sb0 · (ΨSb0 ,Sb0 )−1 · ΨSb0 ,Sbc ,
0

b0 |)
n1 ×(p−|S

and where G ∈ R
(1) e (1) (1)
(XSb , X
). Then
b ,
S
0

0

0

0

has i.i.d. standard normal entries drawn independently from (X(0) , y(0) ) and from

0



1
(1)
y(1) = X(1) β + (1) = XSb · βSb0 + ΓβSbc + G · Λ /2 βSbc + (1) ,

and note that

0

0

0

βSb0 + ΓβSbc = β partial
0

by definition of the expected partial regression coefficients. This is because from (21), we get that

h
i


(1)
(1)
(1)
(1)
(1)
(1)
(1)
β partial = XSb † XSb βSb0 + E XSbc βSbc XSb ; X(0) , y(0) = XSb † XSb βSb0 + XSb ΓβSbc
0

0

0

0

0

0

0

0

0

= βSb0 + ΓβSbc ,
0

(1) †

where XSb

0

(1) >

is a short hand for (XSb

0

y|

(1)

(1) >

XSb )−1 XSb

(1) e (1) 
X(0) , y(0) , XSb , X
b0
S
0

0

0

∼N

"

. Therefore,

# "
0
y(0)
(1) partial ,
0
XSb β
0



0



2
βS>
I n1
bc + σ
bc ΛβS
0

0

#!

.

We are now in the same setting as in the proof of Theorem 2 (see (27)), except with β partial in place of β and with
a new variance level that has increased due to the randomness in the missed signals. The remainder of the proof
thus proceeds identically.

A.2

Key lemmas

The proof of Theorem 4 relies on two key lemmas: Lemma 1 earlier which treats non-i.i.d. Bernoulli sequences,
and the following result regarding the statistics W arising in the knockoff procedure.
Lemma 2. With the assumptions from Theorem 4, let V be the σ-algebra generated by the random variables
e > y, |(X − X)
e > y|}, and define the vector S = (S1 , . . . , Sp ) of signs as
{(X + X)
Then the following holds:


e j )> y · sign(Wj ) .
Sj = sign (Xj − X

• |W| and S belong to V.
• The signs of Wj are mutually independent after conditioning on V.

e j )> µ , then
• For each j = 1, . . . , p, if Sj 6= sign (Xj − X

P {sign(Wj ) = −1 | V} ≥ 1/2.

29

02/2016; Revised 09/2017

Before presenting the formal proof, we give an intuitive explanation for why the probabilities for sign(Wj ) work
in our favor; that is, why we are able to obtain P {sign(Wj ) = −1 | V} ≥ 1/2, which ensures that we are more
likely to over- rather than under-estimate the number of false positives. This is similar to the discussion for the
low-dimensional setting as in (16). In this case, for the high-dimensional setting, we see that if Sj 6= sign (Xj −

e j )> µ , then:
X
e j )> µ ≥ 0, then sign(Wj ) = −1 whenever (Xj − X
e j )> y > 0;
• If (Xj − X
e j )> µ ≤ 0, then sign(Wj ) = −1 whenever (Xj − X
e j )> y < 0.
• If (Xj − X

e j )> y is a normal random variable whose mean is given by (Xj − X
e j )> µ. If this mean is
Of course, (Xj − X
>
e
nonnegative, then the random variable (Xj − Xj ) y is more likely to be positive than negative; if this mean is
nonpositive, then the opposite is true. Either way, we see that sign(Wj ) = −1 has probability at least 1/2.


 



e > X X
e and X X
e > y or, equivalently, is
Proof. First, by the sufficiency property, W is a function of X X
of the form



 
e > X X
e , (X + X)
e > y, |(X − X)
e > y|, sign((X − X)
e > y)
W=w X X

for some function w. Second, the antisymmetry property implies that flipping the signs of any subset of the
last argument of the function w flips the signs of the corresponding outputs. (Of course, for any j such that
e j )> y = 0, the antisymmetry property implies that Wj = 0; we ignore these features in our proof from
(Xj − X
this point on as they can never be selected by our method.) This property makes clear that V determines W up to
a sign change, so that |W| ∈ V. Furthermore, define a vector of signs S such that


 

e > X X
e , (X + X)
e > y, |(X − X)
e > y|, S .
|W| = w X X
Clearly, S ∈ V and, by the antisymmetry property, we also see that


 

e > X X
e , (X + X)
e > y, |(X − X)
e > y|, S · sign(W) = W.
w X X

e > y) or, expressed differently, S = sign(W) ·
This equality means that S · sign(W) must be equal to sign((X − X)
>
e y) thereby coinciding with the definition given in the statement of the lemma. This establishes the
sign((X − X)
first claim.
Next, an elementary calculation shows that
"
#
"
# 
e >y
e >µ
e > Θ(X + X)
e
(X + X)
(X + X)
(X + X)
∼
N
,
>
>
e y
e µ
0
(X − X)
(X − X)

!
0
,
diag{d}

(30)

where d = (d1 , . . . , dp ) ∈ Rp is some nonnegative vector. The form of the covariance term arises from the fact
e > Θ(X + X)
e = 0 and (X − X)
e > Θ(X − X)
e is diagonal, according to the pairwise exchangeability
that (X − X)

> 

e
e
e j )> y are mutually independent, and are inassumption on X X Θ X X . Therefore, the terms (Xj − X
>
e y. Also, since S ∈ V the terms sign(Wj ) = Sj · sign((Xj − X
e j )> y) are mutually
dependent from (X + X)
independent after conditioning on V, proving our second claim.

Finally, observe that if Z ∼ N (θ, σ 2 ), then it follows from a simple calculation that


P {sign(Z) = −1 | |Z|}
2θ|z|
= exp − 2
.
P {sign(Z) = +1 | |Z|}
σ
n
o
e j )> y) = ±1 V , applying the above formula gives
Now since P {sign(Wj ) = ±1 | V} = P Sj · sign((Xj − X
e > µ · |(Xj − X
e j )> y)|
P {sign(Wj ) = −1 | V}
2Sj · (Xj − X)
= exp −
P {sign(Wj ) = +1 | V}
dj

!

.

30

02/2016; Revised 09/2017


e j )> µ , then the argument in the exponential in the right-hand side of this last
Clearly, if Sj 6= sign (Xj − X
inequality is positive and thus, P {sign(Wj ) = −1 | V} ≥ 1/2.

The original knockoff paper [1] worked in the simpler setting where the signs of Wj for the nulls are i.i.d. unbiased.
If H0 is the set of nulls and T > 0 the adaptive threshold of the knockoff filter, [1] established that the i.i.d. property
for the signs gives

 
|{j ∈ H0 : Wj ≥ T }|
|W| ≤ 1.
E E
1 + |{j ∈ H0 : Wj ≤ −T }|
Here, we use Lemma 1 to develop a similar bound when the signs are no longer i.i.d. unbiased but are, instead,
independent with at most a 50% chance of being positive conditionally on V.

Corollary 1. Suppose M ⊂ {1, . . . , p} is a random set belonging to V and chosen such that for each j ∈ M,
|Wj | > 0 and P {sign(Wj ) = −1 | V} ≥ ρ. Letting T > 0 be the adaptive threshold of either the knockoff or
knockoff+ filter, then


|{j ∈ M : Wj ≥ T }|
E
V ≤ ρ−1 − 1.
(31)
1 + |{j ∈ M : Wj ≤ −T }|

Proof. The proof is straightforward and only consists in rewriting the left-hand side of (31) in such a way that we
recognize the formulation from Lemma 1. To begin with, we can treat the ordering of |W| as fixed since Lemma 2
gives |W| ∈ V. Now assume without loss of generality that M = {1, . . . , m} and reorder the indices of the Wj ’s
in M so that |W(1) | ≥ |W(2) | ≥ · · · ≥ |W(m) |. Set Bj = 1W(j) <0 . Next we will condition also on sign(Wj ) for
e be the larger σ-algebra generated by V and {sign(Wj ) : j 6∈ M}. By definition of M, when we
all j 6∈ M. Let V
⊥
⊥

condition on V, we have Bj ∼ Bernoulli(ρj ) with ρj ≥ ρ for all j ∈ M. By Lemma 2, the same is still true even
after conditioning on {sign(Wj ) : j 6∈ M} also, since the signs of W are mutually independent conditional on V.
To summarize, we have
⊥
⊥
e∼
(Bj )j∈M | V
Bernoulli(ρj ),
with ρj ≥ ρ for all j ∈ M.
Next, we can write

#{j ∈ M : Wj ≥ T }
(1 − B1 ) + · · · + (1 − BJ )
1+J
=
=
−1,
1 + #{j ∈ M : Wj ≤ −T }
1 + B1 + · · · + BJ
1 + B1 + · · · + BJ

where J is the index such that |W(1) | ≥ · · · ≥ |W(J) | ≥ T > |W(J+1) | ≥ · · · ≥ |W(m) |. After conditioning on
e so that we can treat |W| and {sign(Wj ) : j 6∈ M} as fixed, we observe that the index J can be expressed as a
V,
stopping time, in reverse time, with respect to the filtration {Fj } given by
Fj = {B1 + · · · + Bj , Bj+1 , . . . , Bm }.

The conclusion now follows from Lemma 1.

A.3

Proof of Theorem 4

First, consider control of mFDRdir under the knockoff. Since j ∈ Sb if and only if Wj ≥ T where T is our adaptive
threshold, the modified directional FDR is given by
o 
 n
e j )> µ)
d j 6= sign((Xj − X
j : Wj ≥ T, sign

mFDRdir = E 
|{j : Wj ≥ T }| + q −1
o
 n

e j )> µ)
d j 6= sign((Xj − X
j : Wj ≥ T, sign
1
+
|{j
:
W
≤
−T
}|
j

= E
·
1 + |{j : Wj ≤ −T }|
|{j : Wj ≥ T }| + q −1
o 
 n
e j )> µ)
d j 6= sign((Xj − X
j : Wj ≥ T, sign

≤ q · E
1 + |{j : Wj ≤ −T }|

31

02/2016; Revised 09/2017

since |{j : Wj ≤ −T }| ≤ q · |{j : Wj ≥ T }| by definition of T . Now rewrite the set of sign errors as
n
o n
o

e j )> µ) = j : Wj ≥ T, sign (Xj − X
e j )> y 6= sign((Xj − X
e j )> µ)
d j 6= sign((Xj − X
j : Wj ≥ T, sign
n
o
e j )> µ ,
= j : Wj ≥ T, Sj 6= sign (Xj − X

d j in (16), and the second from the definition Sj =
where the first equality follows from the definition of sign

>
e
c0 = {j :
sign (Xj − Xj ) y · sign(Wj ) as in Lemma 2 (since sign(Wj ) = +1 for all j in this set). Define H

e j )> µ } and let V be as in Lemma 2. By the tower law we obtain
Sj 6= sign (Xj − X

  n
o
e j )> µ
j : Wj ≥ T, Sj 6= sign (Xj − X
V 
mFDRdir ≤ q · E E 
1 + |{j : Wj ≤ −T }|
o

  n
c0 : Wj ≥ T
j∈H
= q · E E 
V 
1 + |{j : Wj ≤ −T }|
n
o
 

c0 : Wj ≥ T
j∈H
n
o V  .
≤ q · E E 
c0 : Wj ≤ −T
1+ j ∈H

c0 ∈ V since S ∈ V by Lemma 2, and so inside of the conditional expectation, we can treat the set
Next, note that H
c
H0 as fixed. Lemma 1 proves that conditional on V, the signs of W are independent, with P {sign(Wj ) = −1 | V} ≥
c0 . Therefore, applying Corollary 1 with ρ = 1/2 gives
1/2 for all j ∈ H
n
o


c0 : Wj ≥ T
j∈H
n
o V  ≤ 1,
E
c0 : Wj ≤ −T
1+ j ∈H
which concludes the argument.

The proof for FDRdir , when using knockoff+ instead of knockoff, follows similarly: we have
o

 n
e j )> µ)
d j 6= sign((Xj − X
j : Wj ≥ T+ , sign
1 + |{j : Wj ≤ −T+ }| 
·
FDRdir = E 
1 + |{j : Wj ≤ −T+ }|
1 ∨ |{j : Wj ≥ T+ }|
n
o


e j )> µ)
d j 6= sign((Xj − X
j : Wj ≥ T+ , sign
 ,
≤ q · E
1 + |{j : Wj ≤ −T+ }|

where the last step uses the (slightly more conservative) definition of the threshold T+ for the knockoff+ procedure.
From then on, the argument is the same as before.

A.4

Proof of Lemma 1

We first give a slight generalization of Lemma 4 in [1]:
iid

Lemma 3. Suppose that B1 , . . . , Bn ∼ Bernoulli(ρ). Let J be a stopping time in reverse time with respect to the
filtration {Fj }, where Fj 3 B1 + · · · + Bj , Bj+1 , . . . , Bn , and where the variables B1 , . . . , Bj are exchangeable
with respect to Fj . Then


1+J
E
≤ ρ−1 .
1 + B1 + · · · + BJ

32

02/2016; Revised 09/2017

Proof of Lemma 3. Define the sum
Sj = B1 + · · · + Bj ∈ Fj ,

and define the process
Mj =

1+j
1+j
=
∈ Fj .
1 + B1 + · · · + Bj
1 + Sj

In [1] it is shown that E [Mn ] ≤ ρ−1 , therefore, by the optional stopping time theorem it suffices to show that
{Mj } is a supermartingale with respect to {Fj }. First, since {B1 , . . . , Bj+1 } are exchangeable with respect to
Fj+1 , we have
Sj+1
.
P {Bj+1 = 1 | Fj+1 } =
1+j
Therefore, if Sj+1 > 0,
1+j
1+j
· P {Bj+1 = 0 | Fj+1 , Sj+1 } +
· P {Bj+1 = 1 | Fj+1 , Sj+1 }
1 + Sj+1
1 + Sj+1 − 1
1 + j − Sj+1
1+j
Sj+1
1+j
·
+
·
=
1 + Sj+1
1+j
1 + Sj+1 − 1 1 + j
1 + (j + 1)
1 + j − Sj+1
+1=
= Mj+1 .
=
1 + Sj+1
1 + Sj+1

E [Mj | Fj+1 ] =

If instead Sj+1 = 0, then trivially Sj = 0 also, and so Mj = 1 + j < 2 + j = Mj+1 . This proves that {Mj } is a
supermartingale with respect to {Fj }, as desired.
iid

Corollary 2. Suppose that A ⊆ [n] is fixed, while B1 , .P
. . , Bn ∼ Bernoulli(ρ). Let J be a stopping time in reverse
time with respect to the filtration {Fj }, where Fj 3 i≤j,i∈A Bi , and the variables {Bi : i ≤ j, i ∈ A} are
exchangeable with respect to Fj . Then
"
#
1 + |{i ≤ J : i ∈ A}|
P
E
≤ ρ−1 .
1 + i≤J,i∈A Bi
Proof of Corollary 2. Let A = {i1 , . . . , im } where 1 ≤ i1 < · · · < im ≤ n. Then by considering the i.i.d.
sequence
B i1 , . . . , B im
in place of B1 , . . . , Bn , we see that this result is equivalent to Lemma 3.
Proof of Lemma 1. We may assume ρ < 1 to avoid the trivial case. We now create a different construction for the
Bi ’s. First, generate a random set A ⊆ [n] where for each i, independently,
P {i ∈ A} =

1 − ρi
.
1−ρ

Next, define variables
iid

Q1 , . . . , Qn ∼ Bernoulli(ρ),

which are generated independently of the random set A. Finally, define
Bi = Qi · 1i∈A + 1i6∈A .

(32)

Then, clearly, the Bi ’s are mutually independent with P {Bi = 1} = ρi , as required by the lemma. Next, since
Bi = Qi · 1i∈A + 1i6∈A for all i, we have
1+J
1 + |{i ≤ J : i ∈ A}| + |{i ≤ J : i 6∈ A}|
1 + |{i ≤ J : i ∈ A}|
P
P
=
≤
,
1 + B1 + · · · + BJ
1 + i≤J,i∈A Qi + |{i ≤ J : i 6∈ A}|
1 + i≤J,i∈A Qi

33

02/2016; Revised 09/2017

a
where the last step uses the identity a+c
b+c ≤ b whenever 0 < b ≤ a and c ≥ 0. Therefore, it will be sufficient to
bound the right-hand side. We will use Corollary 2 to prove that
"
#
1 + |{i ≤ J : i ∈ A}|
P
A ≤ ρ−1 ,
(33)
E
1 + i≤J,i∈A Qi

which will be sufficient to prove the lemma by applying the tower law of expectations.
e i = Qi · 1i∈A , and define a filtration {F 0 } where F 0 is the σ-algebra generated as
To prove (33), first let Q
j

j

n
o
0
e
e
e
e
Fj = σ Q1 + · · · + Qj , Qj+1 , . . . , Qn , A
.

Next, for any j, by (32) we see that

B1 + · · · + Bj , Bj+1 , . . . , Bn ∈ Fj0

⇒

Fj0 ⊇ Fj ,

and so T is a stopping time (in reverse time) with respect to {Fj0 } also. We also see that trivially, the variables
{Qi : i ≤ j, i ∈ A} are exchangeable with respect to Fj0 (since the Qi ’s are i.i.d. and are independent from A).
Finally, since the Qi ’s are independent from A, the desired bound (33) follows directly from Corollary 2 after
conditioning on A.

