Year: 2018

Abstract:
A new statistical procedure (Model-X Candès et al. (2018)) has provided a way to identify important factors using any supervised learning method controlling for FDR. This line of research has shown great potential to expand the horizon of machine learning methods beyond the task of prediction, to serve the broader needs in scientific researches for interpretable findings. However, the lack of a practical and flexible method to generate knockoffs remains the major obstacle for wide application of Model-X procedure. This paper fills in the gap by proposing a model-free knockoff generator which approximates the correlation structure between features through latent variable representation. We demonstrate our proposed method can achieve FDR control and better power than two existing methods in various simulated settings and a real data example for finding mutations associated with drug resistance in HIV-1 patients.

heading: 1 INTRODUCTION
text: In the past decades, the machine learning methods have achieved great advancement in improving prediction accuracy. However, prediction accuracy is not the sole quest of ’big data’ analysis. A universal aim across a lot of scientific disciplines is to identify the causes of certain outcome. For example, in new drug development, a personalized medicine strategy can be formed by identifying the genetic markers that is biologically linked to certain drug response or resistance. In electronic health record analysis, the goal is to identify ’actionable’ factors for the purpose of reducing cost and improving the quality of care. In epidemiology or sociology studies, it is of interest to identify protection or risk factors, especially those that can be changed through public policy or civil planning to improve public welfare. In these scientific endeavors, the need is to identify the important factors associated with a certain outcome, so that further confirmatory investigation (e.g. randomized international studies) could be conducted to expand knowledge or change future actions for a better outcome. For this purpose, we are interested in procedures that control the Type I error, which is the chance of a false discovery. With a large amount of factors to test, controlling the chance of any false discovery (family-wise error) is too stringent. Therefore the objective is usually relaxed to control the False Discovery Rate (FDR) (Benjamini & Hochberg, 1995).
A recent breakthrough in the statistical theory, i.e. the Model-X framework(Barber & Candès, 2015; Candès et al., 2018), provides a general solution. It can be incorporated with any machine learning models to select true signals associated with the outcome, with rigorous control for FDR. This line of research sheds light on expanding the horizon of supervised learning methods beyond prediction. However, the implementation of Model-X requires the generation of the so called ’knockoffs’, which has limited existing methods. The goal of our paper is to fill in the gap by proposing a model-free method for generating knockoffs with relaxed assumptions for the data, leveraging the power of the recent development in deep generative models. Specifically, we propose a framework to generate knockoffs based on latent variable Z which captures the correlation structure of X . This provides a great generalization from the existing knockoff generation methods for hidden Markov model (HMM) and mixture Gaussian to (but not limited to) Hidden Markov Random Field. The contribution of this paper is it propose a deep learning framework to generate Model-X knockoffs through latent variable. We provide theoretical justification of the approach and also demonstrate the state of art variational auto encoder (VAE) can achieve promising results in the task of FDR controlled variable selection. This paper is among the first works for the new application of representative learning
in generating model-X knockoffs. Our discussion and preliminary results for goodness-of-fit also shed lights on future improvement of VAE for generating knockoffs.

heading: 2 BACKGROUND
text: The model-X framework (Candès et al., 2018) is proposed for the following problem. There are N i.i.d. samples (X(i), Y (i)) from a population, where the predictors are random variables X = (X1, . . . , Xp). The outcome Y only depends on a subset of predictors S ∈ {1, 2, . . . , p}, i.e. given {Xj , j ∈ S}, Y is independent of {Xj , j /∈ S}. The smallest set S that satisfies this requirement is considered as the true predictors. The goal is to find a good estimator Ŝ for S with control of the FDR = E ( #{j:j∈Ŝ\S} #{j:j∈Ŝ}∨1 ) . Barber & Candès (2015) first proposed the method for linear regression with Gaussian errors, where the design matrix X is assumed to be fixed. Candès et al. (2018) greatly generalized it to the Model-X framework. The traditional statistical models focusing on the conditional distribution Y |X1, . . . , Xp. In contrast, Model-X makes no assumption for the conditional distribution, but assumeX has a known distribution for the purpose of generating knockoffs X̃ that satisfies Definition 2.1.
Definition 2.1 Model-X knockoffs for the family of random variables X = (X1, . . . , Xp) are a new family of random variables X̃ = (X̃1, . . . , X̃p) constructed with the following two properties: (1) for any subset S ⊂ {1, . . . , p}, (X, X̃)swap(S) =d (X, X̃); (2) X̃ ⊥⊥ Y |X if there is a response

heading: Y. (2) is guaranteed if X̃ is constructed without looking at Y.
text: The following proposition 1 provide an equivalent condition for checking Model-X knockoffs.
Proposition 1 (Candès 2018) The random variables X̃ = (X̃1, . . . , X̃p) are Model-X knockoffs for X = (X1, . . . , Xp) if and only if the joint distributions are equal after swapping Xj and X̃j , i.e. (Xj , X−j , X̃j , X̃−j) =d (X̃j , X−j , Xj , X̃−j), for any j ∈ {1, . . . , p}, where X−j denote the other variables. And Y |= X̃|X .
Suppose one can generate Model-X knockoffs, the features are selected through the model-specific feature statistics Wj = wj([X, X̃], y). The requirement for the feature statistic is that the distribution for the Wj’s from the null variables (Xj 6∈ S) need to be symmetric. This is formalized as the flip-sign property,
Wj([X, X̃]swap(S), y) = (1− 21{j∈S})wj([X, X̃], y).
According to Thm 3.4 in Candès et al. (2018), the following procedure (Knockoff) can estimate the Ŝ controlling for modified FDR, i.e. E
( #{j:j∈Ŝ\S}
#{j:j∈Ŝ}+q−1
) ≤ q. First, compute the threshold
τ = min { t > 0 :
#{j:Wj≤−t} #{j:Wj≥t} ≤ q
} , then the selected set is
(Knockoff): Ŝ = {j :Wj ≥ τ}. (1)
A slightly more conservatively procedure (Knockoff+) increase the negative number in numerator by 1, and is able to control for the usual FDR. The threshold for Knockoff+ is τ+ = min { t > 0 : 1+#{j:Wj≤−t} #{j:Wj≥t} ≤ q } , and the selected set is
(Knockoff+): Ŝ = {j :Wj ≥ τ+}. (2)
Notice that the Knockoff approach controls for the modified FDR, which might have more power than Knockoff+, but it does not guarantee the control for the usual FDR. In the experiments of this paper, we used the Knockoff+.
An example of the feature statistics is the signed max lambda statistics in L1 penalized regressions Barber & Candès (2015): for the concatenated design matrix [X, X̃], define β̂(λ) = argminb { Loss(y, [XX̃]b) + λ‖b‖L1 } . Intuitively, the true predictors with strong signals would
have entered the LASSO selection with large λ and Xj is entering the model with larger λ than its
knockoff X̃j . The null variables will enter the model with smaller penalty. Thus the signed max lambda statistic is defined as Wj = (Zj ∨ Z̃j)sign(Zj − Z̃j). Where Zj = sup{λ : β̂2j−1(λ) 6= 0} and Z̃j = sup{λ : β̂2j(λ) 6= 0}. Various feature statistics has been proposed for common learning methods (Candès et al., 2018; Gimenez et al., 2018), however the current knockoff generation methods are still limited for real data from unknown distributions. Although Candès et al. (2018) provides a general algorithm for generating Model-X knockoffs, i.e. the Sequential Conditional Independent Pairs algorithm, it needs to sample from the conditional distribution of X̃j from f(Xj |X−j , X̃1:j−1) which becomes intractable with slightly complex distributional assumption for X . Candès et al. (2018) also proposes a second-order approach by matching first two moments of (X, X̃)swap(S) and (X, X̃). The second-order approach generate the true knockoff when X is from Gaussian distribution which have no high-order moments. Sesia et al. (2018) proposes an algorithm to sample knockoffs when X is from Hidden Markov Model. Gimenez et al. (2018) proposes algorithm that applicable for X from simple Bayesian network models (such as mixture Gaussian) with distributional assumptions for certain conditional probabilities which need to be estimated and sampled from.

heading: 3 A MODEL-FREE KNOCKOFF GENERATOR WITH LATENT ENCODING VARIABLES.

heading: 3.1 EXISTENCE OF DESIRED LATENT VARIABLE Z
text: Conceptually, the Model-X aims at controlling for the FDR of testing the null hypotheses Xj |= Y |X−j , where j = 1, . . . , p, and X−j denotes the variables except for Xj . Thus, a necessary condition for Xj to be chosen by Model-X is that it contains signals not contained in the other variables, i.e. Xj can not be a determinant function of X−j . In the knockoff generating process, any component of Xj’s that is independent from X−j is easy to be generated. The difficulty and necessity to consider the complex conditional distributions come from the correlated component. Thus we formed a general framework based on the assumption that there exists some latent variable Z so that Xj’s are mutually independent conditional on Z.
Such Z always exists for any X , one can take Z = X or Z = f(X), such that f is an one-onone mapping between X and Z. However, a larger difference in X̃j and Xj leads to more powerful knockoffs. In our framework, to generate X̃ different from originalX , our framework requires there exists Z such that X|Z is mutually independent and X|Z is not degenerated, i.e. V ar(X|Z) > 0 at least for a subset of X with non-zero measure. In the following, we provide motivational examples followed by a lemma to prove such Z exists for general X .
A general class of models satisfies this assumption is the Hidden Markov Random Field as illustrated by the undirected graphical models in Figure 1. The latent variables are from a Markov Random Field and the observed random variables Xj’s are generated from Z ′s with independent errors. However, our method is not limited for Hidden Markov Random Field, Z could exist for the non
Gibbs Random Field (X is not from Markov Random Field). For example, f(x1, x2) = 21{x1 + x2 ≤ 1, x1 ≥ 0 and x2 ≥ 0}. Consider the partition of the support {x1 + x2 ≤ 1, x1 ≥ 0 and x2 ≥ 0} into countable many squares. The desired Z can be defined as the indicator which squarer X belong to. X|Z is then the uniform distribution on that square, X1 and X2 are independent given Z. The desired latent variable Z always exists for non-degenerate categorical X . If X is from binomial (for binary data) or multinomial distribution (for categorical data) with positive probability mass, the desired latent variable Z always exists as a categorical variable with probabilities from the tensor decomposition of multi-way contingency array of probability forX Bhattacharya & Dunson (2012),i.e., p(x) = ∑ k πk ∏ j p(xj |z = k). Then Xi|Z’s are generated from Bernoulli or multinomial with parameters as functions of Z and thus the conditional density p(Xi|Z)’s are not degenerate.
For continuous or infinite discrete state X , De Finetti’s theorem states that if the distribution (X1, X2, . . . , Xn) is exchangeable then there exists a one dimensional non-degenerate latent variable Z that satisfies our requirement. Thus the desired Z exists for the multivariate normal distribution or trans-elliptical distribution (i.e., any distribution that can be transformed into multivariate normal distribution by element-wise transformation). The mixture of trans-elliptical distribution also has the desired latent Z as the combination of the single latent variable which indicate the cluster and the latent Z for each cluster.
More generally, we can prove the desired Z exists with the following Lemma 3.1. In the lemma, we consider X from a non-degenerate distribution with density f(x) with respect to a measure µ for a sigma field defined on Rp. The following lemma showed that one can decompose f into mixture of densities fk, where at least one fk has non-zero variance.
Lemma 3.1 For any non-negative density function f : Rp → R whose support contains an interior point, there exist non-negative functions fk : Rp → R, k = 1, 2, · · · such that f(x) = ∑ k πkfk(x) with πk > 0 and fk satisfies the following: 1)fk’s can be factorized into functions of each coordinates of x, i.e., fk(x) = ∏ j fkj(xj); 2) at least one component fk has non-zero variance, i.e.
V (fk) = ∫ x2fk(x)µ(dx)− [ ∫ xfk(x)µ(dx)] 2 > 0.
Proof: Since f(x)’s discontinuous points are at most countable, there exists an interior point of the support x∗ at which f(x) is continuous. Therefore there exists a cubic area B(x∗, ) = {x : |xj − x∗j | ≤ , j = 1, · · · , p} such that for any x ∈ B(x∗, ), f(x) ≥ f(x∗)/2 > 0 and π = P (x ∈ B(x∗, )) > 0. So we have f(x) = f(x
∗) 2 ∏ j I(xj ∈ [x∗j − , x∗j + ]) + g(x) where g(x) is a non-
negative function. This can be viewed as a mixture of a uniform distribution on the cubic area with probability f(x
∗)π 2 and a distribution with density g(x)
1− f(x ∗)π 2
with probability 1− f(x ∗)π 2 . Notice that
the uniform component is mutually independent. The remaining non-negative g(x) can be further decomposed into ∑∞ k=2 πkfk, where fk’s are the sequence of step functions to approximate the g−1(x). So now we have f(x) = f(x ∗)
2 ∏ j I(xj ∈ [x∗j − , x∗j + ])+ ∑∞ k=2 πkfk and the first term
has positive variance, which finish the proof.
Remark: The desired Z thus exist as the latent indicator of which component fk, X is from. Where X|Z and Z have the following densities (X|Z = k) ∼ fk(x)∫
fk(x)µ(dx) and P (Z = k) = πk ∫ fk(x)µ(dx). At least for the uniform component, X|Z is not degenerated.
To sum up, the desired latent variable Z generally exists for non degenerate continuous and discrete X .

heading: 3.2 METHOD AND ALGORITHM
text: Our proposed procedure has two steps, the first step is to generate Z̃ from the conditional distribution Z|X and then generate X’s knockoff from X|Z = Z̃. Specifically, consider the non-degenerate random vector of covariates X = (X1, X2, . . . , Xp) are from pX(x). And let Z be the vector of latent random variables Z from distribution pZ(z). Each component Xj is from p j (fi(Z)), j = 1, . . . , p, where j’s are mutually independent given Z. For example, for modeling continuous data, Xj ∼ N (fj(Z), σ2j ); and for modeling binary data, Xj ∼ Bernoulli(fj(Z)). Given the observed
samples {x(i)}Ni=1 from pX(x), to estimate pZ(z) and pX|Z may not be an identifiable problem. However, we just need to find one Z that approximately satisfies the above condition in practice for our purpose of knockoff generation.
In practice, we assume Z and the conditional distributions are from parametric families and the parameters can be approximated through variational inference. We denote two working models from parametric families: a) an encoder QZ|X(z|x; θ) to approximate pZ|X(z|x); b) a decoder QX|Z(x|z; f, q ) to approximate pX|Z(x|z). Where θ denotes the parameter for the encoder, and f denotes the parameters for decoder. And q denotes the working noise distribution to generate X̃ given f̂(Z̃), and is element-wise independent q = ∏ j q j . And θ̂, f̂ denote the estimates for these
parameters. Then we propose to generate knockoff X̃ as Algorithm 1.
Algorithm 1 Auto-encoding Knockoff Generator.
1: Train model parameters to get estimates θ̂, f̂ . 2: Generate Z̃ from QZ|X(z|X; θ̂). 3: Generate X̃ from QX|Z(x|Z̃; f̂ , q ).
Notice that the Step 2 in Algorithm 2 is different from the usual VAE implementation. In the latter, Z̃ is usually generated from a prior distribution pZ(z), however in our algorithm Z̃ is generated from the conditional distribution QZ|X . The following Theorem 3.2 justifies why our procedure is so.
By fitting a working model from parametric families with X|Z mutually independent, the good reconstruction of X from this model is a surrogate objective for X̃ is the knockoff of X . Thus the state of art Variational Auto Encoders (VAEs) (Kingma & Welling, 2014; Rezende et al., 2014; Maddison et al., 2017; Jang et al., 2017) provide a way to estimate the working models in Step 1. It also worth noting that Gao et al. (2018) point out that under certain conditions, the ELBO for VAE is equivalent to a lower bound for the Total Correlation Explanation, which is optimized when Z fully decorrelate (disentangle) the X (i.e. conditional of Z, X is mutually independent).
We explain the connection between our notation and the VAE methods in the following two examples. In the common variational auto encoder, Z ∼ N(0, Iq), the working model QZi|X is assumed to be from N (µi(X), σ2i (X)), and QXi|Z is N (fi(Z), σ2), where µi(•) and σ2i (•) and f(•) are functions to be approximated by neural networks. The q in this case is i.i.d. normal distribution with infinitesimal variance. So in practice people generate X̃ as f̂(Z̃). Another example is to approximate discrete Z’s by the concrete distribution Jang et al. (2017); Maddison et al. (2017). In this case X̃ is generated from Bernoulli(f̂(Z̃)) and = X̃ − f̂(Z̃) are element-wise independent. Remark: The objective function of our learning framework targets directly on the reconstruction of X , the data’s compatibility with the working models for Z|X and X|Z. Alternatively, if one construct the knockoffs directly from deep neural network to approximate X̃ = f(X,E) where E represents the generated random noise. The objective is then to minimize the the discrepancy between two distributions (X, X̃) and (X, X̃)swap(S). Since S is any subset of {1,. . . , p}, there are 2p − 1 ways of swapping. Thus the objective should be a composite function that incorporate all the swapping. This approach may have advantage over our approach when the ’working’ variational auto encoding model is not compatible with the true data distribution. Our proposed framework may have advantage in less computational burden and easiness to implement. The pros and cons of our method and various other deep learning based knockoff generator await to be evaluated and compared case by case in various real data sets.

heading: 3.3 THEORETICAL RESULTS
text: Theorem 3.2 shows that if X is mutually independent conditional on Z, if given the true conditional distributions Z|X and X|Z, one can first generate Z̃ from Z|X and then generate X̃ from X|Z = Z̃. The X̃ is the model-X knockoff of X . Theorem 3.3 provides the FDR bound for estimated conditional distribution Z|X and X|Z.
Theorem 3.2 For any vector of random variables Z, such that conditional on Z, Xi’s are mutually independent. If Z̃ is from pZ|X(•|X), then X̃ from pX|Z(•|Z̃) is model-X knockoffs.
Proof: For any j ∈ 1, . . . , p and sets Aj , Bj , A−j , B−j , P (Xj ∈ Aj , X̃j ∈ Bj , X−j ∈ A−j , X̃−j ∈ B−j) = ∫ pX|Z(Aj , A−j |Z)pX|Z(Bj , B−j |Z)dF (Z)
= ∫ pXj |Z(Aj |Z)pX−j |Z(A−j |Z)pXj |Z(Bj |Z)pX−j |Z(B−j |Z)dF (Z)
= ∫ pXj |Z(Bj |Z)pX−j |Z(A−j |Z)pXj |Z(Aj |Z)pX−j |Z(B−j |Z)dF (Z)
= ∫ pX|Z(Bj , A−j |Z)pX|Z(Aj , B−j |Z)dF (Z) = P (Xj ∈ Bj , X̃j ∈ Aj , X−j ∈ A−j , X̃−j ∈ B−j)
We have shownXj and X̃j are exchangeable in the joint distribution, then X̃ is the Model-X knockoff by Proposition 1.
Now we consider the FDR bound for the knockoffs generated from the estimated encoder and decoder. We assume the working models and observed distribution ofX are approximately compatible as in the conditions of the following theorem. We then apply the results per Barber et al. (2018) to show FDR control, where they showed the FDR derived from approximate knockoffs is bounded by a function of the observed KL divergence between (Xj , X̃j , X−j , X̃−j) and (X̃j , Xj , X−j , X̃−j) as below:
K̂Lj = log
( P (Xj , X̃j , X−j , X̃−j)
P (X̃j , Xj , X−j , X̃−j)
) =
n∑ i=1 log
( P (Xij , X̃ij , Xi,−j , X̃i,−j)
P (X̃ij , Xij , Xi,−j , X̃i,−j)
)
Theorem 3.3 Assume the observed marginal distribution pX(x) and the working models (QZ|X , QX|Z , q estimated in Algorithm 1) are approximately compatible as following: there exists a random vector Z from certain distribution pZ(z) and an → 0, such that
sup x | log
( q (x)
p̂(x)
) | ≤ an and sup
z,x | log ( pZ|X̂(z|x) QZ|X(z|x) ) | ≤ an.
Where the density of ̂ = X − f̂(Z) is denoted as p̂, considering f̂ as a fixed function. And denote X̂ = f̂(Z) + as a random variable generated from the decoder QX|Z(x|z = Z; f̂ , q ).
Then the FDR can be controlled at q exp{8na2n + 8 √ n log(p)an}.
Proof of Theorem 3.3 is in Appendix A.
Remark: Theorem 3.3 provides FDR bound given the variational models is able to approximate the Z|X and X|Z sufficiently close. However, there is no guarantees for deep generative model to converge with the conditions as in Theorem 3.2. And generally, deep learning knockoff generators gain the flexibility of no distributional assumptions at the cost of no current theoretical guarantee directly from convergence of algorithm to control of FDR. We will discuss this further in section 7. And we provide two ways for evaluating knockoff generated from our proposed framework in the following section 4.

heading: 4 GOODNESS OF FIT
text: One evaluation approach proposed by Candès et al. (2018) is to estimate the empirical FDR for simulated Y from the X of the real data, and compare the power and FDR control performance of knockoffs from various generator. Another method is to compose some goodness of fit metrics based on the distributional discrepancy between (X, X̃) and (X, X̃)swap. We will discuss these metrics further in section 7. In addition to the metrics shared with other knockoff generator, our framework has a unique measurement of goodness of fit based on the mutual independence condition
for X given Z. Specifically, if the q is assumed to be additive error, i.e. X̃j = f̂j(Z̃) + j . The mutual independence of the residuals X − f̂(Z̃), which is the difference of the real data X and its fitted means, indicate good model fitting. In the following numerical studies, we focus on the first evaluation strategy which is based on FDR control for simulated Y . We include preliminary graphical check for the independence criterion in Appendix C and further discussion in Section 7.

heading: 5 SIMULATION
text: In the simulation, we compare three knockoff generation schemes: a) our proposed method implemented through VAE b) the fixed knockoff generation method in Barber & Candès (2015) c) the second order matching in Candès et al. (2018). We demonstrate the performance of these methods in two simulation scenarios. In both scenarios, the sample size is 200, and number of potential predictors is 100, where the first m variables are the true signals. The coefficients β for the true predictors are alternating ρ and −ρ, where ρ is the magnitude of the signal. In Figures 23, the power and FDR are drawn as curves with respect to ρ.
We generated two types of outcomes: the continuous outcomes are generated with errors from standard normal distribution; the binary outcomes are generated from Bernoulli with probabilities 1/(1 + exp(−Xβ)). We use the signed max lambda statistics in Barber & Candès (2015) (see definition on page 2). We implemented the L1 penalized linear and logistic regression with R package ’glmnet’ correspondingly. We used the Knockoff+ as in (2) to control for FDR. Results based on 100 replications are presented for two different numbers of true signals m = 10, 20.
Setting 1. In this setting, we simulate a set of lightly correlated non-normal distributed continuous predictors: a) Generate independent uniform (0, 1) distributed 200 by 200 random matrix. b) Let C be the Cholesky decomposition of the correlation matrix with 0.1 on the off diagonal entries and compute Z = UC. c) Even column X2i is generated as Z4i+0.5∗Z34i+1, and odd column X2i+1 is generated as Z4i+2− 0.5Z24i+2+0.5 exp(Z4i+3), where i = 0, . . . , 49. d) Rescale X’s to be within range [0, 1] by subtracting the column min and divided by column max−min. Empirically, the pairwise correlation X ranges from −0.3 to 0.6 and the absolute value of the correlation has mean 0.07. The signal to noise ratios of the Gaussian case approximately has the range of 0.25 − 0.3ρ2 for m = 10 and 0.7− 0.85ρ2 for setting with m = 20. Setting 2. The second setting generates correlated categorical predictors :a) Generate independent standard Gaussian distributed 200 by 100 random matrix. b) Let C be the Cholesky decomposition of the correlation matrix with 0.1 on the off diagonal entries. Let Z = UC. c) X is the categorized Z matrix, where Xij = 1(Zij > 0). Approximately, the signal to noise ratio for the Gaussian case is 2ρ2 for m = 10 and 4ρ2 for m = 20. We present results for FDR controlled at 0.1.
Results According to Figures 2 3, our proposed method implemented through VAE has FDR controlled below the predetermined threshold, and it demonstrates higher power than its two competitors. Cases with more true signals m = 20 demonstrate better power for Gaussian outcomes, since we control for the proportion of false discoveries, the criterion can tolerate more false discoveries with larger number of true discoveries it detect. The second-order matching method proposed for the gaussian predictors has very low power due to the violation of its assumption that our generated X are not from Gaussian distribution. The assumptions for Fixed Knockoff generator is satisfied for the linear regression settings(with gaussian errors) but violated for the logistic regression setting. In the simulation, the FDR is still controlled for cases these assumptions are violated. However the powers are low.
Implementation details. For both settings, we implemented our method with VAE Kingma & Welling (2014). We defer the training details to Appendix B.

heading: 6 REAL DATA EXPERIMENT
text: In the real data experiment, we demonstrate the performance of our proposed method on generating knockoffs for sparse genetic mutation data. The dataset and scientific problem is from a study aiming at detecting mutations associated with drug resistance in patients with Human Immunodeficiency Virus Type 1 (HIV-1) Rhee et al. (2006). Due to the space limit, here we present results
for the 7 protease inhibitors. There is no ground truth for which subset of mutations caused the drug resistance. Nevertheless, there is a set of 73 treatment-selected mutations (TSMs) identified from a separate study Rhee et al. (2005), which are selected as the mutations marginally correlated with the patient treatment history with protease inhibitors. Thus to evaluate the performance of the knockoffs, we first simulates the outcome for which we know the true signal, and evaluate the FDR control and power as shown in Figure 4, and then we used the real data outcome and compared our selected mutations with the TSM. Notice that the TSM is not drug specific and can not be considered as ground truth, by showing number of selected mutations in TSM set we are not showing the FDR control but demonstrate reproducibility across studies.
The normal distributed Z does not fit the discrete and sparse nature of the genetic mutation data. Thus we adopted the Categorical Variation Auto-Encoder (CAT-VAE) Jang et al. (2017); Maddison et al. (2017), where in our implementation, the latent variablesZ are 20 Gumbel-Softmax distributed variable with temperature 1, each with 10 categories; both the encoder and decoder has one hidden layer of dimension 200.
Figure 4 demonstrates the results for simulated outcomes using the real data X . Overall speaking, our proposed method still achieves the FDR control and demonstrates the highest power among the three. The correlation is smaller than the previous simulation settings (since the X is a sparse matrix), so the other two methods has better power. Controlled for a FDR level of 0.2, all three methods achieve a high power of 80% very large signal at ρ = 10 in the Gaussian case. For the binary case, our proposed method shows greater advantage in achieving more than 2-times the power of the other two methods.
Table 1 presents the number of selected and matched selection of mutations with the 73 TSMs. Here we present results controlled by both Knockoff+ (as in (2)) and Knockoff (as in (1)) controlled at FDR level 0.2. Similar to the simulated settings, our proposed method select more mutations
than the other methods, the average proportion of selected mutation that is in the TSM set for our proposed method is 75.3%, which demonstrate reproducibility across studies. With less conservative Knockoff procedure, Fixed and Second-Order methods are able to select more mutations, where our proposed method seems be less sensitive to choice of Knockoff or Knockoff+.
Implementation details For these analysis, the CAT-VAE knockoffs are generated for 186 mutations with (≥ 4) occurrence in the data set and the number of patients is 846. Since each drug resistance outcome is observed in a subset of samples, when conduct analysis for each drug, we only consider mutations with≥ 2 occurrence. Following analysis in Rhee et al. (2006), the outcome is a continuous variable, which is the log transformed drug susceptibility. For results in Figure 4, the linear signal is simulated as following: in each replication, randomly chose 50 predictors from 186 mutations , and compute µ as sum of ρxj’s with alternating signs. The Gaussian outcome is µ plus a standard normal noise and the binomial outcome is generated from Bernoulli( 11+exp(−µi) ).

heading: 7 CONCLUSION AND DISCUSSION
text: In this paper, we propose a knockoff generating framework based on a latent variable Z to decorrelate X . We noticed there are two independent works on deep learning generators for knockoffs while writing and revising the paper. (For double blind review consideration, we exclude the citations at this time.) Compared with the other deep learning frameworks for generating knockoffs, our framework is the only framework directly targeting on the reconstruction of X from Z. We provide theoretical results for general existence of such Z. Our framework offers a computational light solution, and if the data and the parametric families for the working models are compatible, our framework may offer a computationally efficient solution to be easily implemented by domain scientists.
Since deep learning approaches mainly target on real data from complex distributions, which does not have tractable conditional distributions. There is a lack of theoretical FDR guarantee for all the current deep knockoff generators. One common obstacle is the non-convexity of deep learning objectives. To solve this problem, an arXiv paper (name omitted for double blind review) proposed several goodness of fit statistics to compare deep generative models. It calls for future research to compare these deep knockoff generators in terms of their power and FDR control performance in simulated Y ’s and also compare them via various goodness-of-fit metrics.
In addition to these metrics, the structure of our framework provides an unique diagnostic metric, which is the mutual independence of X|Z as described in section 4. We include the some result for graphical check in Appendix C, which demonstrates the estimated Z does attenuateX’s correlation. In both simulation settings and the real HIV dataset, the attenuation for the binary X is closer to independence compared with the continuous simulatedX . It remains for further investigation that if we can adjust the objective function of VAE to improve the goodness-of-fit in terms of independence. For example, adding a penalization term with the Maximum Mean Discrepancy metric between distribution ofX−f(Z̃) and the element wise independent distribution from the product of marginal distributions.
To sum up, this paper proposes a framework of deep knockoff generator based on latent variables. This framework offers a computationally light solution to generate knockoffs. The promising numerical results and the easiness of implementation show make it a promising tool among the choices of deep generating models for FDR controlled variables selection.

heading: A. PROOF OF THEOREM 3.2
text: Proof: Consider ˜̂X which is another independent sample from QX|Z(x|z = Z; q ) when sampling X̂ using the same Z. Since q is element-wise independent, conditions in Theorem 2.1 holds, thus ˜̂ X is a model-X knockoff of X̂ and
log
( P (X̂j ∈ Aj , ˜̂Xj ∈ Bj , X̂−j ∈ A−j , ˜̂X−j ∈ B−j)
P (X̂j ∈ Bj , ˜̂Xj ∈ Aj , X̂−j ∈ A−j , ˜̂X−j ∈ B−j)
) = 0.
Notice that by our assumptions
sup x | log
( pX̂(x)
pX(x)
) | ≤ sup
x,z | log ( pX̂|Z(x|z) pX|Z(x|z) ) | = sup x,z | log ( q (x− f̂(z)) p̂(x− f̂(z)) ) | ≤ an
sup x,x̃ | log
( p ˜̂ X|X̃ (x̃|x)
pX̃|X(x̃|x)
) | = sup
x,z,x̃ | log ( ∫ q (x̃− f̂(z))pZ|X̂(z, x)dz∫ q (x̃− f̂(z))QZ|X(z, x)dz ) | ≤ sup z,x | log ( pZ|X̂(z|x) QZ|X(z|x) ) | ≤ an
So supx,x̃ | log ( p (X̂, ˜̂ X) (x,x̃)
p(X,X̃)(x,x̃)
) | ≤ 2an and we can bound the log likelihood ratio by
log ( P (Xj ∈ Aj , X̃j ∈ Bj , X−j ∈ A−j , X̃−j ∈ B−j) P (Xj ∈ Bj , X̃j ∈ Aj , X−j ∈ A−j , X̃−j ∈ B−j) )
= log
( P (X̂j ∈ Aj , ˜̂Xj ∈ Bj , X̂−j ∈ A−j , ˜̂X−j ∈ B−j)
P (X̂j ∈ Bj , ˜̂Xj ∈ Aj , X̂−j ∈ A−j , ˜̂X−j ∈ B−j)
) + log ( P (Xj ∈ Aj , X̃j ∈ Bj , X−j ∈ A−j , X̃−j ∈ B−j)
P (X̂j ∈ Aj , ˜̂Xj ∈ Bj , X̂−j ∈ A−j , ˜̂X−j ∈ B−j)
)
+ log ( P (X̂j ∈ Bj , ˜̂Xj ∈ Aj , X̂−j ∈ A−j , ˜̂X−j ∈ B−j)) P (Xj ∈ Bj , X̃j ∈ Aj , X−j ∈ A−j , X̃−j ∈ B−j) ) ≤ 0 + 2an + 2an = 4an
So we can bound
max j=1,··· ,p
K̂Lj ≤ 8na2n + 8 √ n log(p)an.
And FDR has the bound in the statement by applying Lemma 2 in Barber et al. (2018). Specifically, if an = o((n log(p))−1/2), then FDR will be asymptotically bounded by q + o(1).
B. IMPLEMENTATION DETAILS FOR NEURAL NETWORKS IN SIMULATION
The latent variables are multivariate normal with 300 dimensions. The models were trained with batch size 25 and for 20 epochs with the default ‘adam’ optimizer in ‘keras’Chollet et al. (2015). The architecture of VAE for setting 1 is as following, the encoder networks for µz and log(σ) have two hidden layers with 500 and 400 neurons with ‘tanh’ activation and L2 regularization with tuning parameter 0.2. The activation for the output layer is ‘linear’. The decoder network has two hidden layers with the same specifications, followed by a batch-normalization layer and output with linear activation. The architecture of VAE for setting 2 is as following, the encoder networks has two hidden layers with 500 and 400 neurons with ‘relu’ activation and L2 regularization with tuning parameter 0.3. The activation for the output layer is ‘linear’. The decoder is an one layer network with the ‘sigmoid’ activation function and then threshold by 0.5.
C. ILLUSTRATION FOR DISENTANGLEMENT OF Z .
Here we present the Q-Q plots for the p-values for pairwise correlation test for X and X|Z. Where the correlation ofX|Z is estimated as the pairwise correlation ofX−f̂(Z̃). WhenX|Z are mutually independent, the p-values should be from the uniform distribution and the Q-Q plot will coincide with the diagonal line y = x. For the simulated settings, the figure is the result from one replication.
We see for the HIV real data example in Figure A1 (e) (f), the correlation is greatly attenuated by conditioning on Z. For the simulation studies, the Q-Q plot in setting 2 (Figure A1 (c) and (d)) almost coincide with the diagonal line after conditioning on Z, which indicates good fit of the model. Interestingly, in the above two settings, X’s are binary. In simulation setting 1 (Figure A1 (a) and (b)), when X’s are continuous, the disentanglement effect is still visible but not as strong as setting 2 and the HIV example.
(a) Q-Q plot for p-values of correlation of X in simulation setting 1
(b) Q-Q plot for p-values of correlation of X − f̂(Z̃) in simulation setting 1
(c) Q-Q plot for p-values of correlation of X in simulation setting 2
(d) Q-Q plot for p-values of correlation of X − f̂(Z̃) in simulation setting 2
(e) Q-Q plot for p-values of correlation of X in the HIV example (f) Q-Q plot for p-values of correlation of X − f̂(Z̃) in the HIV example
Figure A1. Illustration of Attenuation Effect of Z for X’s Correlation.

Review summaries:
This manuscript tackles an important problem, namely, generating the knockoff procedure for FDR-controlled feature selection so that it can work on any data, rather than only for data generated from a Gaussian (as in the original work) or several specific other cases (mixtures of Gaussians or HMMs).  The basic idea is to exploit ideas from variational autoencoders to create a generic knockoff generation mechanism. Specifically, the main idea of the paper is to map the input covariates X into latent variable Z using a variational auto-encoder, generate the knockoffs \tilde{Z} in the latent space, and then map \tilde{Z} back to the input space to get the knockoffs \tilde{X}. The authors claim that their contribution in threefold: