Pan for gold

## Abstract

Abstract Many contemporary large-scale applications involve building
interpretable models linking a large set of potential covariates to
a response in a nonlinear fashion, such as when the response is
binary. Although this modeling problem has been extensively studied,
it remains unclear how to effectively control the fraction of false
discoveries even in high-dimensional logistic regression, not to
mention general high-dimensional nonlinear models. To address such a
practical problem, we propose a new framework of model- X knockoffs, which reads from a different
perspective the knockoff procedure (Barber and Candès,, 2015 ) originally
designed for controlling the false discovery rate in linear
models. Whereas the knockoffs procedure is constrained to
homoscedastic linear models with n ≥ p 𝑛 𝑝 n\geq p , the key innovation here
is that model-X knockoffs provide valid inference from
finite samples in settings in which the conditional distribution
of the response is arbitrary and completely unknown. Furthermore,
this holds no matter the number of covariates. Correct inference
in such a broad setting is achieved by constructing knockoff
variables probabilistically instead of
geometrically. To do this, our approach requires the covariates be
random (independent and identically distributed rows) with a
distribution that is known, although we provide preliminary experimental evidence that our procedure is robust to unknown/estimated distributions. To our knowledge, no
other procedure solves the controlled variable selection
problem in such generality, but in the restricted settings where
competitors exist, we demonstrate the superior power of knockoffs
through simulations. Finally, we apply our procedure to data from
a case-control study of Crohn’s disease in the United Kingdom,
making twice as many discoveries as the original analysis of the
same data. Keywords. False discovery rate (FDR), knockoff
filter, testing for conditional independence in nonlinear models,
generalized linear models, logistic regression, Markov blanket,
genome-wide association study (GWAS)

## 1 Introduction

### 1.1 Panning for gold

Certain diseases have a genetic basis, and an important biological
problem is to find which genetic features (e.g., gene expressions or
single nucleotide polymorphisms) are important for determining a given
disease. In health care, researchers often want to know which
electronic medical record entries determine future medical
costs. Political scientists study which demographic or socioeconomic
variables determine political opinions. Economists are similarly
interested in which demographic/socioeconomic variables affect future
income. Those in the technology industry seek out specific software
characteristics they can change to increase user engagement. In the
current data-driven science and engineering era, a list of such
problems would go on and on. The common theme in all these instances
is that we have a deluge of explanatory variables, often many more
than the number of observations, knowing full well that the outcome we
wish to understand better only actually depends on a small fraction of
them. Therefore, a primary goal in modern “big data analysis” is to
identify those important predictors in a sea of noise
variables. Having said this, a reasonable question is why do we have
so many covariates in the first place? The answer is twofold: first,
because we can. To be sure, it may be fairly easy to measure
thousands if not millions of attributes at the same time. For
instance, it has become relatively inexpensive to genotype an
individual, collecting hundreds of thousands of genetic variations at
once. Second, even though we may believe that a trait
or phenotype depends on a comparably small set of genetic
variations, we have a priori no idea about which ones are
relevant and therefore must include them all in our search for
those nuggets of gold, so to speak. To further complicate matters, a
common challenge in these big data problems, and a central focus of
this paper, is that we often have little to no knowledge of how the
outcome even depends on the few truly important variables.

To cast the ubiquitous *(variable) selection* problem in statistical terms, call $Y$ the random variable representing the
outcome whose determining factors we are interested in, and
$X_{1},\dots,X_{p}$ the set of $p$ potential determining factors or explanatory variables. The
object of study is the *conditional* distribution of the outcome
$Y$ given the covariates $X=(X_{1},\dots,X_{p})$, and we shall denote
this conditional distribution function by $F_{Y|X}$. Ideally we would
like to estimate $F_{Y|X}$, but in general, this is effectively impossible
from a finite sample. For instance, even knowing that the conditional
density depends upon 20 known covariates makes the problem
impossible unless either the sample size $n$ is astronomically large,
and/or we are willing to impose a very restrictive model.
However, in most problems $F_{Y|X}$ may realistically be assumed to
depend on a small fraction of the $p$ covariates; that is, the
function $F_{Y|X}(y|x_{1},\ldots,x_{p})$ only depends upon a small number
of coordinates $x_{i}$ (or is well approximated by such a
lower-dimensional function).
Although this assumption does not magically make the estimation of
$F_{Y|X}$ easy, it does suggest consideration of the simpler problem:
which of the many variables does $Y$ depend upon? Often, finding
a few of the important covariates—in other words, teasing out the
relevant factors from those which are not—is already scientifically
extremely useful and can be considered a first step in understanding
the dependence between an outcome and some interesting variables; we
regard this as a crucial problem in modern data science.

### 1.2 A peek at our contribution

This paper addresses the selection problem by considering a very
general conditional model, where the response $Y$ can depend in an
arbitrary fashion on the covariates $X_{1},\dots,X_{p}$. The only
restriction we place on the model is that the observations
$(X_{i1},\dots,X_{ip},Y_{i})$ are independently and identically
distributed (i.i.d.), which is often realistic in high-dimensional
applications such as genetics, where subjects may be drawn randomly
from some large population, or client behavioral modeling, where
experiments on a service or user interface go out to a random subset
of users. Therefore, the model is simply

$$ $(X_{i1},\dots,X_{ip},Y_{i})\stackrel{{\scriptstyle\emph{i.i.d.}}}{{\sim}}F_{XY},\;\;\;i=1,\dots,n$ (1.1) $$

for some arbitrary $(p+1)$-dimensional joint distribution $F_{XY}$. We
will assume *no knowledge* of the conditional distribution of
$Y\,|\,X_{1},\dots,X_{p}$, but we do assume the joint distribution of the
covariates is known, and we will denote it by $F_{X}$. As a
concrete example, consider a case-control experiment to determine
the genetic factors which contribute to a rare disease, with
diseased subjects oversampled to 50% to increase power. Then the
joint distribution of features and disease status obeys the model
([1.1](#S1.E1)), where $F_{X}$ is a 50%–50% mixture of the genetic
distributions of diseased and healthy subjects, and $Y$ is the
subjects’ binary disease state.

In Section [1.3](#S1.SS3) below we shall discuss the merits of this
model but we would like to immediately remark on an important benefit:
namely, one can pose a meaningful problem. To do this, observe that
when we say that the conditional distribution of $Y$ actually depends
upon a (small) subset $\mathcal{S}\subset\{1,\ldots,p\}$ of the
variables $X_{1},\ldots,X_{p}$, which we would like to identify, we mean
that we would like to find the “smallest” subset $\mathcal{S}$ such
that conditionally on $\{X_{j}\}_{j\in\mathcal{S}}$, $Y$ is
independent of all other variables. Another way to say this is that
the other variables do not provide additional information about $Y$. A
minimal set $\mathcal{S}$ with this property is usually called a
Markov blanket or Markov boundary for $Y$ in the literature on
graphical models (Pearl,, , Section 3.2.1). Under very mild
conditions about the joint distribution $F_{XY}$, the Markov blanket
is well defined and unique (see Section [2](#S2)
for details) so that we have a cleanly stated selection
problem. Note also that the Markov Blanket can be defined purely
in terms of $F_{Y|X}$ without any reference to $F_{X}$, so that in our
case-control example the problem is defined in exactly the same way
as if $F_{X}$ were the true population genetic distribution instead of
the oversampled mixture of diseased and healthy genetic
distributions.

In most problems of interest, even with the knowledge of $F_{X}$, it is
beyond hope to recover the blanket $\mathcal{S}$ with no error. Hence,
we are naturally interested in procedures that control a Type I error;
that is, we would like to find as many variables as possible while at
the same time not having too many false positives. In this paper, we
focus on controlling the false discovery rate (FDR)
(Benjamini and Hochberg,, ), which we can define here as follows: letting
$\hat{\mathcal{S}}$ be the outcome of a selection procedure operating
on the sampled data (we put a hat because $\hat{\mathcal{S}}$ is
random), the FDR is

$$ $\text{FDR}:=\mathbb{E}[\text{FDP}],\qquad\text{FDP}=\frac{\#\{j:j\in\hat{\mathcal{S}}\setminus\mathcal{S}\}}{\#\{j:j\in\hat{\mathcal{S}}\}}$ (1.2) $$

with the convention 0/0 = 0. Procedures that control the FDR are
interpretable, as they roughly bound what fraction of discoveries are
false ones, and they can be quite powerful as well.

One achievement of this paper is to show that we can design rather
powerful procedures that rigorously control the FDR ([1.2](#S1.E2)) in
finite samples. This holds no matter the unknown relationship between
the explanatory variables $X$ and the outcome $Y$. We achieve this by
rethinking the conceptual framework of Barber and Candès, (), who
originally introduced the knockoff procedure (throughout this paper,
we will sometimes use “knockoffs” as shorthand for the knockoff
framework or procedure). Their salient idea was to construct a set of
so-called “knockoff” variables which were not (conditionally on the
original variables) associated with the response, but whose structure
mirrored that of the original covariates. These knockoff variables
could then be used as controls for the real covariates, so that only
real covariates which appeared to be considerably more associated with
the response than their knockoff counterparts were selected. Their
main result was achieving exact finite-sample FDR control in the
homoscedastic Gaussian linear model when $n\geq 2p$ (along with a
nearly-exact extension to when $p\leq n<2p$). By reading the knockoffs
framework from a new perspective, the present paper places *no
restriction* on $p$ relative to $n$, in sharp constrast to the
original knockoffs work which required the low-dimensional setting of
$n\geq p$. The conceptual difference is that the original knockoff
procedure treats the $X_{ij}$ as fixed and relies on specific
stochastic properties of the linear model, precluding consideration of
$p>n$ or nonlinear models. By treating the $X_{ij}$ as *random*
and relying on that stochasticity instead, the model-X
perspective allows treatment of the high-dimensional setting which is
increasingly the norm in modern applications. We refer to the new
approach as *model-X* (MX) knockoffs, and by
contrast we refer to the original knockoffs approach of
Barber and Candès, () as *fixed-X* (FX) knockoffs. In a nutshell:

- •
We propose a new knockoff construction amenable to the random
covariate setting ([1.1](#S1.E1)).
- •
As in Barber and Candès, () and further reviewed in
Section [3](#S3), we shall use the knockoff variables
as controls in such a way that we can tease apart important
variables from noise while controlling the FDR. We place no
restriction on the dimensionality of the data or the conditional
distribution of $Y\,|\,X_{1},\dots,X_{p}$.
- •
We apply the new procedure to real data from a case-control
study of Crohn’s disease in the United Kingdom, please see Section
[7](#S7). There, we show that the new knockoff method
makes twice as many discoveries as the original analysis of the
same data.

Before turning to the presentation of our method and results, we pause
to discuss the merits and limitations of our model, the relationships
between this work and others on selective inference, and the larger
problem of high-dimensional statistical testing.

### 1.3 Relationship with the classical setup for inference

It may seem to the statistician that our model appears rather
different from what she is used to. Our framework is, however,
not as exotic as it looks.

##### Classical setup

The usual setup for inference in conditional models is to assume a
strong parametric model for the response conditional on the
covariates, such as a homoscedastic linear model, but to assume as
little as possible about, or even condition on, the covariates. We do
the exact opposite by assuming we know *everything* about the
covariate distribution but *nothing* about the conditional
distribution $Y|X_{1},\dots,X_{p}$. Hence, we merely shift the
*burden of knowledge.* Our philosophy is, therefore,
to model $X$, not $Y$, whereas
classically, $Y$ (given $X$) is modeled and $X$ is not. In
practice, the parametric model in the classical approach is just an
approximation, and does not need to hold exactly to produce useful
inference. Analogously, we do not need to know the covariate
distribution exactly for our method to be useful, as we will
demonstrate in Sections [6](#S6) and [7](#S7).

##### When are our assumptions useful?

We do not claim our assumptions will always be appropriate, but there
are important cases when it is reasonable to think we know much more
about the covariate distribution than about the conditional
distribution of the response, including:

- •
When we in fact know exactly the covariate distribution because
we control it, such as in gene knockout experiments
(Cong et al.,, ; Peters et al.,, ), genetic crossing experiments (Haldane and Waddington,, ),
or sensitivity analysis of numerical models (Saltelli et al.,, ) (for
example climate models). In some cases we may also essentially know
the covariate distribution even when we do not control it, such as
in admixture mapping (Tang et al.,, ).
- •
When we have a large amount of unsupervised data (covariate data without
corresponding responses/labels) in addition to the $n$ labeled
observations. This is not uncommon in genetic or economic studies,
where many other studies will exist that have collected the same
covariate information but different response variables.
- •
When we simply have considerably more prior information about
the covariates than about the response. Indeed, the point of many
conditional modeling problems is to relate a poorly-understood
response variable to a set of well-understood covariates. For
instance, in genetic case-control studies, scientists seek to
understand the genetic basis of an extremely biologically-complex disease
using many comparatively simple single nucleotide polymorphisms
(SNPs) as covariates.

##### Payoff

There are substantial payoffs to our
framework. Perhaps the main advantage is the ability to use the
knockoff framework in high dimensions, a setting that was impossible
using the original approach. Even in high dimensions, previous inference results
rely not only on a parametric model that is often linear and
homoscedastic, but also on the sparsity or ultra-sparsity of the
parameters of that model in order to achieve some asymptotic
guarantee. In contrast, our framework can accommodate *any* model
for both the response and the covariates, and our guarantees are exact
in finite samples (non-asymptotic). In particular, our setup
encompasses any regression, classification, or survival model, including any
generalized linear model (GLM), and allows for arbitrary
nonlinearities and heteroscedasticity, such as are found in many
machine learning applications.

### 1.4 Relationship with work on inference after selection

There is a line of work on inference after selection, or
post-selection inference, for high-dimensional regression, the goal
of which is to first perform selection to make the problem
low-dimensional, and then produce p-values that are valid
*conditional* on the selection step
(Berk et al.,, ; Lockhart et al.,, ; Lee et al.,, ). These works differ from
ours in a number of ways so that we largely see them as
complementary activities.

- •
First, our focus is on selecting the right variables, whereas
the goal of this line of work is to adjust inference after some
selection has taken place. In more detail, these works presuppose a
selection procedure has been chosen (for reasons that may have
nothing to do with controlling Type I error) and then compute
p-values for the selected variables, taking into account the
selection step. In contrast, MX knockoffs is by itself a selection
procedure that controls Type I error.
- •
Second, inference after selection relies heavily on parametric
assumptions about the conditional distribution, namely that
$Y|X_{1},\dots,X_{p}\sim\mathcal{N}(\mu(X_{1},\dots,X_{p}),\sigma^{2}),$
making it unclear how to extend it to the more general setting of
the present paper.
- •
The third difference stems from their objects of inference. In
the selection step, a subset of size $m\leq n$ of the original $p$
covariates is selected, say $X_{j_{1}},\dots,X_{j_{m}}$, and the objects
of inference are the coefficients of the $X_{j_{k}}$’s in the
projection of $\mu$ onto the linear subspace spanned by the
$n\times m$ matrix of observed values of these $X_{j_{k}}$’s. That is,
the $k$th null hypothesis is that the aforementioned coefficient on
$X_{j_{k}}$ is zero—note that whether or not inference on the $j$th
variable is produced at all, and if it is, the object of that
inference, both depend on the initial selection step. In contrast,
if MX knockoffs were restricted to the homoscedastic Gaussian
model above, the $j$th null hypothesis would be that $\mu$ does not
depend on $X_{j}$, and there would be $p$ such null hypotheses, one
for each of the original variables.

### 1.5 Obstacles to obtaining p-values

Our procedure does not follow the canonical approach to FDR control
and multiple testing in general. The canonical approach is to plug
p-values into the BHq procedure, which controls the FDR under p-value
independence and certain forms of dependence
(Benjamini and Hochberg,, ; Benjamini and Yekutieli,, ).
Although these works have seeded a wealth of methodological
innovations over the past two decades (Benjamini,, ), all
these procedures act on a set of valid p-values (or equivalent
statistics), which they assume can be computed.(^1^11 Benjamini and Gavrilov, () and Bogdan et al., () transform the p-value cutoffs of common FDR-controlling procedures into penalized regression analogues to avoid p-values altogether. They only provably control the FDR in homoscedastic linear regression when the design matrix has orthogonal columns (necessitating, importantly, that $n\geq p$) but the latter work empirically retains control more generally whenever the signal obeys sparsity constraints. In a very different setting with spatial hypotheses, Li et al., () use approximations from Gaussian random field theory to heuristically control the FDR.)
The requirement of having valid p-values is quite constraining
for general conditional modeling problems.

#### 1.5.1 Regression p-value approximations

In low-dimensional ($n\geq p$) homoscedastic Gaussian linear
regression, p-values can be computed exactly even if the error
variance is unknown, although the p-values will not in general have
any simple dependence properties like independence or positive
regression dependency on a subset (PRDS). Already
for just the slightly broader class of low-dimensional GLMs, one must
resort to asymptotic p-values derived from maximum-likelihood theory,
which we will show in Section [5.1](#S5.SS1) can be far from valid
in practice.
In high-dimensional ($n<p$) GLMs, it is not clear how to get
p-values at all. Although some work (see for example
van de Geer et al., ()) exists on computing asymptotic p-values
under strong sparsity assumptions (usually the number of
important variables must be $o(\sqrt{n}/\log(p))$), like
their low-dimensional maximum-likelihood counterparts, these
methods suffer from highly non-uniform null p-values in many
finite-sample problems (see, for example, simulations in
Dezeure et al., ()). For binary covariates, the causal
inference literature uses matching and propensity scores for
approximately valid inference, but extending these methods to
high dimensions is still a topic of current research,
requiring similar assumptions and asymptotic approximations to
the aforementioned high-dimensional GLM literature
(Athey et al.,, ). Moving beyond generalized linear models
to the nonparametric setting, there exist measures of feature
importance, but no p-values.(^2^22In their online description of random forests ([http://www.math.usu.edu/~adele/forests/](http://www.math.usu.edu/~adele/forests/)), Leo Breiman and Adele Cutler propose a way to obtain a “z-score” for each variable, but without any theoretical distributional justification, and Strobl and Zeileis, () find “that the suggested test is not appropriate for statements of significance.”)

#### 1.5.2 Marginal testing

Faced with the inability to compute p-values for hypothesis tests of
conditional independence, one solution is to use *marginal*
p-values, i.e., p-values for testing *un*conditional (or
marginal) independence between $Y$ and $X_{j}$. This simplifies the
problem considerably, and many options exist for obtaining valid
p-values for such a test. However, marginal p-values are in general
invalid for testing conditional independence, and replacing
tests of conditional independence with tests of unconditional
independence is often undesirable. Indeed when $p\ll n$, so that
classical (e.g., maximum-likelihood) inference techniques for
regression give valid p-values for parametric tests of conditional
independence, it would be very unusual to resort to marginal testing
to select important covariates, and we cannot think of a textbook
that takes this route. Furthermore, the class of conditional test
statistics is far richer than that of marginal ones, and includes the
most powerful statistical inference and prediction methodology
available. For example, in compressed sensing, the signal recovery
guarantees for state-of-the-art $\ell_{1}$-based (joint) algorithms are
stronger than any guarantees possible with marginal methods. To
constrain oneself to marginal testing is to completely ignore the vast
modern literature on sparse regression that, while lacking
finite-sample Type I error control, has had tremendous success
establishing other useful inferential guarantees such as model
selection consistency under high-dimensional asymptotics in both
parametric (e.g., Lasso (Zhao and Yu,, ; Candès and Plan,, )) and
nonparametric (e.g., random forests (Wager and Athey,, ))
settings. Realizing this, the statistical genetics community has
worked on a number of multivariate approaches to improve power in
genome-wide association studies using both penalized
(Wu et al.,, ; He and Lin,, ) and Bayesian regression
(Guan and Stephens,, ; Li et al.,, ), but both approaches still suffer from a
lack of Type I error control (without making strong assumptions on
parameter priors). We will see that the MX knockoff procedure
is able to leverage the power of any of these techniques while adding
rigorous finite-sample Type I error control when the covariate distribution is known.

Some specific drawbacks of using marginal p-values are:

Even when the covariates are independent, so that
the hypotheses of conditional and unconditional independence
coincide, p-values resulting from marginal testing procedures may be
less powerful than those from conditional testing procedures. This
phenomenon has been reported previously, for example in statistical
genetics by Hoggart et al., () and many others. Intuitively, this is
because a joint model in $X_{1},\dots,X_{p}$ for $Y$ will have less
residual variance than a marginal model in just $X_{j}$.
There are exceptions, for instance if there is only one important
variable, then its marginal model is the correct joint model and a
conditional test will be less powerful due to the uncertainty in how
the other variables are included in the joint model. But in general,
marginal testing becomes increasingly underpowered relative to
conditional testing as the *absolute* number of important
covariates increases (Frommlet et al.,, ), suggesting particular
advantage for conditional testing in modern applications with
complex high-dimensional models.

There are also cases in which important variables are in fact fully
marginally independent of the response. As a toy example, if
$X_{1},X_{2}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\text{Bernoulli}(0.5)$ and $Y={\mathbbm{1}}_{{\{X_{1}+X_{2}=1\}}}$,
then $Y$ is marginally independent of each of $X_{1}$ and $X_{2}$, even
though together they determine $Y$ perfectly. $X_{1}$ and $X_{2}$ are
certainly *conditionally* dependent on $Y$, however, so a
conditional test can have power to discover them.

When the covariates are not independent,
marginal and conditional independence do not correspond, and we end
up asking the wrong question. For example, in a model with a few
important covariates which determine the outcome, and many
unimportant covariates which have no influence on the outcome but
are correlated with the important covariates, marginal testing
will treat such unimportant covariates as important. Thus, because
marginal testing is testing the wrong hypotheses, there will be many
“discoveries” which have no influence on the outcome.

The argument is often made, especially in genetics, that although
discovering an unimportant covariate just because it was
correlated with an important one is technically incorrect, it
can still be useful as it suggests that there is an important
covariate correlated with the discovered one. While this is indeed
useful, especially in genetics where correlated SNPs tend to be very
close to one another on the genome, this comes at a price since it
significantly alters the meaning of the FDR. Indeed, if we adopt
this viewpoint, the unit of inference is no longer a SNP but,
rather, a region on the genome, yet FDR is still being tested at the
SNP level. A consequence of this mismatch is that the result of the
analysis may be completely misleading, as beautifully argued in
Pacifico et al., (); Benjamini and Heller, (); Siegmund et al., (), see also
Chouldechova, () and Brzyski et al., () for later
references.

Marginal p-values will, in general, have
quite complicated joint dependence, so that BHq does not control FDR
exactly. Although procedures for controlling FDR under arbitrary
dependence exist, their increased generality tends to make them
considerably more conservative than BHq. In practice, however, the
FDR of BHq applied to dependent p-values is usually below its
nominal level, but the problem is that it can have highly
*variable* FDP. Recall that FDR is the expectation of FDP, the
latter being the random quantity we actually care about but cannot
control directly. Therefore, FDR control is only useful if the
realized FDP is relatively concentrated around its expectation, and
it is well-established (Efron,, , Chapter 4) that under
correlations BHq can produce highly skewed FDP distributions. In
such cases, with large probability, $\text{FDP}=0$ perhaps because
no discoveries are made, and when discoveries are made, the FDP may
be much higher than the nominal FDR, making it a misleading error
bound.

### 1.6 Getting valid p-values via conditional randomization
testing

If we insist on obtaining p-values for each $X_{j}$, there is in fact a
simple method when the covariate distribution is assumed known, as it
is in this paper. This method is similar in spirit to both propensity
scoring (where the distribution of a binary $X_{j}$ conditional on the
other variables is often estimated) and randomization/permutation
tests (where $X_{j}$ is either the only covariate or fully independent
of the other explanatory variables). Explicitly, a conditional
randomization test for the $j$th variable proceeds by first computing
some feature importance statistic $T_{j}$ for the $j$th variable. Then
the null distribution of $T_{j}$ can be computed through simulation by
independently sampling $X^{*}_{j}$’s from the *conditional*
distribution of $X_{j}$ given the others (derived from the known $F_{X}$)
and recomputing the same statistic $T^{*}_{j}$ with each new $X^{*}_{j}$ in
place of $X_{j}$, see Section [4](#S4) for details. Despite its
simplicity, we have not seen this test proposed previously in the
literature, although it nearly matches the usual randomization test
when the covariates are independent of one another.

### 1.7 Outline of the paper

The remainder of the paper is structured as follows:

- •
Section [2](#S2) frames the controlled
selection problem in rigorous mathematical terms.
- •
Section [3](#S3) introduces the model-X knockoff
procedure, examines its relationship with the earlier proposal
of Barber and Candès, (), proposes knockoff constructions and
feature statistics, and establishes FDR control.
- •
Section [4](#S4) introduces the conditional randomization test.
- •
Section [5](#S5) demonstrates through simulations
that the MX knockoff procedure controls the FDR in a number
of settings where no other procedure does, and that when
competitors exist, knockoffs is more powerful.
- •
Section [6](#S6) gives some preliminary simulations using
artificial and real data regarding the robustness of MX
knockoffs to unknown/misspecified covariate distributions.
- •
Section [7](#S7) applies our procedure to a
case-control study of Crohn’s disease in the United Kingdom.
- •
Section [8](#S8) concludes the paper with
extensions and potential lines of future research.

## 2 Problem statement

To state the controlled variable selection problem carefully, suppose
we have $n$ i.i.d. samples from a population, each of the form
$(X,Y)$, where $X=(X_{1},\ldots,X_{p})\in\mathbb{R}^{p}$ and $Y\in\mathbb{R}$. If the conditional distribution of $Y$ actually depends upon
a smaller subset of these variables, we would like to classify each
variable as relevant or not depending on whether it belongs to this
subset or not. Mathematically speaking, we are looking for the Markov
blanket $\mathcal{S}$, i.e. the “smallest” subset $\mathcal{S}$ such
that conditionally on $\{X_{j}\}_{j\in\mathcal{S}}$, $Y$ is
independent of all other variables.
For almost all joint
distributions of $(X,Y)$, there exists a unique Markov blanket but
there are pathological cases where it does not. An example is this:
suppose that $X_{1}$ and $X_{2}$ are independent Gaussian variables and
that $X_{3}=X_{1}-X_{2}$. Further assume that the distribution of $Y$
depends upon the vector $X$ only through $X_{1}+X_{2}$,
e.g., $Y\,|X\sim\mathcal{N}(X_{1}+X_{2},1)$. Then the set of
relevant variables—or equivalently, the Markov blanket—is ill
defined since we can say that the likelihood of $Y$ depends upon $X$
through either $(X_{1},X_{2})$, $(X_{1},X_{3})$, or $(X_{2},X_{3})$, all these
subsets being equally good. In order to define a unique set of
relevant variables, we shall work with the notion of conditional
*pairwise* independence.

###### Definition 2.1 .

A variable $X_{j}$ is said to be “null” if and only if $Y$ is
independent of $X_{j}$ conditionally on the other variables
$X_{{\text{-}j}}=\{X_{1},\ldots X_{p}\}\setminus\{X_{j}\}$. The subset of null
variables is denoted by $\mathcal{H}_{0}\subset\{1,\ldots p\}$ and
we call a variable $X_{j}$ “nonnull” or relevant if
$j\notin\mathcal{H}_{0}$.

From now on, our goal is to discover as many relevant
(conditionally dependent) variables as possible while keeping the
FDR under control.(^3^33Using the methods of Janson and Su, (), other error rates such as the $k$-familywise error rate can also be controlled using MX knockoffs, but we focus on FDR for this paper.) Formally, for a selection rule that
selects a subset $\hat{\mathcal{S}}$ of the covariates,

$$ $\text{FDR}:=\mathbb{E}\,\,\left[\frac{|\hat{\mathcal{S}}\cap\mathcal{H}_{0}|}{|\hat{\mathcal{S}}|}\right].$ (2.1) $$

In the example above, because of the perfect functional relationship
$X_{3}=X_{2}-X_{1}$, all three variables $X_{1},X_{2},X_{3}$ would be
classified as nulls. Imagine, however, breaking this relationship by
adding a bit of noise, e.g., $X_{3}=X_{2}-X_{1}+Z$, where $Z$ is
Gaussian noise (independent of $X_{1}$ and $X_{2}$) however small. Then
according to our definition, $X_{1}$ and $X_{2}$ are both nonnull while
$X_{3}$ is null—and everything makes sense. Having said this, we
should not let ourselves be distracted by such subtleties. In the
literature on graphical models there, in fact, exist weak regularity
conditions that guarantee that the (unique) set of relevant variables
defined by pairwise conditional independence, exactly coincides with
the Markov blanket so that there is no ambiguity. In this field,
researchers typically assume these weak regularity conditions hold
(examples would include the local and global Markov properties), and
proceed from there. For example, the textbook Edwards, ()
describes these properties on page 8 as holding “under quite general
conditions” and then assumes them for the rest of the book.

Our definition is very natural to anyone working with parametric
GLMs. In a GLM, the response $Y$ has a
probability distribution taken from an exponential family, which
depends upon the covariates only through the linear combination
$\eta=\beta_{1}X_{1}+\cdots+\beta_{p}X_{p}$. The relationship between
$Y$ and $X$ is specified via a link function $g$ such that
$\mbox{$\mathbb{E}\left(Y|X\right)$}=g^{-1}(\eta)$. In such models and under broad conditions,
$Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{j}\,|\,X_{\text{-}j}$ if and only if
$\beta_{j}=0$. In this context, testing the hypothesis that $X_{j}$ is a
null variable is the same as testing $H_{j}:\,\beta_{j}=0$.

###### Proposition 2.2 .

Take a family of random variables $X_{1},\ldots,X_{p}$ such that one
cannot perfectly predict any one of them from knowledge of the
others. If the likelihood of $Y$ follows a GLM, then
$Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{j}\,|\,X_{\text{-}j}$ if and only if
$\beta_{j}=0$. Hence, $\mathcal{H}_{0}$ from Definition [2.1](#S2.Thmtheorem1)
is exactly the set $\{j:\beta_{j}=0\}$.

###### Proof.

We prove this in the case of the logistic regression model as the
general case is similar. Here, the conditional distribution of $Y$ is
Bernoulli with

$$ $\mbox{$\mathbb{E}\left(Y|X\right)$}=\mathbb{P}(Y=1|X)=\frac{e^{\eta}}{1+e^{\eta}}=g^{-1}(\eta),\quad\eta=\beta_{1}X_{1}+\cdots+\beta_{p}X_{p},$ $$

and please note that the assumption about the covariates implies that
the model is identifiable. Now assume first that $\beta_{j}=0$. Then

$$ $\displaystyle p_{Y,X_{j}|X_{{\text{-}j}}}(y,x_{j}\,|\,x_{{\text{-}j}})$ $\displaystyle=p_{Y|X_{j},X_{{\text{-}j}}}(y\,|\,x_{j},x_{{\text{-}j}})\,p_{X_{j}|X_{{\text{-}j}}}(x_{j}\,|\,x_{{\text{-}j}})$ (2.2) $$

and since the first factor in the right-hand side does not depend on
$X_{j}$, we see that the conditional probability distribution function
factorizes. This implies conditional independence. In the other
direction, assume that $Y$ and $X_{j}$ are conditionally
independent. Then the likelihood function

$$ $\frac{\exp(Y{(\beta_{1}X_{1}+\cdots+\beta_{p}X_{p})})}{1+\exp({\beta_{1}X_{1}+\cdots+\beta_{p}X_{p}})}$ $$

must, conditionally on $X_{{\text{-}j}}$, factorize into a function of $Y$
times a function of $X_{j}$. A consequence of this is that conditionally
on $X_{{\text{-}j}}$, the odds ratio must not depend on $X_{j}$ (it must be
constant). However, this ratio is equal to $\exp(\beta_{j}X_{j})$ and is
constant only if $\beta_{j}=0$ since, by assumption, $X_{j}$ is not
determined by $X_{{\text{-}j}}$.
∎

The assumption regarding the covariates is needed. Indeed, suppose
$X_{1}\sim\mathcal{N}(0,1)$, $X_{2}=1\{X_{1}>0\}$, and $Y$ follows a
logistic model as above with $\eta=X_{1}+X_{2}$. Then
$Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{2}\,|\,X_{1}$ even though $\beta_{2}=1$. In this
example, the conditional distribution of $Y$ depends on $(X_{1},X_{2})$
only through $X_{1}$. Therefore, for the purpose of identifying
important variables (recall our task is to find important variables
and not to learn exactly how the likelihood function depends upon these
variables), we would like to find $X_{1}$ and actually do not care about
$X_{2}$ since it provides no new information.

###### Definition 2.1 .

###### Proposition 2.2 .

###### Proof.

## 3 Methodology

### 3.1 Model- X knockoffs

#### 3.1.1 Definition

###### Definition 3.1 .

Model-X knockoffs for the family of random
variables $X=(X_{1},\ldots,X_{p})$ are a new family of random
variables $\tilde{X}=(\tilde{X}_{1},\ldots,\tilde{X}_{p})$ constructed
with the following two properties: (1) for any subset
$S\subset\{1,\ldots,p\}$,(^4^44$\,{\buildrel d\over{=}}\,$ denotes equality in distribution, and the definition of the swapping operation is given just below.)

$$ $(X,\,\tilde{X})_{\textnormal{swap}{(S)}}\,\,{\buildrel d\over{=}}\,\,(X,\,\tilde{X});$ (3.1) $$

(2) $\tilde{X}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y\,|\,X$ if there is a response $Y$.
(2) is guaranteed if $\tilde{X}$ is constructed without looking at $Y$.

Above, the vector $(X,\,\tilde{X})_{\textnormal{swap}{(S)}}$ is obtained from
$(X,\,\tilde{X})$ by swapping the entries $X_{j}$ and $\tilde{X}_{j}$
for each $j\in S$; for example, with $p=3$ and $S=\{2,3\}$,

$$ $(X_{1},X_{2},X_{3},\tilde{X}_{1},\tilde{X}_{2},\tilde{X}_{3})_{\textnormal{swap}{(\{2,3\})}}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\,{\buildrel d\over{=}}\,}\,(X_{1},\tilde{X}_{2},\tilde{X}_{3},\tilde{X}_{1},X_{2},X_{3}).$ $$

We see from ([3.1](#S3.E1)) that original and knockoff variables
are pairwise exchangeable: taking any subset of variables and
swapping them with their knockoffs leaves the joint distribution
invariant. Note that our exchangeability condition is on
the covariates, and thus bears little resemblance to
exchangeability conditions for closed permutation testing
(see, e.g., Westfall and Troendle, ()). To give an example of MX knockoffs, suppose that
$X\sim\mathcal{N}(0,\bm{\Sigma})$. Then a joint distribution obeying
([3.1](#S3.E1)) is this:

$$ $(X,\,\tilde{X})\sim\mathcal{N}(0,\bm{G}),\quad\text{where}\quad\bm{G}=\begin{bmatrix}\bm{\Sigma}&\bm{\Sigma}-\operatorname{diag}\{s\}\\ \bm{\Sigma}-\operatorname{diag}\{s\}&\bm{\Sigma}\end{bmatrix};$ (3.2) $$

here, $\operatorname{diag}\{s\}$ is any diagonal matrix selected in such a
way that the joint covariance matrix $\bm{G}$ is positive
semidefinite. Indeed, the distribution obtained by swapping
variables with their knockoffs is Gaussian with a covariance
given by $\bm{P}\bm{G}\bm{P}$, where $\bm{P}$ is the
permutation matrix encoding the swap. Since
$\bm{P}\bm{G}\bm{P}=\bm{G}$ for any swapping operation, the
distribution is invariant. For an interesting connection with
the invariance condition in Barber and Candès, (), see
Appendix [A](#A1).

We will soon be interested in the problem of constructing knockoff
variables, having observed $X$. In the above example, a possibility is
to sample the knockoff vector $\tilde{X}$ from the conditional
distribution

$$ $\tilde{X}\,|X\,{\buildrel d\over{=}}\,\mathcal{N}(\mu,\bm{V}),$ $$

where $\mu$ and $\bm{V}$ are given by classical regression formulas,
namely,

$$ $\displaystyle\mu$ $\displaystyle=X-X{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\Sigma}^{-1}\operatorname{diag}\{s\}},$ $\displaystyle\bm{V}$ $\displaystyle=2\operatorname{diag}\{s\}-\operatorname{diag}\{s\}\bm{\Sigma}^{-1}\operatorname{diag}\{s\}.$ $$

There are, of course, many other ways of constructing knockoff
variables, and for the time being, we prefer to postpone the
discussion of more general constructions.

In the setting of the paper, we are given i.i.d. pairs
$(X_{i1},\ldots,X_{ip},Y_{i})\in\mathbb{R}^{p}\times\mathbb{R}$ of covariates and
responses, which we can assemble in a data matrix $\bm{X}$ and a data
vector $y$ in such a way that the $i$th row of $\bm{X}$ is
$(X_{i1},\ldots,X_{ip})$ and the $i$th entry of $y$ is $Y_{i}$. Then
the MX knockoff matrix $\tilde{\bm{X}}$ is constructed in such a
way that for each observation label $i$,
$(\tilde{X}_{i1},\ldots,\tilde{X}_{ip})$ is a knockoff for
$(X_{i1},\ldots,X_{ip})$ as explained above; that is to say, the
joint vector
$(X_{i1},\ldots,X_{ip},\tilde{X}_{i1},\ldots,\tilde{X}_{ip})$
obeys the pairwise exchangeability property ([3.1](#S3.E1)).

###### Definition 3.1 .

#### 3.1.2 Exchangeability of null covariates and their knockoffs

A crucial property of MX knockoffs is that we can swap null
covariates with their knockoffs without changing the joint
distribution of the original covariates $X$ and their knockoffs
$\tilde{X}$, conditional on the response $Y$. From now on, $X_{i:j}$ for
$i\leq j$ is a shorthand for $(X_{i},\ldots,X_{j})$.

###### Lemma 3.2 .

Take any subset $S\subset\mathcal{H}_{0}$ of nulls. Then

$$ ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}[\bm{X},\,\tilde{\bm{X}}]\mid y\,\,{\buildrel d\over{=}}\,\,[\bm{X},\,\tilde{\bm{X}}]_{\textnormal{swap}{(S)}}\mid y.}$ $$

###### Proof.

Since $y$’s marginal distribution is the same on both
sides of the equation, it is equivalent to show that
$([\bm{X},\,\tilde{\bm{X}}],\,y)\,\,{\buildrel d\over{=}}\,\,([\bm{X},\,\tilde{\bm{X}}]_{\textnormal{swap}{(S)}},\,y)$, which
is how we proceed. Assume without loss of generality
that $S=\{1,2,\ldots,m\}$. By row independence, it
suffices to show that $((X,\tilde{X}),\,Y)$ $\,{\buildrel d\over{=}}\,$
$((X,\tilde{X})_{\textnormal{swap}{(S)}},\,Y)$, where $X$ (resp. $Y$)
is a row of $\bm{X}$ (resp. $y$). Furthermore, since
$(X,\tilde{X})$ $\,{\buildrel d\over{=}}\,$ $(X,\tilde{X})_{\textnormal{swap}{(S)}}$, we only
need to establish that

$$ $Y\,|\,(X,\tilde{X})_{\textnormal{swap}{(S)}}\,\,{\buildrel d\over{=}}\,\,Y\,|\,(X,\tilde{X}).$ (3.3) $$

Letting $p_{Y|X}(y|x)$ be the conditional
distribution of $Y$, observe that

$$ $\displaystyle p_{Y|(X,\tilde{X})_{\textnormal{swap}{(S)}}}(y|(x,\tilde{x}))$ $\displaystyle=p_{Y|(X,\tilde{X})}(y|(x,\tilde{x})_{\textnormal{swap}{(S)}})$ $\displaystyle=p_{Y|X}(y|x^{\prime}),$ $$

where $x^{\prime}_{i}=\tilde{x}_{i}$ if $i\in S$ and $x^{\prime}_{i}=x_{i}$
otherwise. The second equality above comes from the fact that $Y$ is
conditionally independent of $\tilde{X}$ by property (2) in the
definition of MX knockoffs. Next, since $Y$ and $X_{1}$ are independent
conditional on $X_{2:p}$, we have

$$ $\displaystyle p_{Y|X_{1:p}}(y|\tilde{x}_{1},x^{\prime}_{2:p})$ $\displaystyle=p_{Y|X_{2:p}}(y|x^{\prime}_{2:p})$ $\displaystyle=p_{Y|X_{1:p}}(y|x_{1},x^{\prime}_{2:p}).$ $$

This shows that

$$ $Y\,|\,(X,\tilde{X})_{\textnormal{swap}{(S)}}\,\,{\buildrel d\over{=}}\,\,Y\,|\,(X,\tilde{X})_{\textnormal{swap}{(S\setminus\{1\})}}.$ $$

We can repeat this argument with the second variable, the third, and so on until $S$ is empty. This proves ([3.3](#S3.E3)).
∎

###### Lemma 3.2 .

###### Proof.

### 3.2 Feature statistics

In order to find the relevant variables, we now compute statistics
$W_{j}$ for each $j\in\{1,\ldots,p\}$, a large positive value of
$W_{j}$ providing evidence against the hypothesis that $X_{j}$ is null.
This statistic depends on the response and the original variables
but also on the knockoffs; that is,

$$ $W_{j}=w_{j}([\bm{X},\,\tilde{\bm{X}}],\,y)$ $$

for some function $w_{j}$. As in Barber and Candès, (), we impose a flip-sign property, which says that swapping the $j$th variable
with its knockoff has the effect of changing the sign of
$W_{j}$. Formally, if $[\bm{X},\,\tilde{\bm{X}}]_{\textnormal{swap}(S)}$ is the matrix
obtained by swapping columns in $S$,

$$ $w_{j}([\bm{X},\,\tilde{\bm{X}}]_{\textnormal{swap}(S)},\,y)=\begin{cases}w_{j}([\bm{X},\,\tilde{\bm{X}}],y),&j\not\in S,\\ -w_{j}([\bm{X},\,\tilde{\bm{X}}],y),&j\in S.\end{cases}$ (3.4) $$

In contrast to the aforementioned work, we do not require the sufficiency
property that $w_{j}$ depend on $\bm{X}$, $\tilde{\bm{X}}$, and $y$ only
through $[\bm{X},\,\tilde{\bm{X}}]^{\top}[\bm{X},\,\tilde{\bm{X}}]$ and
$[\bm{X},\,\tilde{\bm{X}}]^{\top}y$.

At this point, it may help the reader unfamiliar with the knockoff
framework to think about knockoff statistics $W=(W_{1},\ldots,W_{p})$
in two steps: first, consider a statistic $T$ for each original and
knockoff variable,

$$ $T\triangleq(Z,\,\tilde{Z})=(Z_{1},\ldots,Z_{p},\tilde{Z}_{1},\ldots,\tilde{Z}_{p})=t([\bm{X},\,\tilde{\bm{X}}],\,y),$ $$

with the idea that $Z_{j}$ (resp. $\tilde{Z}_{j}$) measures the importance
of $X_{j}$ (resp. $\tilde{X}_{j}$). Assume the natural property that
switching a variable with its knockoff simply switches the
components of $T$ in the same
way, namely, for each $S\subset\{1,\ldots,p\}$,

$$ $(Z,\,\tilde{Z})_{\textnormal{swap}{(S)}}=t([\bm{X},\,\tilde{\bm{X}}]_{\textnormal{swap}{(S)}},\,y).$ (3.5) $$

Then one can construct a $W_{j}$ obeying the flip-sign condition
([3.4](#S3.E4)) by setting

$$ $W_{j}=f_{j}(Z_{j},\tilde{Z}_{j}),$ $$

where $f_{j}$ is any anti-symmetric function.(^5^55An anti-symmetric function $f$ is such that $f(v,u)=-f(u,v)$.) (Conversely, any
statistic $W_{j}$ verifying the flip sign condition can be constructed
in this fashion.) Adopting this approach, consider a regression
problem and run the Lasso on the original design augmented with
knockoffs,

$$ $\text{min}_{b\in\mathbb{R}^{2p}}\quad\textstyle{\frac{1}{2}}\|y-[\bm{X},\,\tilde{\bm{X}}]b\|_{2}^{2}+\lambda\|b\|_{1}$ (3.6) $$

and denote the solution by $\hat{b}(\lambda)$ (the first $p$
components are the coefficients of the original variables and the last
$p$ are for the knockoffs). Then the Lasso
coefficient-difference (LCD) statistic sets
$Z_{j}=|\hat{b}_{j}(\lambda)|$, $\tilde{Z}_{j}=|\hat{b}_{j+p}(\lambda)|$, and

$$ $W_{j}=Z_{j}-\tilde{Z}_{j}=|\hat{b}_{j}(\lambda)|-|\hat{b}_{j+p}(\lambda)|.$ (3.7) $$

A large positive value of $W_{j}$ provides some evidence that the
distribution of $Y$ depends upon $X_{j}$, whereas under the null $W_{j}$
has a symmetric distribution and, therefore, is equally likely to take on
positive and negative values, as we shall see next. Before moving on,
however, please carefully observe that the value of $\lambda$ in
([3.7](#S3.E7)) does not need to be fixed in advance, and can be
computed from $y$ and $[\bm{X},\,\tilde{\bm{X}}]$ in any data-dependent fashion as long as
permuting the columns of $\bm{X}$ does not change its value; for
instance, it can be selected by cross-validation.

###### Lemma 3.3 .

Conditional on $(|W_{1}|,\dots,|W_{p}|)$, the signs of the null $W_{j}$’s,
$j\in\mathcal{H}_{0}$, are i.i.d. coin flips.

###### Proof.

Let $\epsilon=(\epsilon_{1},\ldots,\epsilon_{p})$ be a sequence of
independent random variables such that $\epsilon_{j}=\pm 1$ with
probability $1/2$ if $j\in\mathcal{H}_{0}$, and $\epsilon_{j}=1$
otherwise. To prove the claim, it suffices to establish that

$$ $W\,\,{\buildrel d\over{=}}\,\,\epsilon\odot W,$ (3.8) $$

where $\odot$ denotes pointwise multiplication, i.e.
$\epsilon\odot W=(\epsilon_{1}W_{1},\ldots,\epsilon_{p}W_{p})$. Now,
take $\epsilon$ as above and put
$S=\{j:\epsilon_{j}=-1\}\subset\mathcal{H}_{0}$. Consider
swapping variables in $S$:

$$ $W_{\textnormal{swap}(S)}\triangleq w([\bm{X},\,\tilde{\bm{X}}]_{\textnormal{swap}(S)},y).$ $$

On the one hand, it follows from the flip-sign property that
$W_{\textnormal{swap}(S)}=\epsilon\odot W$. On the other hand, Lemma
[3.2](#S3.Thmtheorem2) implies that $W_{\textnormal{swap}(S)}\,\,{\buildrel d\over{=}}\,\,W$ since
$S\subset\mathcal{H}_{0}$. These last two properties give ([3.8](#S3.E8)).
∎

In fact, since the pairwise exchangeability property of $[\bm{X},\,\tilde{\bm{X}}]$
holds conditionally on $y$ according to Lemma [3.2](#S3.Thmtheorem2),
the coin-flipping property also holds conditionally on $y$.

###### Lemma 3.3 .

###### Proof.

### 3.3 FDR control

From now on, our methodology follows that of Barber and Candès, () and we
simply rehearse the main ingredients while referring to their paper
for additional insights. It follows from Lemma [3.3](#S3.Thmtheorem3) that the
null statistics $W_{j}$ are symmetric and that for any fixed threshold
$t>0$,

$$ $\#\{j:W_{j}\leq-t\}\geq\#\{\text{null }j:W_{j}\leq-t\}\,{\buildrel d\over{=}}\,\#\{\text{null }j:W_{j}\geq t\}.$ $$

Imagine then selecting those variables such that $W_{j}$ is sufficiently
large, e.g., $W_{j}\geq t$, then the false discovery proportion (FDP)

$$ $\textnormal{FDP}(t)=\frac{\#\{\text{null }j:W_{j}\geq t\}}{\#\{j:W_{j}\geq t\}}$ (3.9) $$

can be estimated via the statistic

$$ $\widehat{\textnormal{FDP}}(t)=\frac{\#\{j:W_{j}\leq-t\}}{\#\{j:W_{j}\geq t\}}$ $$

since the numerator is an upward-biased estimate of the unknown
numerator in ([3.9](#S3.E9)). The idea of the knockoff procedure is to
choose a data-dependent threshold as liberal as possible while having
an estimate of the FDP under control. The theorem below shows that
estimates of the FDR process can be inverted to give tight FDR
control.

###### Theorem 3.4 .

Choose a threshold $\tau>0$ by setting(^6^66When we write $\min\{t>0:\ldots\}$, we abuse notation since we actually mean $\min\{t\in\mathcal{W}_{+}:\ldots\}$, where $\mathcal{W}_{+}=\{|W_{j}|:|W_{j}|>0\}$.)

$$ $\tau=\min\left\{t>0:\frac{\#\{j:W_{j}\leq-t\}}{\#\{j:W_{j}\geq t\}}\leq q\right\}\quad(\textbf{Knockoffs}),$ (3.10) $$

where $q$ is the target FDR level (or $\tau=+\infty$ if the set above
is empty). Then the procedure selecting the variables

$$ $\hat{S}=\{j:W_{j}\geq\tau\}\;$ $$

controls the
modified FDR defined as

$$ $\textnormal{mFDR}=\mathbb{E}\left[\frac{|\{j\in\hat{S}\cap\mathcal{H}_{0}\}|}{|\hat{S}|+1/q}\right]\leq q.$ $$

The slightly more conservative procedure, given by incrementing the
number of negatives by one,

$$ $\tau_{+}=\min\left\{t>0:\frac{1+\#\{j:W_{j}\leq-t\}}{\#\{j:W_{j}\geq t\}}\leq q\right\}\quad(\textbf{Knockoffs+})$ (3.11) $$

and setting $\hat{S}=\{j:W_{j}\geq\tau_{+}\}$, controls the usual FDR,

$$ $\mathbb{E}\left[\frac{|\{j\in\hat{S}\cap\mathcal{H}_{0}\}|}{|\hat{S}|\vee 1}\right]\leq q.$ $$

These results are non-asymptotic and hold no matter the dependence
between the response and the covariates—in
fact, they hold *conditionally* on the response $y$.

The proof is the same as that of Theorems 1 and 2 in
Barber and Candès, ()—and, therefore, omitted—since all we
need is that the null statistics have signs distributed as
i.i.d. coin flips (even conditionally on $y$). Note that
Theorem [3.4](#S3.Thmtheorem4) only tells one side of the story: Type
I error control; the other very important side is power, which
leads us to spend most of the remainder of the paper
considering how best to construct knockoff variables and
statistics.

###### Theorem 3.4 .

### 3.4 Constructing model- X knockoffs

#### 3.4.1 Exact constructions

We have seen in Section [3.1.1](#S3.SS1.SSS1) one way of constructing
MX knockoffs in the case where the covariates are Gaussian.
How should we proceed for non-Gaussian data? In this regard, the
characterization below may be useful.

###### Proposition 3.5 .

The random variables $(\tilde{X}_{1},\ldots,\tilde{X}_{p})$ are
model-X knockoffs for $(X_{1},\ldots,X_{p})$ if and only if for any
$j\in\{1,\ldots,p\}$, the pair $(X_{j},\tilde{X}_{j})$ is
exchangeable conditional on all the other variables and their
knockoffs (and, of course, $\tilde{X}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y\,|\,X$).

The proof consists of simple manipulations of the definition and is,
therefore, omitted. Our problem can thus also be posed as constructing pairs that are
conditionally exchangeable. If the components of the vector $X$ are
independent, then any independent copy of $X$ would work; that is, any
vector $\tilde{X}$ independently sampled from the same joint
distribution as $X$ would work. With dependent coordinates, we may
proceed as follows:

Figure: Algorithm 1 Sequential Conditional Independent Pairs.

Above, $\mathcal{L}(X_{j}\,|\,X_{{\text{-}j}},\,\tilde{X}_{1:j-1})$ is the
conditional distribution of $X_{j}$ given $(X_{{\text{-}j}},\tilde{X}_{1:j-1})$.
When ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}p}=3$, this would work as follows: sample $\tilde{X}_{1}$ from
$\mathcal{L}(X_{1}\,|\,X_{2:3})$. Once this is done,
$\mathcal{L}(X_{1:3},\tilde{X}_{1})$ is available and we, therefore,
know $\mathcal{L}(X_{2}\,|\,X_{1},X_{3},\tilde{X}_{1})$. Hence, we can
sample $\tilde{X}_{2}$ from this distribution. Continuing,
$\mathcal{L}(X_{1:3},\tilde{X}_{1:2})$ becomes known and we can
sample $\tilde{X}_{3}$ from
$\mathcal{L}(X_{3}\,|\,X_{1:2},\tilde{X}_{1:2})$.

It is not immediately clear why Algorithm [1](#algorithm1) yields
a sequence of random variables obeying the exchangeability property
([3.1](#S3.E1)), and we prove this fact in
Appendix [B](#A2).
There is, of course, nothing special about the ordering in which
knockoffs are created and equally valid constructions may be
obtained by looping through an arbitrary ordering of the
variables. For example, in a data analysis application where we
would need to build a knockoff copy for each row of the design,
independent (random) orderings may be used.

To have power or, equivalently, to have a low Type II error rate, it
is intuitive that we would like to have original features $X_{j}$ and
their knockoff companions $\tilde{X}_{j}$ to be as “independent” as
possible.

We do not mean to imply that running Algorithm [1](#algorithm1) is
a simple matter. In fact, it may prove rather complicated since we
would have to recompute the conditional distribution at each step;
this problem is left for future research. Instead, in this paper we
shall work with approximate MX knockoffs and will demonstrate
empirically that for models of interest, such constructions yield FDR
control.

###### Proposition 3.5 .

#### 3.4.2 Approximate constructions: second-order model- X knockoffs

Rather than asking that $(X,\,\tilde{X})_{\textnormal{swap}{(S)}}$ and
$(X,\,\tilde{X})$ have the same distribution for any subset $S$, we
can ask that they have the same first two moments, i.e., the same mean
and covariance. Equality of means is a simple matter. As far as the
covariances are concerned, equality is equivalent to

$$ $\operatorname{cov}(X,\,\tilde{X})=\bm{G},\quad\text{where}\quad\bm{G}=\begin{bmatrix}\bm{\Sigma}&\bm{\Sigma}-\operatorname{diag}\{s\}\\ \bm{\Sigma}-\operatorname{diag}\{s\}&\bm{\Sigma}\end{bmatrix}.$ (3.12) $$

We, of course, recognize the same form as in ([3.2](#S3.E2))
where the parameter $s$ is chosen to yield a positive semidefinite
covariance matrix. (When $(X,\tilde{X})$ is Gaussian, a matching of
the first two moments implies a matching of the joint distributions so
that we have an exact construction.) Furthermore,
Appendix [A](#A1) shows that the same problem was already solved
in Barber and Candès, (), as the same constraint on $s$ applies but with
the empirical covariance replacing the true covariance. This means
that the same two constructions proposed in Barber and Candès, () are
just as applicable to second-order model-X knockoffs.

For the remainder of this section, we will assume the covariates have
each been translated and rescaled to have mean zero and variance
one. To review, the *equicorrelated* construction uses

$$ $s^{\text{EQ}}_{j}=2\lambda_{\text{min}}(\bm{\Sigma})\wedge 1\text{ for all }j,$ $$

which minimizes the correlation between variable-knockoff pairs
subject to the constraint that all such pairs must have the same
correlation. The *semidefinite program (SDP)* construction solves
the convex program

$$ $\begin{array}[]{rl}\text{minimize}&\quad\sum_{j}|1-s^{\text{SDP}}_{j}|\\ \text{subject to}&\quad s^{\text{SDP}}_{j}\geq 0\\ &\quad\text{diag}\{s^{\text{SDP}}\}\preceq 2\bm{\Sigma},\end{array}$ (3.13) $$

which minimizes the sum of absolute values of variable-knockoff
correlations among all suitable $s$.

In applying these constructions
to problems with large $p$, we run into some new difficulties:

- •
Excepting very specially-structured matrices like the identity,
$\lambda_{\text{min}}(\bm{\Sigma})$ tends to be extremely small as $p$
gets large. The result is that constructing equicorrelated
knockoffs in high dimensions, while fairly computationally easy,
will result in very low power, since all the original variables will
be nearly indistinguishable from their knockoff counterparts.
- •
For large $p$, ([3.13](#S3.E13)), while convex, is prohibitively
computationally expensive. However, if it could be computed, it
would produce much larger $s_{j}$’s than the equicorrelated
construction and thus be considerably more powerful.

To address these difficulties, we first generalize the two knockoff
constructions by the following two-step procedure, which we call
the approximate semidefinite program (ASDP) construction:

Choose an approximation $\bm{\Sigma}_{\text{approx}}$ of
$\bm{\Sigma}$ and solve:

$$ $\begin{array}[]{rl}\text{minimize}&\quad\sum_{j}|1-\hat{s}_{j}|\\ \text{subject to}&\quad\hat{s}_{j}\geq 0\\ &\quad\text{diag}\{\hat{s}\}\preceq 2\bm{\Sigma}_{\text{approx}}.\end{array}$ (3.14) $$

Solve:

$$ $\begin{array}[]{rl}\text{maximize}&\quad\gamma\\ \text{subject to}&\quad\text{diag}\{\gamma\hat{s}\}\preceq 2\bm{\Sigma},\end{array}$ (3.15) $$

and set $s^{\text{ASDP}}=\gamma\hat{s}$. Note this problem can be solved quickly
by, e.g., bisection search over $\gamma\in[0,1]$.

ASDP with $\bm{\Sigma}_{\text{approx}}=\bm{I}$ trivially gives $\hat{s}_{j}=1$
and $\gamma=2\lambda_{\text{min}}(\bm{\Sigma})\wedge 1$, reproducing the
equicorrelated construction. ASDP with $\bm{\Sigma}_{\text{approx}}=\bm{\Sigma}$
clearly gives $\hat{s}_{j}=s^{\text{SDP}}$ and $\gamma=1$, reproducing
the SDP construction. Note that the ASDP step 2 is always fast, so the
speed of the equicorrelated construction comes largely because the
problem *separates* into $p$ computationally independent SDP
subproblems of
$\text{min}|1-\hat{s}_{j}|\text{ s.t. }0\leq\hat{s}_{j}\leq 2$. However,
power is lost due to the very naïve approximation
$\bm{\Sigma}_{\text{approx}}=\bm{I}$ which results in a very small $\gamma$.

In general, we can choose $\bm{\Sigma}_{\text{approx}}$ to be an
$m$-block-diagonal approximation of $\bm{\Sigma}$, so that the ASDP from
Step 1 separates into $m$ smaller, more computationally tractable, and
trivially parallelizable SDP subproblems. If the approximation is
fairly accurate, we may also find that $\gamma$ remains large, so that
the knockoffs are nearly as powerful as if we had used the SDP
construction. We demonstrate the ASDP construction in
Section [7](#S7) when we analyze the Crohn’s disease data.

## 4 The conditional randomization test

This section presents an alternative approach to the controlled
variable selection problem. To describe our approach, it may be best
to consider an example. Assume we are in a regression setting and let
$\hat{b}_{j}(\lambda)$ be the value of the Lasso estimate of the $j$th
regression coefficient. We would like to use the statistic
$\hat{b}_{j}(\lambda)$ to test whether $Y$ is conditionally independent
of $X_{j}$ since large values of $|\hat{b}_{j}(\lambda)|$ provide evidence
against the null. To construct a test, however, we would need to know
the sampling distribution of $\hat{b}_{j}(\lambda)$ under the null
hypothesis that $Y$ and $X_{j}$ are conditionally independent, and it is
quite unclear how one would obtain such knowledge.

### 4.1 The test

A way out is to sample the covariate $X_{j}$ conditional on all the
other covariates (but not the response), where by “sample” we
explicitly mean to draw a new sample from the conditional distribution
of $X_{j}\,|\,X_{\text{-}j}$ using a random number generator. We then
compute the Lasso statistic $\hat{b}^{*}_{j}(\lambda)$, where the
^∗ superscript indicates that the statistic is computed from the artificially
sampled value of the covariate $X_{j}$. Now, under the null hypothesis
of conditional independence between $Y$ and $X_{j}$, it happens that
$\hat{b}^{*}_{j}(\lambda)$ and $\hat{b}_{j}(\lambda)$ are identically
distributed and that, furthermore, this statement holds true
conditional on $Y$ and all the other covariates. This claim is proved
in Lemma [4.1](#S4.Thmtheorem1) below. A consequence of this is that by
simulating a covariate conditional on the others, we can sample at
will from the conditional distribution of any test statistic and compute
p-values as described in Algorithm [2](#algorithm2).

Figure: Algorithm 2 Conditional Randomization Test.

###### Lemma 4.1 .

Let $(Z_{1},Z_{2},Y)$ be a triple of random variables, and
construct another triple $(Z_{1}^{*},Z_{2},Y)$ as

$$ $Z_{1}^{*}\,|\,(Z_{2},Y)\quad\,{\buildrel d\over{=}}\,\quad Z_{1}\,|\,Z_{2}.$ $$

Then under the null hypothesis $Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z_{1}\,|\,Z_{2}$, any test
statistic $T=t(Z_{1},Z_{2},Y)$ obeys

$$ $T\,|\,(Z_{2},Y)\quad\,{\buildrel d\over{=}}\,\quad T^{*}\,|\,(Z_{2},Y),$ $$

where $T^{*}=t(Z_{1}^{*},Z_{2},Y)$.

###### Proof.

To prove the claim, it suffices to show that $Z_{1}$ and $Z_{1}^{*}$ have
the same distribution conditionally on $(Z_{2},Y)$. This follows from

$$ $\displaystyle Z_{1}^{*}\,|\,(Y,Z_{2})\quad\,{\buildrel d\over{=}}\,\quad Z_{1}\,|\,Z_{2}\quad\,{\buildrel d\over{=}}\,\quad Z_{1}\,|\,(Z_{2},Y).$ $$

The first equality comes from the definition of $Z_{1}^{*}$ while the
second follows from the conditional independence of $Y$ and $Z_{1}$,
which holds under the null.
∎

The consequence of Lemma [4.1](#S4.Thmtheorem1) is that we can compute the 95%
percentile, say, of the conditional distribution of $T^{*}$ denoted
by $t^{*}_{0.95}(Z_{2},Y)$. Then by definition, under the null,

$$ $\mathbb{P}(T>t^{*}_{0.95}(Z_{2},Y)\,|\,(Z_{2},Y))\leq 0.05.$ $$

Since this equality holds conditionally, it also holds marginally.

###### Lemma 4.1 .

###### Proof.

### 4.2 Literature review

The conditional randomization test is most closely related to the propensity score
(Rosenbaum and Rubin,, ), which also uses the conditional distribution
$X_{j}\,|\,X_{\text{-}j}$ to perform inference on the conditional
relationship between $Y$ and $X_{j}$ given $X_{\text{-}j}$. However, propensity scores require
$X_{j}$ be binary, and the propensity score itself is normally
estimated, although Rosenbaum, () shows that when all the covariates
jointly take a small number of discrete values, propensity score
analysis can be done exactly. Doran et al., () also rely on
the data containing repeated observations of $X_{\text{-}j}$ so that
certain observations can be permuted nonparametrically while
maintaining the null distribution. In fact, the exact term
“conditional randomization test” has also been used in randomized
controlled experiments to test for independence of $Y$ and $X_{j}$
conditioned more generally on some function of $X_{\text{-}j}$ (such
as a measure of imbalance in $X_{\text{-}j}$ if $X_{j}$ is binary),
again relying on discreteness of the function so that there exist
permutations of $X_{j}$ which leave the function value
unchanged. Despite the similar name, our conditional
randomization test is quite distinct from these, as it does not rely
on discreteness or experimental control in any of the covariates.

Another line of work exists within the linear model regime, whereby
the null (without $X_{j}$) model is estimated and then the empirical
residuals are permuted to produce a null distribution for $Y$
(Freedman and Lane,, ). Because this approach is only exact when the
empirical residuals match the true residuals, it explicitly relies on
a parametric model for $Y\,|\,X_{\text{-}j}$, as well as the ability to
estimate it quite accurately.

### 4.3 Comparisons with knockoffs

One major limitation of the conditional randomization method is its
computational cost. It requires computing randomization p-values for
many covariates and to a high-enough resolution for
multiple-comparisons correction. Clearly, this requires samples in
the extreme tail of the p-value distribution. This means computing a
very large number of feature importance statistics $T_{j}$, each of
which can be expensive since for reasons outlined in the drawbacks
associated with marginal testing, powerful $T_{j}$’s will take into
account the full dimensionality of the model, e.g., absolute value of
the Lasso-estimated coefficient. In fact, the number of computations
of $T_{j}$, tallied over all $j$, required by the conditional
randomization method is $\Omega(p)$.(^7^77$a(N)\in\Omega(b(N))$ means that there exist $N_{0}$ and $C>0$ such that $a(N)\geq Cb(N)$ for $N\geq N_{0}$.) To see this, suppose for simplicity that all $R$
rejected p-values take on the value of half the BHq cutoff equal to
$\tau=qR/p$, and all we need to do is upper-bound them below
$\tau$. This means there are $R$ p-values $P_{j}$ for which plugging
$K=\infty$ into Algorithm [2](#algorithm2) would yield $P_{j}=\tau/2$.
After $K<\infty$ samples, the approximate p-value (ignoring the $+1$
correction) is distributed as $K^{-1}\,\operatorname{Bin}(K,P_{j})$.
We could then use this binomial count to construct a confidence
interval for $P_{j}$. A simple calculation shows that to be reasonably
confident that $P_{j}\leq\tau$, $K$ must be on the order of
at least $1/\tau$. Since there are $R$ such p-values, this justifies
the claim.

Note that for knockoffs, the analogous computation of $T$ need only
be done exactly once. If, for instance, each $T_{j}$ requires a Lasso
computation, then the conditional randomization test’s computational
burden is very challenging for medium-scale $p$ in the thousands
and prohibitive for large-scale (e.g., genetics) $p$ in the hundreds of
thousands or millions. We will see in Section [5.3.1](#S5.SS3.SSS1) that
there are power gains, along with huge computational costs, to be had
by using conditional randomization in place of knockoffs, and
Section [7](#S7) will show that the MX knockoff
procedure easily scales to large data sets.

Another advantage of MX knockoffs is its guaranteed control of
the FDR, whereas the BHq procedure does not offer strict control when
applied to arbitrarily dependent p-values.

## 5 Numerical simulations

In this section we demonstrate the importance, utility, and
practicality of MX knockoffs for high-dimensional
nonparametric conditional modeling.

### 5.1 Logistic regression p-values

Asymptotic maximum likelihood theory promises valid p-values for each
coefficient in a GLM only when $n\gg p$. However, these approximate
p-values can usually be computed as long as $n>p$, so a natural
question arising from high-dimensional applications is whether such
asymptotic p-values are valid when $n$ and $p$ are both large with
$p/n\geq 0.1$, for example. We simulated $10^{4}$ independent design
matrices ($n=500$, $p=200$) and binary responses from a logistic
regression for the following two settings:

- (1)
$(X_{1},\dots,X_{p})$ is an AR(1) time series with AR
coefficient $0.5$ and
$Y\,|\,X_{1},\dots,X_{p}\sim\text{Bernoulli}(0.5)$
- (2)
$(X_{1},\dots,X_{p})$ is an AR(1) time series with AR
coefficient $0.5$ and
$Y\,|\,X_{1},\dots,X_{p}\sim\text{Bernoulli}\left(\text{logit}\left(0.08(X_{2}+\cdots+X_{21})\right)\right)$

Figure: Figure 1: Distribution of null logistic regression p-values with $n=500$ and $p=200$; 10,000 replications.
Refer to caption: /html/1610.02351/assets/x1.png

Histograms for the p-values for $\beta_{1}$ (null in all cases) are
shown in Figure [1](#S5.F1). Both histograms are far from uniform,
and Table [1](#S5.T1) shows each distribution’s concentration
near zero. We see that the small quantiles have extremely inflated
probabilities—over 20 times nominal for
$\mathbb{P}\left\{{p\text{-value}\leq 0.1\%}\right\}$ in setting (2). We also see that the
exact null distribution depends on the unknown coefficient sequence
$\beta_{2},\dots,\beta_{p}$, since the probabilities between settings
differ statistically significantly at all three cutoffs.

To confirm that this non-uniformity is not just a finite-sample
effect, we also simulated $10^{4}$ i.i.d. $\mathcal{N}(0,1)$ design
matrices with independent Bernoulli$(0.5)$ responses for $n=500$,
$p=200$ and $n=5000$, $p=2000$ as settings (3) and (4),
respectively. Table [1](#S5.T1) shows that the distribution does
not really change as $n$ and $p$ are increased with constant
proportion.

**Table 1: Inflated p-value probabilities with estimated Monte Carlo standard errors in parentheses. See text for meanings of settings (1), (2), (3), (4).**
|  | (1) | (2) | (3) | (4) |
| --- | --- | --- | --- | --- |
| $\mathbb{P}\left\{{p\text{-value}\leq 5\%}\right\}$ | $16.89\%$ $(0.37\%)$ | $19.17\%$ $(0.39\%)$ | $16.88\%$ $(0.37\%)$ | $16.78\%$ $(0.37\%)$ |
| $\mathbb{P}\left\{{p\text{-value}\leq 1\%}\right\}$ | $6.78\%$ $(0.25\%)$ | $8.49\%$ $(0.28\%)$ | $7.02\%$ $(0.26\%)$ | $7.03\%$ $(0.26\%)$ |
| $\mathbb{P}\left\{{p\text{-value}\leq 0.1\%}\right\}$ | $1.53\%$ $(0.12\%)$ | $2.27\%$ $(0.15\%)$ | $1.87\%$ $(0.14\%)$ | $2.04\%$ $(0.14\%)$ |

These results show that the usual logistic regression p-values one
might use when $n\geq p$ can have null distributions that are quite far
from uniform, and even if one wanted to correct that distribution, it
depends in general on unknown problem parameters, further complicating
matters. When $n<p$ the problem becomes even more challenging, with
existing methods similarly asymptotic as well as requiring stringent sparsity assumptions
(van de Geer et al.,, ). Thus, despite the wealth of research on
controlling FDR, without a way to obtain valid p-values, even the problem of controlling
FDR in medium-to-high-dimensional GLMs remains unsolved.

### 5.2 Alternative knockoff statistics

As mentioned in Section [3.2](#S3.SS2), the new model-X
knockoffs framework allows for a wider variety of $W$
statistics to be used than in the fixed-X
framework. Choices of $Z_{j}$ include well-studied statistical
measures such as the coefficient estimated in a GLM, but can
also include much more ad-hoc/heuristic measures such as
random forest bagging feature importances or sensitivity
analysis measures such as the Monte-Carlo-estimated total
sensitivity index. By providing variable selection with rigorous
Type I error control for general models and statistics,
knockoffs can be used to improve the interpretability of
complex black-box supervised/machine learning models. There
are also many available choices for the anti-symmetric
function $f_{j}$, such as ${|Z_{j}|-|\tilde{Z}_{j}|}$,
${\operatorname*{sign}(|Z_{j}|-|\tilde{Z}_{j}|)\max\{|Z_{j}|,|\tilde{Z}_{j}|\}}$, or
${\log(|Z_{j}|)-\log(|\tilde{Z}_{j}|)}$.

The main point of this subsection is that knockoffs can be used as a wrapper around essentially *any* data-fitting or prediction algorithm, and regardless of the chosen algorithm still provides rigorous error control for variable selection. We discuss here a few appealing new options for statistics $W$, but defer full exploration of these very extensive possibilities to future work.

#### 5.2.1 Adaptive knockoff statistics

The default statistic suggested in Barber and Candès, () is the Lasso
Signed Max (LSM), which corresponds to $Z_{j}$ being the largest penalty
parameter at which the $j$th variable enters the model in the Lasso
regression of $y$ on $[\bm{X},\tilde{\bm{X}}]$, and
$f_{j}={\operatorname*{sign}(|Z_{j}|-|\tilde{Z}_{j}|)\max\{|Z_{j}|,|\tilde{Z}_{j}|\}}$. In addition to the
LSM statistic, Barber and Candès, () suggested alternatives such as the
difference in absolute values of estimated coefficients for
a variable and its knockoff:

$$ $W_{j}=|\hat{b}_{j}|-|\hat{b}_{j+p}|,$ $$

where the $\hat{b}_{j},\hat{b}_{j+p}$ are estimated so that $W$ obeys
the sufficiency property required by the FX knockoff procedure,
e.g., by ordinary least squares or the Lasso with a pre-specified tuning
parameter. The removal of the sufficiency requirement for MX
knockoffs allows us to improve this class of statistics by adaptively
tuning the fitted model. The simplest example is the LCD statistic
introduced in Section [3.2](#S3.SS2), which uses cross-validation to choose
the tuning parameter in the Lasso. Note the LCD statistic can be
easily extended to any GLM by replacing the first term in
([3.6](#S3.E6)) by a non-Gaussian negative log-likelihood, such as in
logistic regression; we will refer to all such statistics generically
as LCD. The key is that the tuning and cross-validation is done on the
augmented design matrix $[\bm{X},\,\tilde{\bm{X}}]$, so that $W$ still obeys the
flip-sign property.

More generally, MX knockoffs allows us to construct statistics that
are highly adaptive to the data, as long as that adaptivity does not
distinguish between original and knockoff variables. For instance, we could compute the
cross-validated error of the ordinary Lasso (still of $y$ on $[\bm{X},\tilde{\bm{X}}]$) and compare it to that of a random forest, and choose $Z$ to
be a feature importance measure derived from whichever one has smaller
error. Since the Lasso works best when the true model is close to
linear, while random forests work best in non-smooth models, this
approach gives us high-level adaptivity to the model smoothness, while
the MX knockoff framework ensures strict Type I error control.

Returning to the simpler example of adaptivity, we found the LCD
statistic to be uniformly more powerful than the LSM statistic across
a wide range of simulations (linear and binomial GLMs, ranging
covariate dependence, effect size, sparsity, sample size, total number
of variables), particularly under covariate dependence. We note,
however, the importance of choosing the penalty parameter that
minimizes the cross-validated error, as opposed to the default in some
computational packages of using the “one standard error” rule, as the latter
causes LCD to be underpowered compared to LSM in low-power
settings. Figure [2](#S5.F2) shows a simulation with $n=3000$,
$p=1000$ of a binomial linear model (with statistics computed from
Lasso logistic regression) that is representative of the power
difference between the two statistics. In all our simulations, unless otherwise specified,
MX knockoffs is always run using the LCD
statistic. Explicitly, when the response variable is continuous, we
use the standard Lasso with Gaussian linear model likelihood, and when
the response is binary, we use Lasso-penalized logistic regression.

Figure: Figure 2: Power and FDR (target is 10%) for knockoffs with the LCD and LSM statistics. The design matrix has i.i.d. rows and AR(1) columns with autocorrelation coefficient specified by the x-axes of the plots, and marginally each $X_{j}\sim\mathcal{N}(0,1/n)$. Here, $n=3000$, $p=1000$, and $y$ comes from a binomial linear model with logit link function with 60 nonzero regression coefficients of magnitude 3.5 and random signs. Each point represents 200 replications.
Refer to caption: /html/1610.02351/assets/x3.png

#### 5.2.2 Bayesian knockoff statistics

Another very interesting source of knockoff statistics comes from
Bayesian procedures. If a statistician has prior knowledge about the
problem, he or she can encode it in a Bayesian model and use the
resulting estimators to construct a statistic (e.g., difference of
absolute posterior mean coefficients, or difference or log ratio of
posterior probabilities of nonzero coefficients with a sparse
prior). What makes this especially appealing is that the statistician
gets the power advantages of incorporating prior information, while
maintaining a strict frequentist guarantee on the Type I error,
*even if the prior is false*!

As an example, we ran knockoffs in an experiment with a Bayesian
hierarchical regression model with $n=300$, $p=1000$, and
$\mbox{$\mathbb{E}\left(\|\beta\|_{0}\right)$}=60$ ($\|\cdot\|_{0}$ denotes the $\ell_{0}$ norm, or the
number of nonzero entries in a vector); see Appendix [C](#A3) for
details. We chose a simple canonical model with Gaussian response to
demonstrate our point, but the same principle applies to more complex,
nonlinear, and non-Gaussian Bayesian models as well. The statistics we
used were the LCD and a Bayesian variable selection (BVS) statistic,
namely $Z_{j}-\tilde{Z}_{j}$ where $Z_{j}$ and $\tilde{Z}_{j}$ are the
posterior probabilities that the $j$th original and knockoff
coefficients are nonzero, respectively (George and McCulloch,, ); again see
Appendix [C](#A3) for details. Figure [3](#S5.F3) shows
that the accurate prior information supplied to the Bayesian knockoff
statistic gives it improved power over LCD which lacks such
information, but that they have the same FDR control (and they would
even if the prior information were incorrect).

Figure: Figure 3: Power and FDR (target is 10%) for knockoffs with the LCD and BVS statistics. The design matrix is i.i.d. $\mathcal{N}(0,1/n)$, $n=300$, $p=1000$, and $y$ comes from a Gaussian linear model with $\beta$ and the noise variance randomly chosen (see Appendix [C](#A3) for the precise model). Here, the nonzero entries of $\beta$ are Gaussian with mean zero and standard deviation given on the x-axis; the expected number of nonzero components is 60; the expected variance of the noise is 1. Each point represents 200 replications.
Refer to caption: /html/1610.02351/assets/x5.png

### 5.3 Alternative procedures

To assess the relative power of knockoffs, we compare to a number of
alternatives in settings in which they are valid:

- 1.
The FX knockoff procedure with settings recommended in
Barber and Candès, (). This method can only be
applied in homoscedastic Gaussian linear regression when $n\geq p$.
- 2.
BHq applied to asymptotic GLM p-values. This method can only
be applied when $n\geq p$, and although for linear regression exact
p-values can be computed (when the MLE exists), for any other GLM
these p-values can be far from valid unless $n\gg p$, as shown in
Section [5.1](#S5.SS1).
- 3.
BHq applied to marginal test p-values. The correlation
between the response and each covariate is computed and compared to its null
distribution, which under certain Gaussian assumptions is
closed-form, but in general can at least be simulated exactly by
conditioning on $y$ and using the known marginal distribution of
$X_{j}$. Although these tests are valid for testing hypotheses of *marginal* independence
(regardless of $n$ and $p$), such hypotheses only agree with the
desired *conditional* independence hypotheses when the
covariates are exactly independent of one another.
- 4.
BHq applied to the p-values from the conditional
randomization test described in
Section [1.6](#S1.SS6) and Section [4](#S4).

Note that we are using knockoffs, not knockoffs+, in all simulations, and thus we are technically controlling a slightly modified version of FDR. The FDR is nevertheless effectively controlled in all simulations except in extremely low-power settings, and even then the violations are small. We could have sacrificed a small amount of power and used knockoffs+ (both MX and FX) for exact FDR control, but then a more fair comparison in settings 2–4 would replace BHq with the conservative procedure in Benjamini and Yekutieli, (), since the joint distribution of the p-values will not in general satisfy the assumptions for BHq to control the FDR exactly. However that conservative procedure had extremely noncompetitive power, so we prefer instead to compare knockoffs and regular BHq, which are more powerful and still effectively control the FDR.

#### 5.3.1 Comparison with conditional randomization

We start by comparing MX knockoffs with procedure 4, BHq
applied to conditional randomization test p-values, for computational
reasons. We simulated $n=400$ i.i.d. rows of $p=600$ AR(1) covariates
with autocorrelation 0.3, and response following a logistic regression
model with 40 nonzero coefficients of random
signs. Figure [4](#S5.F4) shows the power and FDR curves as the
coefficient amplitude was varied.

Figure: Figure 4: Power and FDR (target is 10%) for MX knockoffs and BHq applied to conditional randomization test p-values. The design matrix has i.i.d. rows and AR(1) columns with autocorrelation 0.3, $n=400$, $p=600$, and $y$ comes from a binomial linear model with logit link function with $\|\beta\|_{0}=40$, and all nonzero entries of $\beta$ having equal magnitudes and random signs; each point represents 200 replications.
Refer to caption: /html/1610.02351/assets/x7.png

We see that the conditional randomization test gives higher
power with similar FDR control, but this comes at a hugely
increased computational cost. This simulation has considerably
smaller $n$ and $p$ than any other simulation in the paper,
and we still had to apply a number of computational
speed-ups/shortcuts, described in Appendix [D](#A4), to keep the computation time within
reason.

With these speed-ups, Figure [4](#S5.F4) took roughly three
years of serial computation time, while the MX knockoffs
component took only about six hours, or about 5000 times less (all
computation was run in Matlab 2015b, and both methods used glmnet to
compute statistics). Because of the heavy computational burden, we
were unable to include the conditional randomization test in our
further, larger simulations—we show in Section [4.3](#S4.SS3) that
the number of $T_{j}$ computations scales optimistically linearly in
$p$. To summarize, conditional randomization testing appears somewhat
more powerful than MX knockoffs, but is computationally
infeasible for large data sets (like that in
Section [7](#S7)).

#### 5.3.2 Effect of signal amplitude

Our first simulation comparing MX knockoffs to procedures 1–3 is by
necessity in a Gaussian linear model with $n>p$ and independent
covariates—the only setting in which all procedures approximately
control the FDR. Specifically, the left side of
Figure [5](#S5.F5) plots the power and FDR for the four
procedures when
$X_{ij}\stackrel{{\scriptstyle\emph{i.i.d.}}}{{\sim}}\mathcal{N}(0,1/n)$, $n=3000$, $p=1000$,
$\|\beta\|_{0}=60$, the noise variance $\sigma^{2}=1$, and the nonzero
entries of $\beta$ have random signs and equal magnitudes, varied
along the x-axis. All methods indeed control the FDR, and MX
knockoffs is the most powerful, with as much as 10% higher power than
its *nearest* alternative. The right side of Figure [5](#S5.F5) shows the same setup but in high dimensions:
$p=6000$. In the high-dimensional regime, neither maximum likelihood
p-values nor FX knockoffs can even be computed, and the MX
knockoff procedure has considerably higher power than BHq applied
to marginal p-values.

Figure: Figure 5: Power and FDR (target is 10%) for MX knockoffs and alternative procedures. The design matrix is i.i.d. $\mathcal{N}(0,1/n)$, $n=3000$, $p=$ (left): $1000$ and (right): $6000$, and $y$ comes from a Gaussian linear model with 60 nonzero regression coefficients having equal magnitudes and random signs. The noise variance is 1. Each point represents 200 replications.
Refer to caption: /html/1610.02351/assets/x9.png

Next we move beyond the Gaussian linear model to a binomial linear
model with logit link function, precluding the use of the original
knockoff procedure. Figure [6](#S5.F6) shows the same simulations
as Figure [5](#S5.F5) but with $Y$ following the binomial
model. The results are similar to those for the Gaussian linear model, except
that BHq applied to the asymptotic maximum likelihood p-values now has
an FDR above 50% (rendering its high power meaningless), which can be
understood as a manifestation of the phenomenon from Section [5.1](#S5.SS1). In
summary, MX knockoffs continues to have the highest power among
FDR-controlling procedures.

Figure: Figure 6: Power and FDR (target is 10%) for MX knockoffs and alternative procedures. The design matrix is i.i.d. $\mathcal{N}(0,1/n)$, $n=3000$, $p=$ (left): $1000$ and (right): $6000$, and $y$ comes from a binomial linear model with logit link function, and 60 nonzero regression coefficients having equal magnitudes and random signs. Each point represents 200 replications.
Refer to caption: /html/1610.02351/assets/x11.png

#### 5.3.3 Effect of covariate dependence

To assess the relative power and FDR control of MX knockoffs
as a function of covariate dependence, we ran similar simulations as
in the previous section, but with covariates that are AR(1) with varying
autocorrelation coefficient (while the coefficient amplitude remains
fixed). It is now relevant to specify that the locations of the
nonzero coefficients are uniformly distributed on $\{1,\dots,p\}$. In
the interest of space, we only show the low-dimensional ($p=1000$)
Gaussian setting (where all four procedures can be computed) and the
high-dimensional ($p=6000$) binomial setting, as little new
information is contained in the plots for the remaining two
settings. Figure [7](#S5.F7) shows that, as
expected, BHq with marginal testing quickly loses FDR control with
increasing covariate dependence. This is because the marginal tests
are testing the null hypothesis of *marginal* independence
between covariate and response, while recall from
Definition [2.1](#S2.Thmtheorem1) that all conditionally independent
covariates are considered null, even if they are marginally dependent
on the response. Concentrating on the remaining methods and just the
left-hand part of the BHq marginal curves where FDR is controlled,
Figure [7](#S5.F7) shows that MX knockoffs continues to be considerably more
powerful than alternatives as covariate dependence is introduced, in
low- and high-dimensional linear and nonlinear models.

Figure: Figure 7: Power and FDR (target is 10%) for MX knockoffs and alternative procedures. The design matrix has i.i.d. rows and AR(1) columns with autocorrelation coefficient specified by the x-axes of the plots, and marginally each $X_{j}\sim\mathcal{N}(0,1/n)$. (left): $n=3000$, $p=1000$, and $y$ follows a Gaussian linear model. (right): $n=3000$, $p=6000$, and $y$ follows a binomial linear model with logit link function. In both cases, there are 60 nonzero coefficients having magnitudes equal to 3.5 on the left and 10 on the right, random signs, and randomly selected locations. Each point represents 200 replications.
Refer to caption: /html/1610.02351/assets/x13.png

## 6 Robustness

In many real applications, the true joint covariate distribution may
not be known exactly, forcing the user to estimate it from the
available data. As already mentioned, this is a challenging problem by
itself, but often we have considerable outside information or
unsupervised data that can be brought to bear to improve
estimation. This raises the important question of how robust
MX knockoffs is to error in the joint covariate
distribution. Theoretical guarantees of robustness are
beyond the scope of this paper, but we present instead three
compelling simulation studies to demonstrate robustness. The
first study investigates error that biases that distribution toward the empirical
covariate distribution, often referred to as overfitting
error, on simulated data. We generated knockoffs for Gaussian variables,
but instead of using the true covariance matrix, we used
in-sample covariance estimates which ranged in
overfitting error. Figure [8](#S6.F8) shows the power and FDR as the
covariance we use ranges from the true covariance matrix (AR(1) with
autocorrelation 0.3), to a graphical lasso estimator, to
convex combinations of the true and empirical covariance (see
Appendix [E](#A5) for explicit formulae for the estimators). The
plot is indexed on the x-axis by the average relative
Frobenius norm $\|\hat{\bm{\Sigma}}-\bm{\Sigma}\|_{\text{Fro}}/\|\bm{\Sigma}\|_{\text{Fro}}$
of the estimator $\hat{\bm{\Sigma}}$.

Figure: Figure 8: Power and FDR (target is 10%) for knockoffs with the LCD statistic as the covariance matrix used to generate knockoffs ranges from the truth to an estimated covariance; see text for details. The design matrix has i.i.d. rows and AR(1) columns with autocorrelation coefficient 0.3, and the matrix (including knockoffs) is standardized so that each column has mean zero and Euclidean norm 1. Here, $n=800$, $p=1500$, and $y$ comes from a binomial linear model with logit link function with 50 nonzero entries having magnitude 20 and random signs. Each point represents 200 replications.
Refer to caption: /html/1610.02351/assets/x15.png

Although the graphical Lasso is well suited for this problem
since the covariates have a sparse precision matrix, its
covariance estimate is still off by nearly 50%, and yet
surprisingly the resulting power and FDR are nearly
indistinguishable from when the exact covariance is used. The
covariance estimate worsens as the empirical covariance—a
very poor estimate of the true covariance given the high
dimensionality—is combined in increasing proportion with the
truth. At 75% weight on the empirical covariance, the
covariance estimate is nearly 100% off and yet the power and FDR of MX
knockoffs are only slightly decreased. Beyond this point, MX
knockoffs becomes quite conservative, with power and FDR
approaching zero as the estimated covariance approaches the
empirical covariance. This behavior at 100% weight on the
empirical covariance is not
surprising, since $p>n$ and thus the empirical covariance is
rank-deficient, forcing the knockoff variables to be exact
replicas of their original counterparts.(^8^88When the knockoff variables are exact copies of the original variables we are guaranteed zero power and zero FDR since all $W_{j}=0$. Although in principle we could break ties and assign signs by coin flips when $W_{j}=0$, we prefer to only select $X_{j}$ with $W_{j}>0$, as $W_{j}=0$ provides no evidence against the null hypothesis.) The main take-aways from this plot are that (a)
the nominal level of 10% FDR is never violated, even for
covariance estimates very far from the truth, and (b) the
more overfitting done on the covariance, the more conservative
the procedure is, although even at almost 100% relative
error, MX knockoffs had only lost about 20% of the power it
would have had if the covariance were known exactly. Intuitively, instead of treating the variables as coming from their true joint distribution, MX knockoffs with an overfit covariate distribution seems to treat them as coming from their true distribution “conditional” on being similar to their observed values. Thus FDR should be roughly controlled *conditionally* which implies marginal FDR control, while power may be lost if the conditioning is too great, which matches what we see in the simulations.

Our second and third experiments use real covariate data from
a genome-wide association study (GWAS), the details of which
are given in Section [7](#S7) and
Appendix [F](#A6). In brief, it is a
high-dimensional setting with $X_{ij}\in\{0,1,2\}$ and strong
spatial structure, whose covariance we estimate in-sample
using the GWAS-tailored covariance estimator of
Wen and Stephens, (). We check the robustness of constructing
second-order model-X knockoffs by ASDP (from
Section [3.4.2](#S3.SS4.SSS2)) by choosing a reasonable but
artificial model for $Y|X_{1},\dots,X_{p}$ and simulating
artificial response data using the real covariate data. The
exact details of the simulation are given in
Section [7.1](#S7.SS1), but note that this simulation
used the exact same covariance estimation, SNP clustering and
representative selection, knockoff construction, and knockoff
selection procedure as used for the real data analysis of the
next section. Our second experiment varies the signal
amplitude in a binomial linear model, and
Figure [9](#S6.F9) shows the FDR and power. As hoped,
the FDR is consistently controlled over a wide range of
powers.

Figure: Figure 9: Power and FDR (target is 10%) for knockoffs with the LCD statistic applied to subsamples of a real genetic design matrix. Each boxplot represents 10 different logistic regression models with 60 nonzero coefficients with amplitudes given by the x-axis, and for each model, 1000 common observations were used for picking cluster representatives, and the remaining 13,708 observations were divided into 10 disjoint parts and knockoffs run on each part, with power and FDR for that model then computed by averaging the results over those 10 parts.
Refer to caption: /html/1610.02351/assets/x17.png

Our third experiment, instead of varying the signal strength
of our artificial model for $Y|X_{1},\dots,X_{p}$, deliberately
corrupts the covariance estimate of Wen and Stephens, () by
varying a shrinkage parameter that is not meant to be
varied. That parameter is $m$ in Wen and Stephens, (, equation
(2.7)), and we vary it from $1/10$ to $10$ times
its intended value. This variation has a huge effect on how
much shrinkage is applied off the diagonal, with the average
correlation varying by a factor of about 13 over the range of
shrinkage. Figure [10](#S6.F10) shows that, even
as we range from substantial undershrinkage to substantial
overshrinkage, MX knockoffs never significantly violates FDR
control, with only a bit of conservativeness when the
undershrinkage is most drastic (the same phenomenon as the
right side of Figure [8](#S6.F8)).

Figure: Figure 10: The setup is the same as Figure [9](#S6.F9) except with the amplitude fixed at 14 and the amount of off-diagonal shrinkage in the covariance estimate varied on the x-axis.
Refer to caption: /html/1610.02351/assets/x19.png

## 7 Genetic analysis of Crohn’s disease

To test the robustness and practicality of the new knockoff
procedure, we applied it to a data set containing genetic information
on cases and controls for Crohn’s disease (CD). The data is provided
by the Wellcomme Trust Case Control Consortium and has been studied
previously (WTCCC,, ). After preprocessing, there were $p=377,749$ single nucleotide
polymorphisms (SNPs) measured on $n=4,913$ subjects (1,917 CD patients
and 2,996 healthy controls). Although $p\gg n$, the joint dependence
of SNPs has a strong spatial structure, and outside data can be used
to improve estimation. In particular, we approximated the standardized
joint distribution as multivariate Gaussian with covariance matrix
estimated using the methodology of Wen and Stephens, (), which shrinks
the off-diagonal entries of the empirical covariance matrix using
genetic distance information estimated from the HapMap CEU
population. This approximation was used on each chromosome, and SNPs
on different chromosomes were assumed to be independent. The statistic
we use is the LCD. Although the data itself cannot be made available, all
code is available at [http://web.stanford.edu/group/candes/knockoffs/software/knockoff/](http://web.stanford.edu/group/candes/knockoffs/software/knockoff/).

One aspect of SNP data is that it contains some very high
correlations, which presents two challenges to our methodology. The
first is generic to the variable selection problem: it is very hard to
choose between two or more nearly-identical (highly-correlated)
variables if the data supports at least one of them being
selected.(^9^99This is purely a problem of power and would not affect the Type I error control of knockoffs.) To alleviate this, we
clustered the SNPs using the estimated correlations as a similarity
measure with a single-linkage cutoff of 0.5, and settle for
discovering important SNP clusters. To do so we choose one
representative from each cluster and approximate the null hypothesis
that a cluster is conditionally independent of the response given the
other clusters by the null hypothesis that a cluster
*representative* is conditionally independent of the response
given the other cluster *representatives*. To choose the
representatives, we could ignore the observed data altogether and do
something like pick representatives with the highest minor allele
frequency (computed from outside data), and then run knockoffs as
described in the paper. Although this produces a powerful procedure
(about 50% more powerful than the original analysis by
WTCCC, ()), a more powerful approach is to select cluster
representatives using a fraction of the observations, including their
responses, such as by marginal testing. Note that such a
data-splitting approach appears to make our null hypotheses random, as
in the work on inference after selection reviewed in
Section [1.4](#S1.SS4). However, the approximation we are making is
that each representative stands for its cluster, and each cluster has
exactly one associated null hypothesis, no matter how selection is
performed, even if it were nonrandom. That is, the *approximate*
hypotheses being tested do not actually depend on the selection
(unlike Brzyski et al., () where clusters are selected and where the
very definition of a cluster actually depends on the selection), and
our approach remains free of a model for $Y\mid X$, which together should make it clear
that it is still quite different from the literature on inference
after selection.

Explicitly, we randomly chose 20% of our observations and on those
observations only, we ran marginal t-tests between each SNP and the
response. Then from each cluster we chose the SNP with smallest t-test
p-value to be the single representative of that cluster. Because the
observations used for selecting cluster representatives have had their
representative covariate values selected for dependence on the
outcome, if we constructed knockoff variables as usual and included
them in our procedure, the required exchangeability established in
Lemma [3.2](#S3.Thmtheorem2) would be violated. However, taking a page from
Barber and Candès, (), we can still use these observations in our
procedure by making their knockoff variables just identical copies of
the original variables (just for these observations). It is easy to
show (see Appendix [G](#A7)) that constructing knockoffs in
this way, as exact copies for the observations used to pick
representatives and as usual for the remaining observations, the
pairwise exchangeability of null covariates with their knockoffs is
maintained. Of course, the observations with identical original and
knockoff covariate values do not directly help the procedure
distinguish between the original and knockoff variables, but including
them improves the accuracy of the fitting procedure used to produce
the feature importance statistics, so that power is improved
indirectly because the $Z_{j}$ become more accurate measures of feature
importance.

Replacing clusters with single representatives
reduces $p$ down to $71,145$ (so the average cluster size was just
over five SNPs, although there was substantial variance) and, by
construction, upper-bounds pairwise SNP correlations by 0.5. Note that
this is far from removing all dependence among the SNPs, so
considering conditional independence instead of marginal dependence
remains necessary for interpretation. Scientifically, we consider a
selected SNP to be a true discovery if and only if it is the
representative of a cluster containing a truly important SNP.

The second challenge is the one discussed in
Section [3.4](#S3.SS4), and we use the approximate SDP
knockoff construction proposed there. The approximate covariance
matrix was just the estimated covariance matrix with zeros outside of
the block diagonal, with the blocks chosen by single-linkage
clustering on the estimated correlation matrix, aggregating clusters
up the dendrogram as much as possible subject to a maximum cluster
size of 999. In this case, even separating the problem by chromosome,
the SDP construction was computationally infeasible and the
equicorrelated construction produced extremely small $s$:
$\text{mean}(s^{\text{EQ}})=0.08$.
The parallelized approximate SDP construction took just a matter of
hours to run, and increased $s$ on average by almost an order of magnitude,
with $\text{mean}(s^{\text{ASDP}})=0.57$.

Although it incorporates strong prior information, our estimate of the
joint distribution of the SNPs is still an approximation, so our
first step is to test the robustness of knockoffs to this
approximation.

### 7.1 Simulations with genetic design matrix

In Section [6](#S6), the second and third
experiments check the robustness of knockoffs to this
particular joint distribution approximation by taking a
reasonable model for $Y|X_{1},\dots,X_{p}$ and simulating
artificial response data using the real covariate data
itself. We split the rows of our design matrix into 10 smaller
data sets (re-estimating the joint covariate distribution each
time), and for each conditional model we run knockoffs 10
times and compute the realized FDP each time. Then averaging
the 10 FDPs gives an estimate of knockoffs’ FDR for that
conditional model. In an attempt to make each smaller data set
have size more comparable to the $n\approx 5,000$ in our
actual experiment, we combined the healthy and CD genetic
information with that from 5 other diseases from the same data
set:(^10^1010Bipolar disorder was also part of this data set, but a formatting error in the data we were given prevented us from including it.) coronary artery disease,
hypertension, rheumatoid arthritis, type 1 diabetes, and type
2 diabetes. This made for $14,708$ samples altogether. To
further increase the effective sample size and match the
actual experiment, we used a random subset of 1000
observations for choosing cluster representatives, but the
*same* 1000 for each smaller data set, so that each
subsampled data set contained $\approx 1,400$ unique samples,
$+1000$ common samples (for each of these common samples
$\tilde{X}_{ij}=X_{ij}$). For computational reasons, we used
only the first chromosome in this experiment, so each of our
simulations had (pre-clustering) $29,258$ covariates. The
conditional model for the response was chosen to be a logistic
regression model with 60 nonzero coefficients of random signs
and locations uniformly chosen from among the original (not
just representatives) SNPs. In the second experiment (Figure [9](#S6.F9)), the coefficient amplitude was
varied, and for each amplitude value, 10 different conditional
models (different random locations of the nonzero
coefficients) were simulated. Each simulation ran the exact
same covariance estimation, SNP clustering and representative
selection, knockoff construction, and knockoff selection
procedure as used for the real data. In the third
experiment (Figure [10](#S6.F10)), the
shrinkage was amplified or attenuated by varying the value
$m$ in Wen and Stephens, (, Equation (2.7)) from its intended
value, while the coefficient amplitude was fixed at 14. In both of these experiments, knockoffs appears to be quite robust to the approximation we make in estimating the joint distribution of the SNPs.

### 7.2 Results on real data

Encouraged by these robustness simulations, we
proceeded to run knockoffs with a nominal FDR level of 10% on the
full genetic design matrix and real Crohn’s disease outcomes. Since
knockoffs is a randomized procedure, we re-ran knockoffs 10 times
(after choosing the representatives)
and recorded the selected SNPs over all repetitions, summarized in
Table [2](#S7.T2). The serial computation time for a single
run of knockoffs was about 6 hours, but the knockoff generation
process is trivially parallelizable over chromosomes, so with 20
available computation nodes, the total parallelized computation time was
about one hour.

**Table 2: SNP clusters discovered to be important for Crohn’s disease over 10 repetitions of knockoffs. Clusters not found in Franke et al., () represent promising sites for further investigation, especially rs6601764 and rs4692386, whose nearest genes have been independendently linked to CD. See text for detailed description. SNP positions are as listed in the original data, which uses Human Genome Build 35.**
| Selection frequency | Cluster    Representative (Cluster Size) | Chrom. | Position Range (Mb) | Confirmed in Franke et al., ()? | Selected in WTCCC, ()? |
| --- | --- | --- | --- | --- | --- |
| 100% | rs11805303 (16) | 1 | 67.31–67.46 | Yes | Yes |
| 100% | rs11209026 (2) | 1 | 67.31–67.42 | Yes | Yes |
| 100% | rs6431654 (20) | 2 | 233.94–234.11 | Yes | Yes |
| 100% | rs6601764 (1) | 10 | 3.85–3.85 | No | No |
| 100% | rs7095491 (18) | 10 | 101.26–101.32 | Yes | Yes |
| 90% | rs6688532 (33) | 1 | 169.4–169.65 | Yes | No |
| 90% | rs17234657 (1) | 5 | 40.44–40.44 | Yes | Yes |
| 90% | rs3135503 (16) | 16 | 49.28–49.36 | Yes | Yes |
| 80% | rs9783122 (234) | 10 | 106.43–107.61 | No | No |
| 80% | rs11627513 (7) | 14 | 96.61–96.63 | No | No |
| 60% | rs4437159 (4) | 3 | 84.8–84.81 | No | No |
| 60% | rs7768538 (1145) | 6 | 25.19–32.91 | Yes | No |
| 60% | rs6500315 (4) | 16 | 49.03–49.07 | Yes | Yes |
| 60% | rs2738758 (5) | 20 | 61.71–61.82 | Yes | No |
| 50% | rs7726744 (46) | 5 | 40.35–40.71 | Yes | Yes |
| 50% | rs4246045 (46) | 5 | 150.07–150.41 | Yes | Yes |
| 50% | rs2390248 (13) | 7 | 19.8–19.89 | No | No |
| 50% | rs7186163 (6) | 16 | 49.2–49.25 | Yes | Yes |
| 40% | rs10916631 (14) | 1 | 220.87–221.08 | No | No |
| 40% | rs4692386 (1) | 4 | 25.81–25.81 | No | No |
| 40% | rs7655059 (5) | 4 | 89.5–89.53 | No | No |
| 40% | rs7759649 (2) | 6 | 21.57–21.58 | Yes* | No |
| 40% | rs1345022 (44) | 9 | 21.67–21.92 | No | No |
| 30% | rs6825958 (3) | 4 | 55.73–55.77 | No | No |
| 30% | rs9469615 (2) | 6 | 33.91–33.92 | Yes* | No |
| 30% | rs4263839 (23) | 9 | 114.58–114.78 | Yes | No |
| 30% | rs2836753 (5) | 21 | 39.21–39.23 | No | No |
| 10% | rs459160 (2) | 1 | 44.75–44.75 | No | No |
| 10% | rs6743984 (23) | 2 | 230.91–231.05 | Yes | No |
| 10% | rs2279980 (20) | 5 | 57.95–58.07 | No | No |
| 10% | rs4959830 (11) | 6 | 3.36–3.41 | Yes | No |
| 10% | rs13230911 (9) | 7 | 1.9–2.06 | No | No |
| 10% | rs7807268 (5) | 7 | 147.65–147.7 | No | No |
| 10% | rs2147240 (1) | 9 | 71.83–71.83 | No | No |
| 10% | rs10761659 (53) | 10 | 64.06–64.41 | Yes | Yes |
| 10% | rs4984405 (3) | 15 | 93.06–93.08 | No | No |
| 10% | rs17694108 (1) | 19 | 38.42–38.42 | Yes | No |
| 10% | rs3932489 (30) | 20 | 15.01–15.09 | No | No |

Although in this case we certainly do not
know the ground truth, we can get some sort of confirmation by
comparing to the results of studies with newer and much larger data
sets than ours. In particular, we compared with the results of
Franke et al., (), which used roughly 22,000 cases and 29,000
controls, or about 10 times the sample size of the WTCCC data. We also
compare to the WTCCC, () results, where the p-value cutoff
used was justified as controlling the Bayesian FDR at close to 10%—the same level we use. We
consider discovered clusters in different studies to correspond
(“Yes” in Table [2](#S7.T2)) if
their position ranges overlap, and to nearly correspond (“Yes*” in
Table [2](#S7.T2)) if the distance from our discovered cluster to the
nearest cluster in the other study was less than the width of that
cluster in the other study.

One thing to notice in the table is that a small number of the discovered
clusters actually overlap with other clusters, specifically the
clusters represented by rs11805303 and rs11209026 on chromosome 1, and
rs17234657 and rs7726744 on chromosome 5. And although they don’t
overlap one another, the three nearby clusters represented by rs3135503,
rs6500315, and rs7186163 on chromosome 16 all overlap the same
discovered region in Franke et al., (). Although puzzling at first,
this phenomenon is readily explained by one of four possibilities:

- •
By construction, the representatives of overlapping clusters are
not very highly-correlated (less than 0.5), so the fact that
knockoffs chose multiple clusters in the same region may mean there
are multiple important SNPs in this region, with one (or more) in
each cluster. Focusing on the clusters on chromosome 1, the same
region on the IL23R gene was reported in Franke et al., () to have
by far the strongest signal (estimated odds ratio of 2.66, next
highest was 1.53) among all the regions they identified. If we
conclude from this that the region or gene is fundamentally
important for Crohn’s disease, it stands to reason that mutations at
multiple nearby but distinct loci could have important detrimental
effects of their own.
- •
There could be an important SNP located between the two
clusters, but which was not in the data. Then the two clusters on either
side would both be conditionally dependent on the
response, and knockoffs would be correct to reject them.
- •
One or more of the overlapping clusters could be mundane false
discoveries caused by null covariates that happened to take values
more conditionally related to the response than would be
typical. This would be a facet of the data itself, and thus an
unavoidable mistake.
- •
The tandem discoveries could also be due to a breakdown in our
covariate distribution estimate. If the covariance between the two
representatives were substantially underestimated and one had a
large effect while the other was null, then the Lasso would have a
(relatively) hard time distinguishing the two original variables,
but a much easier time separating them from the knockoff of the null
representative, since it is less correlated with the signal
variable. As a result, the null variable and its knockoff would not
be exchangeable as required by knockoffs, and a consistent error
could be made. However, given that the empirical and estimated
correlations between the two representatives on chromosome 1 were
-0.1813 and -0.1811, respectively, and for the two on chromosome 5
were -0.2287 and -0.2286, respectively, it seems unlikely that we
have made a gross estimation error. Also, note that the separation
of nearly all the discovered regions, along with the simulations of
the previous subsection, suggest this effect is at worst very small
among our discoveries.

Overlapping clusters aside, the knockoffs results display a number of
advantages over the original marginal analysis in WTCCC, ():

- •
First, the power is much higher, with WTCCC, ()
making 9 discoveries, while knockoffs made 18 discoveries on
average, doubling the power.
- •
Quite a few of the discoveries made by knockoffs that were confirmed by
Franke et al., () were not discovered in WTCCC, ()’s
original analysis.
- •
Knockoffs made a number of discoveries not found in either
WTCCC, () or Franke et al., (). Of course we expect some
(roughly 10%) of these to be false discoveries, particularly towards
the bottom of the table. However, especially given the evidence from
the simulations of the previous subsection suggesting the FDR is
controlled, it is likely that many of these correspond to true
discoveries. Indeed, evidence from independent studies about
adjacent genes shows some of the top hits to be promising
candidates. For example, the closest gene to rs6601764 is KLF6,
which has been found to be associated with multiple forms of IBD,
including CD and ulcerative colitis (Goodman et al.,, ); and the closest
gene to rs4692386 is RBP-J, which has been linked to CD through its
role in macrophage polarization (Barros et al.,, ).

Note that these benefits required relatively little customization of
the knockoff procedure. For instance, WTCCC, () used marginal tests specifically
tailored to SNP case-control data, while we simply used the LCD
statistic. We conjecture that the careful use of knockoffs by domain
experts would compound the advantages of knockoffs, as such
users could devise more powerful statistics and better model/cluster the
covariates for their particular application.

## 8 Discussion

This paper introduced a novel approach to variable selection
in general non-parametric models, which teases apart important
from irrelevant variables while guaranteeing Type I error
control. This approach is a significant rethinking of the
knockoff filter from Barber and Candès, (). A distinctive feature
of our approach is that selection is achieved without ever
constructing p-values. This is attractive since (1) p-values
are not needed and (2) it is unclear how they could be
efficiently constructed, in general. (The conditional
randomization approach we proposed is one way of getting such
p-values but it comes at a computational cost.)

##### Deployment in highly correlated settings

We posed a simple
question: which variables does a response of interest depend upon? In
many problems, there may not be enough “resolution” in the data to
tell whether $Y$ depends on $X_{1}$ or, instead, upon $X_{2}$ when the two
are strongly correlated. This issue is apparent in our genetic
analysis of Crohn’s disease from Section [7](#S7), where
co-located SNPs may be extremely correlated. In such examples,
controlling the FDR may not be a fruitful question. A more meaningful
question is whether the response appears to depend upon a group of
correlated variables while controlling for the effects of a number of
other variables (e.g., from SNPs located in a certain region of the
genome while controlling for the effects of SNPs elsewhere on the
chromosomes). In such problems, we envision applying our techniques to
grouped variables: one possibility is to develop a model-X group
knockoff approach following Dai and Barber, (). Another is to
construct group representatives and proceed as we have done in Section
[7](#S7). It is likely that there are several other ways to
formulate a meaningful problem and solution.

##### Open questions

Admittedly, this paper may pose more problems than it solves;
we close our discussion with a few of them below.

- •
How do we construct MX knockoffs? Even though we presented
a general strategy for constructing knockoffs, we have essentially
skirted this issue except for the important case of Gaussian
covariates. It would be important to address this problem, and write
down concrete algorithms for some specific distributions of features
of practical relevance.
- •
Which MX knockoffs? Even in the case of Gaussian
covariates, the question remains of how to choose
$\operatorname{corr}(X_{j},\tilde{X}_{j})$ or, equivalently, the
parameter $s_{j}$ from Section [3](#S3) since
$\operatorname{corr}(X_{j},\tilde{X}_{j})=1-s_{j}$. Should we make
the marginal correlations small? Should we make the partial
correlations small? Should we take an information-theoretic approach
and minimize a functional of the joint distribution such as the
mutual information between $X$ and $\tilde{X}$?
- •
What would we do with multiple MX knockoffs? As suggested
in Barber and Candès, (), we could in principle construct multiple
knockoff variables $(\tilde{X}^{(1)},\ldots,\tilde{X}^{(d)})$ in such a way
that the $(d+1)p$-dimensional family
$(X,\tilde{X}^{(1)},\ldots,\tilde{X}^{(d)})$ obeys the following extended
exchangeability property: for any variable $X_{j}$, any permutation in
the list $(X_{j},\tilde{X}_{j}^{(1)},\ldots,\tilde{X}_{j}^{(d)}$) leaves the
joint distribution invariant. On the one hand, such constructions
would yield more accurate information since we could compute, among
multiple knockoffs, the rank with which an original variable enters
into a model. On the other hand, this would constrain the
construction of knockoffs a bit more, perhaps making them less
distinguishable from the original features. What is the right trade
off?
Another point of view is to construct several knockoff matrices
exactly as described in the paper. Each knockoff matrix would yield
a selection, with each selection providing FDR control as described
in this paper. Now an important question is this: is it possible to
combine/aggregate all these selections leading to an increase in
power while still controlling the FDR?
- •
Can we prove some form of robustness? Although our
theoretical guarantees rely on the knowledge of the joint covariate
distribution, Section [6](#S6) showed preliminary examples with remarkable robustness when this
distribution is simply estimated from data. For instance, the
estimation of the precision matrix for certain Gaussian designs
seems to have rather secondary effects on FDR and power levels. It
would be interesting to provide some theoretical insights into this
phenomenon.
- •
Which feature importance statistics should we use? The
knockoff framework can be seen as an inference machine: the
statistician provides the test statistic $W_{j}$ and the machine
performs inference. It is of interest to understand which statistics
yield high power, as well as design new ones.
- •
Can we speed up the conditional randomization testing
procedure? Conditional randomization provides a powerful
alternative method for controlling the false discovery rate in
model-X variable selection, but at a computational cost that is
currently prohibitive for large problems. However, there exist a
number of promising directions for speeding it up, including:
importance sampling to estimate small p-values with fewer
randomizations, faster feature statistics $T_{j}$ with comparable or
higher power than the absolute value of lasso-estimated
coefficients, and efficient computation re-use and warm starts to
take advantage of the fact that each randomization changes
only a single column of the design matrix.

In conclusion, much remains to be done. On the upside, though, we have
shown how to select features in high-dimensional nonlinear models
(e.g., GLMs) in a reliable way. This arguably is a fundamental
problem, and it is really not clear how else it could be achieved.

## Appendix A Relationship with the knockoffs of Barber and Candès

We can immediately see a key difference with the earlier
framework of Barber and Candès, (). In their work, the design matrix is
viewed as being fixed and setting $\bm{\Sigma}=\bm{X}^{\top}\bm{X}$, knockoff
variables are constructed as to obey

$$ $[\bm{X},\,\tilde{\bm{X}}]^{\top}\,[\bm{X},\,\tilde{\bm{X}}]=\bm{G},\quad\bm{G}=\begin{bmatrix}\bm{\Sigma}&\bm{\Sigma}-\operatorname{diag}\{s\}\\ \bm{\Sigma}-\operatorname{diag}\{s\}&\bm{\Sigma}\end{bmatrix}.$ (A.1) $$

Imagine the columns of $\bm{X}$ are centered, i.e., have vanishing means,
so that we can think of $\bm{X}^{\top}\bm{X}$ as a sample covariance
matrix. Then the construction of Barber and Candès is asking that
the sample covariance matrix of the joint set of variables obeys
the exchangeability property—i.e., swapping rows and columns leaves
the covariance invariant—whereas in this paper, it is the population covariance that must be invariant. In particular, the
MX knockoffs will be far from obeying the relationship
$\bm{X}^{\top}\bm{X}=\tilde{\bm{X}}^{\top}\tilde{\bm{X}}$ required in
([A.1](#A1.E1)). To drive this point home, assume
$X\sim\mathcal{N}(0,\bm{I})$. Then we can choose
$\tilde{X}\sim\mathcal{N}(0,\bm{I})$ independently from $X$, in
which case $\bm{X}^{\top}\bm{X}/n$ and $\tilde{\bm{X}}^{\top}\tilde{\bm{X}}/n$ are
two independent Wishart variables. An important consequence of this is
that in the MX approach, the sample correlation between the $j$th
columns of $\bm{X}$ and $\tilde{\bm{X}}$ will typically be far smaller than
that in the original framework of Barber and Candès. For example,
take $n=3000$ and $p=1000$ and assume the equicorrelated
construction from Barber and Candès, (). Then the sample correlation
between any variable $\bm{X}_{j}$ and its knockoff $\tilde{\bm{X}}_{j}$ will
be about 0.65 while in the random case, the average magnitude of the
correlation is about 0.015. This explains the power gain of our
new method that is observed in Section [5](#S5).

## Appendix B Sequential conditional independent pairs algorithm

We will prove by induction on $j$ that Algorithm [1](#algorithm1)
produces knockoffs that satisfy the exchangeability property
([3.1](#S3.E1)). We prove the result for the discrete case; the
general case follows the same argument with a slightly more careful
measure-theoretic treatment using Radon–Nikodym derivatives instead
of probability mass functions. Below we denote the probability mass
function (PMF) of $(X_{1:p},\tilde{X}_{1:j-1})$ by
$\mathcal{L}(X_{\text{-}j},X_{j},\tilde{X}_{1:j-1})$.

###### Induction Hypothesis 1 .

After $j$ steps, every pair $X_{k},\tilde{X}_{k}$ is exchangeable in the
joint distribution of $(X_{1:p},\tilde{X}_{1:j})$ for $k=1,\ldots,j$.

By construction, the induction hypothesis is true after 1 step since
$X_{1}$ and $\tilde{X}_{1}$ are conditionally independent and have the same
marginal distribution (this implies conditional exchangeability).
Assuming the induction hypothesis holds until $j-1$, we prove that it
holds after $j$ steps. Note that by by assumption, $\mathcal{L}$ is
symmetric in $X_{k},\tilde{X}_{k}$ (that is, the function value remains
unchanged when the argument values $X_{k},\tilde{X}_{k}$ are swapped) for
$k=1,\dots,j-1$. Also, the conditional PMF of $\tilde{X}_{j}$ given
$X_{1:p},\tilde{X}_{1:j-1}$ is given by

$$ $\frac{\mathcal{L}(X_{\text{-}j},\tilde{X}_{j},\tilde{X}_{1:j-1})}{\sum_{u}\mathcal{L}(X_{\text{-}j},u,\tilde{X}_{1:j-1})}.$ $$

Therefore, the joint PMF of $(X_{1:p},\tilde{X}_{1:j})$ is given by the
product of the aforementioned conditional PMF with the joint PMF of $(X_{1:p},\tilde{X}_{1:j-1})$:

$$ $\frac{\mathcal{L}(X_{\text{-}j},X_{j},\tilde{X}_{1:j-1})\mathcal{L}(X_{\text{-}j},\tilde{X}_{j},\tilde{X}_{1:j-1})}{\sum_{u}\mathcal{L}(X_{\text{-}j},u,\tilde{X}_{1:j-1})}.$ (B.1) $$

Exchangeability of $X_{j},\tilde{X}_{j}$ follows from the symmetry of
([B.1](#A2.E1)) to those two values. For $k<j$, note that
([B.1](#A2.E1)) only depends on $X_{k},\tilde{X}_{k}$ through the function
$\mathcal{L}$, and that $\mathcal{L}$ is symmetric in $X_{k},\tilde{X}_{k}$. Therefore,
([B.1](#A2.E1)) is also symmetric in $X_{k},\tilde{X}_{k}$, and therefore
the pair is exchangeable in the joint distribution of
$(X_{1:p},\tilde{X}_{1:j})$.

###### Induction Hypothesis 1 .

## Appendix C Bayesian knockoff statistics

The data for the simulation of Section [5.2.2](#S5.SS2.SSS2) was drawn from:

$$ $\displaystyle X_{j}$ $\displaystyle\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathcal{N}(0,1/n),\;j\in\{1,\dots,p\},$ $\displaystyle\beta_{j}$ $\displaystyle\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathcal{N}(0,\tau^{2}),\;j\in\{1,\dots,p\},$ $\displaystyle\delta_{j}$ $\displaystyle\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\text{Bernoulli}(\pi),\;j\in\{1,\dots,p\},$ $\displaystyle\frac{1}{\sigma^{2}}$ $\displaystyle\sim\text{Gamma}(A,B)\qquad\text{(shape/scale parameterization, as opposed to shape/rate)},$ $\displaystyle Y$ $\displaystyle\sim\mathcal{N}\left(\sum_{j:\delta_{j}=1}X_{j}\beta_{j},\;\sigma^{2}\right).$ $$

The simulation used $n=300$, $p=1000$, with parameter values $\pi=\frac{60}{1000}$, $A=5$, $B=4$, and $\tau$ varied along the x-axis of the plot.

To compute the Bayesian variable selection (BVS) knockoff statistic,
we used a Gibbs sampler on the following model (treating $X_{1},\dots,X_{p}$ and $\tilde{X}_{1},\dots,\tilde{X}_{p}$ as fixed):

$$ $\displaystyle\beta_{j}$ $\displaystyle\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathcal{N}(0,\tau^{2}),\;j\in\{1,\dots,p\},$ $\displaystyle\tilde{\beta}_{j}$ $\displaystyle\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathcal{N}(0,\tau^{2}),\;j\in\{1,\dots,p\},$ $\displaystyle\lambda_{j}$ $\displaystyle\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\text{Bernoulli}(\pi),\;j\in\{1,\dots,p\},$ $\displaystyle(\delta_{j},\tilde{\delta}_{j})$ $\displaystyle\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\left\{\begin{array}[]{ll}(0,0)&\text{if }\lambda_{j}=0\\ (0,1)\text{ w.p. }1/2&\text{if }\lambda_{j}=1\\ (1,0)\text{ w.p. }1/2&\text{if }\lambda_{j}=1\end{array}\right\},\;j\in\{1,\dots,p\},$ $\displaystyle\frac{1}{\sigma^{2}}$ $\displaystyle\sim\text{Gamma}(A,B)\qquad\text{(shape/scale parameterization, as opposed to shape/rate)},$ $\displaystyle Y$ $\displaystyle\sim\mathcal{N}\left(\sum_{j:\delta_{j}=1}X_{j}\beta_{j}+\sum_{j:\tilde{\delta}_{j}=1}\tilde{X}_{j}\tilde{\beta}_{j},\;\sigma^{2}\right),$ $$

which requires only a very slight modification of the procedure in
George and McCulloch, (). After computing the posterior probabilities
$\hat{\delta}_{j}$ and $\hat{\tilde{\delta}}_{j}$ with 500 Gibbs samples (after 50 burn-in samples), we
computed the $j$th knockoff statistic as

$$ $W_{j}=\hat{\delta}_{j}-\hat{\tilde{\delta}}_{j}.$ $$

## Appendix D Conditional randomization speedups

In order to make the conditional randomization computation feasible
for Figure [4](#S5.F4), we had to apply a number of
speed-ups/shortcuts. With these, the serial computation time for the figure was
reduced to roughly three years.

- •
As mentioned, the LCD statistic is a powerful one for MX
knockoffs, so we wanted to compare it to the analogue for the
conditional randomization test, namely, the absolute value of the
Lasso-estimated coefficient. However, choosing the penalty parameter
by cross-validation turned out to be too expensive, so instead we
chose a fixed value by simulating repeatedly from the known model
(with amplitude 30), running cross-validation for each
repetition, and choosing the median of the
cross-validation-error-minimizing penalty parameters (the chosen
value was 0.00053). For a fair comparison, we did the same for
MX knockoffs (except the simulations included knockoff
variables too—the chosen value was 0.00077). This speed-up is of
course impossible in practice, as the true model is not known. It
is not clear how one would choose the penalty parameter if not by
cross-validation (at considerable extra computational
expense), although a topic of current research is how to choose tuning
parameters efficiently and without explicit reliance on a model for
the response, e.g., (Lv and Liu,, ).
- •
Because the statistic used was the (absolute value of the)
estimated coefficient from a sparse estimation procedure, many of
the observed statistics were exactly zero, and for these, the
conditional randomization p-values can be set to one without further
computation. This did not require any prior knowledge, although it
will only work for statistics whose distribution has a point mass at
the smallest possible value.
- •
Because the power was calibrated, we knew to expect at least
around 10 discoveries, and thus could anticipate the BHq cutoff
being at least $0.1\times 10/600$. This cutoff gives a sense of
the p-value resolution needed, and we chose the number of
randomizations to be roughly 10 times the inverse of the BHq
cutoff, namely, $10,000$. However, we made sure that all $10,000$
randomizations were only computed for very few covariates, both
using the previous bullet point and also by checking after
periodic numbers of randomizations whether we can reject the
hypothesis that the true p-value is below the approximate BHq
cutoff upper-bound of $0.1\times 44/600$ (the 44 comes from 40
nonzeros with 10% false discoveries). For instance, after just 10
randomizations, if the approximate p-value so far is greater than
0.2, we can reject the hypothesis that the exact p-value is below
$0.1\times 44/600$ at significance 0.0001 (and thus compute no
further randomizations for that covariate). Speed-ups like this
are possible in practice, although they require having an idea of
how many rejections will be made, which is not generally available
ahead of time.

## Appendix E Robustness Simulations

In Figure [8](#S6.F8), the points labeled “$100\alpha\%$
Emp. Cov” represents knockoffs run using the following convex
combination of the true and empirical covariance matrices:

$$ $\bm{\Sigma}_{\text{EC}}=(1-\alpha)\bm{\Sigma}+\alpha\hat{\bm{\Sigma}},$ $$

where $\hat{\bm{\Sigma}}$ is the empirical covariance
matrix. The point labeled “Graph. Lasso” represents knockoffs run
using the following covariance estimate:

$$ $\bm{\Sigma}_{\text{GL}}=\operatorname{diag}(r)\,\hat{\bm{\Theta}}^{-1}\operatorname{diag}(r)$ $$

where the rescaling vector $r$ has
$r_{j}=\sqrt{\Sigma_{jj}/(\hat{\bm{\Theta}}^{-1})_{jj}}$, and
$\hat{\bm{\Theta}}$ is the inverse covariance estimated by the
graphical Lasso with penalty parameter chosen by 2-fold
cross-validation on the log-likelihood.

## Appendix F Genetic analysis of Crohn’s disease

The SNP arrays came from an Affymetrix 500K chip, with calls made by
the BRLMM algorithm Affymetrix, (). SNPs satisfying any
of the following conditions were removed:

- •
minor allele frequency (MAF) $<1\%$,
- •
Hardy–Weinberg equilibrium test p-value $<0.01\%$,
- •
missing $>5\%$ of values (values were considered missing if
the BRLMM score was $>0.5$),
- •
position was listed as same as another SNP (this occured just
for the pair rs16969329 and rs4886982; the former had smaller MAF
and was removed),
- •
position was not in genetic map.

Furthermore, subjects missing $>5\%$ of values were removed. Following
the original paper describing/using this data (WTCCC,, ), we
did not adjust for population structure. Any missing values that
remained after the above preprocessing were replaced by the mean of
the nonmissing values for their respective SNPs.

## Appendix G Knockoffs with selection

First, recall that the results of Theorem [3.4](#S3.Thmtheorem4) hold if for
any subset $S\subset\mathcal{H}_{0}$, we have

$$ $([\bm{X},\,\tilde{\bm{X}}]_{\textnormal{swap}(S)},y)\,\,{\buildrel d\over{=}}\,\,([\bm{X},\,\tilde{\bm{X}}],y).$ (G.1) $$

In fact, MX knockoffs are defined in such a way that this property
holds. Now, the procedure employed in Section [7](#S7) to
construct knockoffs is slightly different from that described in the
rest of the paper. Explicitly, the data looks like

$$ $\bm{X}=\begin{bmatrix}\bm{X}^{(1)}\\ \bm{X}^{(2)}\end{bmatrix},\quad y=\begin{bmatrix}y^{(1)}\\ y^{(2)}\end{bmatrix},$ $$

where $y^{(1)}$ is $983\times 1$, $\bm{X}^{(1)}$ is $983\times 71,145$,
$y^{(2)}$ is $3930\times 1$ and $\bm{X}^{(2)}$ is $3930\times 71,145$;
recall that the set of samples labeled (1) is used to select the
cluster representatives, and that the two sets (1) and (2) are
independent of each other. The knockoffs $\tilde{\bm{X}}^{(2)}$ for
$\bm{X}^{(2)}$ are generated as described in the paper (using the ASDP
construction) and for (1), we set $\tilde{\bm{X}}^{(1)}=\bm{X}^{(1)}$. To
verify ([G.1](#A7.E1)), it suffices to show that for any subset
$S\in\mathcal{H}_{0}$ and each $i\in\{1,2\}$,

$$ $([\bm{X}^{(i)},\,\tilde{\bm{X}}^{(i)}],y^{(i)})\,\,{\buildrel d\over{=}}\,\,([\bm{X}^{(i)},\,\tilde{\bm{X}}^{(i)}]_{\text{swap}(S)},y^{(i)})$ (G.2) $$

since the sets (1) and (2) are independent. ([G.2](#A7.E2)) holds
for $i=2$ because we are following the classical knockoffs
construction. For $i=1$, ([G.2](#A7.E2)) actually holds with
exact equality, and not just equality in distribution.

We can also argue from the perspective of MX knockoffs from Definition
[3.1](#S3.Thmtheorem1). By construction, it is clear that
$\tilde{\bm{X}}^{(2)}$ are valid MX knockoffs for $\bm{X}^{(2)}$. We thus
study $\tilde{\bm{X}}^{(1)}$: since $\tilde{\bm{X}}^{(1)}={\bm{X}}^{(1)}$,
the exchangeability property is trivial; also,
$\tilde{\bm{X}}^{(1)}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}y^{(1)}\,|\,\bm{X}^{(1)}$ since
$\bm{X}^{(1)}$ determines $\tilde{\bm{X}}^{(1)}$.

## Peer Reviews

confidence: 4
contribution: 1
initial_rating: 1
presentation: 2
questions: - can authors discuss the work's limitations and potential impact on future work?
- can authors discuss potential overlap with prior work notably [1] (see above)?
- can authors adjust figure 1 with axis labels and explain why the number of samples (if I understand correctly) varies between epoch 1 and 5?
soundness: 1
strengths: - clear and well-written: the manuscript is well-written and easy to follow
- relevant topic: the authors tackle an interesting finding (i.e., training with random labels leads to substantial performance increase) that is relevant to the community and connected to multiple popular topics like self-supervised representation learning as well as the emerging idea of a universal representation 

[1] Bojanowski, Piotr, and Armand Joulin. "Unsupervised learning by predicting noise." International Conference on Machine Learning. PMLR, 2017.
[2] Reizinger, Patrik, et al. "Cross-Entropy Is All You Need To Invert the Data Generating Process." arXiv preprint arXiv:2410.21869 (2024).
[3] Huh, Minyoung, et al. "The platonic representation hypothesis." arXiv preprint arXiv:2405.07987 (2024).
summary: This paper explores the finding that training neural networks with random labels leads to substantial performance improvements in comparison to randomly initialised neural networks. Authors claim that these experimental findings (supported by empirical evidence & overlapping with prior work -- see below) highlight that the process of learning from data occurs independently of human-imposed structure and inform a novel perspective on the way neural networks work and do not discuss the work's limitations. Authors go on to proposing the use of random labels to fine-tuned pre-trained backbones to improve downstream generalisation.
weaknesses: - strong overlap with non-cited work/lack of novelty: the authors centered their work around the observation that random labels offers substantial performance improvement which they claim is a novel finding ("We completely removed the structure from the learning process by randomizing the class labels, and found that the model actually was able to learn from data despite the complete randomization and even performed better from a generalization perspective."). In fact, this observation has been presented and discussed in several works in the past, including [1], which is not cited by the authors.
- soundness of claims: the paper makes bold claims about "how neural networks learn" and what drives this process ("we present as provocative claim that the process of learning from data happens independently of human-imposed structures. To support this, we introduce the bold alternative hypothesis called the “Pan for Gold”. ").  These claims remain conjectures and hypothesis which are only supported by empirical evidence that the network learns from random labels which does not prove the author's "pan for gold" hypothesis. Additionally, authors further justify the relevance of their work by relying on GradCam visualisation, a method proven to be unreliable -- as also mentioned by authors.
- confidence in empirical findings: while the paper is well-written and clear, there is a lack of polishing of figures and of empirical results which impedes clarity and well as confidence in empirical results (e.g., missing axis labels, randomly masked out portions of curves, single seed experiments, core findings in section one are conducted on two small scale datasets and a single architecture type).
- missing sections: the authors omit important sections to their work including a related work section and a discussion of the paper's limitations.


[1] Bojanowski, Piotr, and Armand Joulin. "Unsupervised learning by predicting noise." International Conference on Machine Learning. PMLR, 2017.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. **Theoretical Justification:** Can the authors provide a rigorous theoretical framework to support the "Pan for Gold" hypothesis? Specifically, how does SGD with unstructured labels in overparameterized models lead to meaningful generalization, and what are the underlying mechanisms?
2. **Experimental Validation on Larger Datasets:** Have the authors considered testing the PUL algorithm on larger and more diverse datasets to validate the generality of their claims? Small-scale datasets may not capture the complexities of modern deep-learning tasks.
3. **Comparative Analysis with Baselines:** How does the PUL algorithm perform compared to existing state-of-the-art methods in unsupervised domain adaptation and object discovery?
4. **Clarity on Pan for Gold Hypothesis:** The hypothesis seems to conflate the effects of noise and regularization in SGD with meaningful learning from unstructured data. Can the authors clarify how their hypothesis differs from existing theories on overparameterization and implicit regularization?
soundness: 1
strengths: * **Challenging Conventional Wisdom:** The paper attempts to question the traditional beliefs about the necessity of structured labels and loss functions in deep learning, which is an interesting and bold endeavor.

* **Novel Hypothesis Introduction:** The "Pan for Gold" hypothesis is a creative metaphor that could inspire new ways of thinking about generalization in deep learning.

* **Exploration of Unstructured Labels:** Investigating the effects of training with unstructured labels is an intresting approach that could uncover overlooked aspects of model training dynamics.
summary: This paper introduces the "Pan for Gold" hypothesis, which challenges the traditional view that structured labels and well-defined loss functions are essential for deep learning models to learn meaningful representations and generalize well. The authors propose that generalization emerges naturally through the stochasticity inherent in SGD when training overparameterized models, even with unstructured (randomized) labels. They suggest that SGD acts like panning for gold, where valuable features are naturally sifted out from noise without relying on human-imposed structures.
To support this hypothesis, the authors conduct experiments where models are trained on datasets with completely randomized labels. Surprisingly, these models still learn meaningful features, as evidenced by improved performance over random initialization. They analyze this phenomenon using the NTK framework and observe a "swing phenomenon," where model outputs fluctuate significantly during early training stages.
Based on these observations, they introduce the PUL algorithm. They demonstrate PUL's effectiveness in tasks like unsupervised domain adaptation and object discovery. Additionally, they suggest that PUL mitigates issues like massive activation in vision transformers, aiding in model quantization.
weaknesses: * **Lack of Theoretical Rigor:** The paper makes strong claims without providing a solid theoretical foundation. The mathematical analysis is superficial and does not rigorously justify the "Pan for Gold" hypothesis or explain why unstructured labels should lead to better generalization.
* **Insufficient Empirical Evidence:** The experimental evaluation is limited and inadequate to support the bold claims made. Experiments are conducted on small datasets like MNIST, CIFAR-10 and SVHN, which are not representative of modern large-scale tasks. The performance improvements reported are marginal and could be due to experimental noise.
* **No Comparison with Baselines:** The paper fails to compare the proposed PUL algorithm with established baselines or state-of-the-art methods in the respective tasks. Without such comparisons, it's impossible to assess the significance of the results or attribute improvements to the proposed method.
* **Overgeneralization of Findings:** The authors make sweeping generalizations about deep learning based on limited and specific experiments. The claim that generalization emerges naturally through SGD in overparameterized models trained with unstructured labels is not convincingly demonstrated.
Methodological Issues: Key details about the experimental setup are missing or unclear, hindering reproducibility. For example, the process of assigning unstructured labels, hyperparameter settings, and specifics of the PUL algorithm are not adequately described.

Weak Analysis of Results: The paper lacks a thorough analysis of the results. It does not explore alternative explanations for the observed phenomena or consider confounding factors. The interpretations often rely on anecdotal observations rather than rigorous investigation.
Ambiguous Writing and Clarity Issues: The paper is difficult to follow in several sections due to ambiguous explanations and poor organization. Key concepts are not clearly defined, and the narrative lacks coherence, making it challenging to understand the proposed ideas fully.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: 1. Are the unstructured labels fixed during training or regenerated each epoch? Have you also experimented with changing labels at each epoch? This would be an important experiment as it could lead to completely different learning dynamics since the network cannot memorize stable image-label pairs.

2. In Sections 4.1 and 4.2, how does the performance change with longer training periods? The paper only shows results with 2-3 epochs, but longer training analysis is necessary to understand the stability and effectiveness of the method.

3. How was the optimal number of random classes determined in the PUL algorithm? Did you perform any experiments with different numbers of classes?
soundness: 2
strengths: 1. The paper presents a novel and interesting perspective on deep learning generalization by proposing a new hypothesis about the role of stochasticity in learning meaningful features
2. The proposed methodology is remarkably simple yet demonstrates effectiveness, requiring only random labels and a few additional training steps
3. The theoretical analysis through Neural Tangent Kernel provides mathematical insights into the learning dynamics and supports the main hypothesis
4. The experimental results show significant performance improvements across various applications including domain adaptation and object discovery tasks, demonstrating the practical utility of the proposed method
summary: This paper presents a new hypothesis about generalization in deep learning, suggesting it's not about learning structured patterns in data but rather a "Pan for Gold" process where SGD naturally filters useful features, and proposes the PUL algorithm utilizing random labels, demonstrating performance improvements in domain adaptation and object discovery tasks.
weaknesses: 1. The paper critically lacks essential experimental details needed for reproduction, including the specific method of generating unstructured labels, exact model architecture, hyperparameter settings, and detailed training procedures, making it difficult to validate the claims independently.

2. The paper lacks clear explanation about whether unstructured labels remain fixed during training. Based on the paper's content, it appears that labels are fixed, in which case neural networks would inevitably learn visual features in the process of memorizing image-label pairs, as they need to recognize some visual patterns to distinguish between images even with random labels. This suggests that learning meaningful features might be a natural consequence of the memorization process rather than the proposed "Pan for Gold" hypothesis.

3. The theoretical analysis is insufficient as the paper lacks in-depth discussion on why the "Pan for Gold" process leads to good generalization, focusing merely on describing phenomena without explaining the underlying mechanisms

4. The experimental validation is limited, lacking analysis of performance with longer training epochs in Sections 4.1 and 4.2, missing ablation studies on the number of unstructured labels, and failing to provide sensitivity analysis for various hyperparameters.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. Could you provide the definition of structure in this paper?
2. How do you distinguish the impact of stochasticity with the impact of unsupervised/semi-supervised representation learning/model-inductive-bias/compression?
3. In Table one, "we applied transfer learning to the frozen encoder." Could you elaborate how is transfer learning carried out?

Minor:
In Fig 2, "red" should be "black"
soundness: 2
strengths: The paper is well-written and presents an interesting analysis of generalization through stochasticity. I like the idea of decoupling the impact of human-imposed labels in generalization.

The paper carries out principled analysis on neural tangent kernel that demonstrates the swing phonemenon in the learning process. The visualizations of learning process through gradCAM and analysis on saliency maps are interesting.
summary: This paper investigates the generalization of deep learning methods when faced with random labels. It argues that generalization is not about learning the structure of data (X, Y), but rather follows a stochastic process that initially fluctuates before converging to a stable function space, akin to "panning for gold." Experiments on unsupervised domain adaptation show that using random labels with KL regularization outperforms the source-only baseline, which does not apply any adaptation. Exploratory analysis also shows that the proposed random labels reduce outliers in the attention map, resulting in a more balanced distribution that is more suitable for quantization.
weaknesses: Although I liked the exploration of stochasticity and its role in learning, I think the paper places too much emphasis on supervised data and stochasticity while neglecting other important factors in unsupervised and semi-supervised learning. The experiments are relatively weak because they rely on a supervised model trained with ground-truth data and KL regularization, which does not fully demonstrate the impact of stochasticity.

**1. Definition of structure, the role of supervised labels, and considerations for unsupervised/semi-supervised learning**

The paper frequently mentions "structure" but does not provide a clear definition. It appears to me that "structure" refers to supervised class labels assigned by humans, and the paper argues that the "structure" itself is not the essence of the 'gold' result and that learning also happens with random labels.

I disagree with the notion that "the goal of deep learning is to learn from data according to structures defined by humans." I think the supervised signal is only one source of information that models use to learn. Other sources of information include data itself as used in unsupervised or self-supervised learning, information from a model decision boundary manifested through unlabeled data as in transductive and semi-supervised learning setups, and assumptions about the world such as convolution for image processing. The paper solely focuses on the structure from supervisory signals and ignores other sources of information, which also plays a crucial role in learning that could be attributed to stochasticity in this paper.

**2. Experiments**

The experiments focus on domain adaptation where a supervised model has been trained on ground-truth labels in the source domain, and the goal is to adapt the model to a closely related target domain. This setup weakens the empirical results because the model is initialized with a learned representation from ground-truth labels and does not fully demonstrate the impact of stochasticity in a learning-from-scratch setup. Moreover, the model only performs small adjustments due to the constraint of KL regularization that penalizes the model for deviating from the source model. It is well-known that, in semi-supervised learning, a model could outperform the source-only baseline without any target labels. It is unclear if the improves is from stochasticity or semi-supervised learning (clustering assumption, KL, transduction through batchnorm).

---

confidence: 4
contribution: 1
initial_rating: 1
presentation: 1
questions: see above
soundness: 1
strengths: The paper presents the observation that training on data with random labels as an unsupervised per-training scheme can learn useful “features”, though the details of how are left very vague.
I am not familiar enough with this area of literature to know how novel this claim is. To me it does not seem surprising that training on random labels leads to a better model initialisation than the no pertaining.
Unfortunately the written quality of the paper is very low and thus it is difficult to determine if the paper really offers any strengths.
summary: The paper presents the observation that training on data with random labels can learn useful “features”, and the author suggest doing this as an unsupervised per-training scheme, though the details of how are left very vague. Unfortunately the written quality of the paper is very low, it is difficult to determine much more than this.
weaknesses: I apologise if English is not the first language of the authors but the written quality of the paper is sadly unacceptable for ICLR. At the moment the paper is almost impossible to follow. While there are a lot of machine learning terms mention the sentences are so unspecific and vague, lacking in concrete definitions and jumping from place to place it makes assessing the ideas of the authors impossible for me. I do not know if LLM’s have been used in the writing of the paper, or in translating from another language but the final quality is just not good enough at the moment.

The issues with the write are as follows:

1. Use of vague terms, which are not defined. For example “structure”, “function speed”, “Noisy features”, “Gold features”

2. Vague handy wavey claims not back up with references: “This aligns with physical intuition, where relaxation of tension leads to the stabilization of the space, and high-energy regions, like artifacts, are naturally eliminated.”

3. Reference to a “Panning through Unstructured Label (PUL) algorithm” which is never really defined, you have to try and glean what the algorithm is from passing comments.

4. The details of the experiment are extremely vague. For example. The full text exampling an experiment is: “Table 5 presents the object discovery performance on various trained models on ResNet50, including ImageNet pretrained, DINO (Caron et al., 2021a), and ImageNet pretrained weights further trained using our method. As can be seen from the results, even with just a three epochs of training using unstructured labels, performance can be easily improved. ”

5. Baselines for the experiments are not explained, the experiments jump from place to place with very little intuition and justification of why important choices were made We would like to recommend the paper is rewritten with a greater leave of specificity, intuition and detail before being resubmitted.


Other.

The experiments are extremely limited, typically only considering a single model, data set with no mention of repeats or the training procedure
Baselines for the experiments are extremely limited a single baseline is used, with little detail of what this baseline was and why it should lead to a fair comparison
In the conclusion the authors claim they introduce a “bold alternative hypothesis called the “Pan for Gold”,” but after reading the paper I’m am left wondering what this current hypothesis that “Pan for Gold” is mean to be an alternative too? This is never explained in any level of detail.