Aggregating Knockoffs for False Discovery Rate Control with an   Application to Gut Microbiome Data

Fang Xie; Johannes Lederer

arXiv:1907.03807·stat.ME·March 2, 2021·Entropy

Aggregating Knockoffs for False Discovery Rate Control with an Application to Gut Microbiome Data

Fang Xie, Johannes Lederer

PDF

1 Repo

TL;DR

This paper introduces an aggregated knockoff method to improve false discovery rate control in high-dimensional microbiome data analysis, leading to new biological insights about obesity-related gut bacteria.

Contribution

The paper proposes a novel aggregation approach for knockoffs that enhances power while maintaining FDR control in microbiome studies.

Findings

01

Aggregated knockoffs increase discovery power.

02

Several overlooked phyla are linked to obesity.

03

Method validated through theory and simulations.

Abstract

Recent discoveries suggest that our gut microbiome plays an important role in our health and wellbeing. However, the gut microbiome data are intricate; for example, the microbial diversity in the gut makes the data high-dimensional. While there are dedicated high-dimensional methods, such as the lasso estimator, they always come with the risk of false discoveries. Knockoffs are a recent approach to control the number of false discoveries. In this paper, we show that knockoffs can be aggregated to increase power while retaining sharp control over the false discoveries. We support our method both in theory and simulations, and we show that it can lead to new discoveries on microbiome data from the American Gut Project. In particular, our results indicate that several phyla that have been overlooked so far are associated with obesity.

Figures40

Click any figure to enlarge with its caption.

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Equations42

\mathrm{FDR}:=\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q}[\boldsymbol{X},\boldsymbol{y}]\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q}[\boldsymbol{X},\boldsymbol{y}]|\vee 1}\Biggr{]}

\mathrm{FDR}:=\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q}[\boldsymbol{X},\boldsymbol{y}]\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q}[\boldsymbol{X},\boldsymbol{y}]|\vee 1}\Biggr{]}

\mathrm{power}:=\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q}[\boldsymbol{X},\boldsymbol{y}]\cap\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q}[\boldsymbol{X},\boldsymbol{y}]|\vee 1}\Biggr{]}\,,

\mathrm{power}:=\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q}[\boldsymbol{X},\boldsymbol{y}]\cap\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q}[\boldsymbol{X},\boldsymbol{y}]|\vee 1}\Biggr{]}\,,

FDR \leq q,

FDR \leq q,

x_{i} ∣ x_{i} \sim N (μ_{i}, V) for each i \in {1, \dots, n}

x_{i} ∣ x_{i} \sim N (μ_{i}, V) for each i \in {1, \dots, n}

μ_{i}

μ_{i}

V

\widehat{\boldsymbol{\beta}}[\tau,\boldsymbol{X},\widetilde{\boldsymbol{X}},\boldsymbol{y}]\in\arg\min\limits_{\boldsymbol{\alpha}\in\mathbb{R}^{2p}}\biggl{\{}\mskip 2.0mul\Bigl{[}\boldsymbol{\alpha}|[\boldsymbol{X}~{}\widetilde{\boldsymbol{X}}],\boldsymbol{y}\Bigr{]}+\sum_{j=1}^{2p}h_{\tau}[\alpha_{j}]\mskip 2.0mu\biggr{\}}\,,

\widehat{\boldsymbol{\beta}}[\tau,\boldsymbol{X},\widetilde{\boldsymbol{X}},\boldsymbol{y}]\in\arg\min\limits_{\boldsymbol{\alpha}\in\mathbb{R}^{2p}}\biggl{\{}\mskip 2.0mul\Bigl{[}\boldsymbol{\alpha}|[\boldsymbol{X}~{}\widetilde{\boldsymbol{X}}],\boldsymbol{y}\Bigr{]}+\sum_{j=1}^{2p}h_{\tau}[\alpha_{j}]\mskip 2.0mu\biggr{\}}\,,

Z_{j} [X, y]

Z_{j} [X, y]

Z_{j} [X, y]

W_{j} := max {Z_{j}, Z_{j}} \cdot sign (Z_{j} - Z_{j}) for j \in {1, \dots, p} .

W_{j} := max {Z_{j}, Z_{j}} \cdot sign (Z_{j} - Z_{j}) for j \in {1, \dots, p} .

T_{q}:=\min\biggl{\{}\mskip 2.0mut\in\mathcal{W}:\frac{\#\{j:W_{j}\leq-t\}}{\#\{j:W_{j}\geq t\}\vee 1}\leq q\mskip 2.0mu\biggr{\}}\,,\\

T_{q}:=\min\biggl{\{}\mskip 2.0mut\in\mathcal{W}:\frac{\#\{j:W_{j}\leq-t\}}{\#\{j:W_{j}\geq t\}\vee 1}\leq q\mskip 2.0mu\biggr{\}}\,,\\

T_{q}^{+}:=\min\biggl{\{}\mskip 2.0mut\in\mathcal{W}:\frac{1+\#\{j:W_{j}\leq-t\}}{\#\{j:W_{j}\geq t\}\vee 1}\leq q\mskip 2.0mu\biggr{\}}\,.

T_{q}^{+}:=\min\biggl{\{}\mskip 2.0mut\in\mathcal{W}:\frac{1+\#\{j:W_{j}\leq-t\}}{\#\{j:W_{j}\geq t\}\vee 1}\leq q\mskip 2.0mu\biggr{\}}\,.

S_{q}

S_{q}

S_{q}^{+}

S_{q, AKO} [k] := \cup_{i = 1}^{k} S_{q_{i}} (or S_{q, AKO}^{+} [k] := \cup_{i = 1}^{k} S_{q_{i}}^{+}) .

S_{q, AKO} [k] := \cup_{i = 1}^{k} S_{q_{i}} (or S_{q, AKO}^{+} [k] := \cup_{i = 1}^{k} S_{q_{i}}^{+}) .

\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q,\operatorname{AKO}}^{+}[k]\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q,\operatorname{AKO}}^{+}[k]|\vee 1}\Biggr{]}\leq q\,.

\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q,\operatorname{AKO}}^{+}[k]\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q,\operatorname{AKO}}^{+}[k]|\vee 1}\Biggr{]}\leq q\,.

\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q_{i}}^{+}\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q_{i}}^{+}|\vee 1}\Biggr{]}\leq q_{i}\,.

\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q_{i}}^{+}\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q_{i}}^{+}|\vee 1}\Biggr{]}\leq q_{i}\,.

\displaystyle\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q,\operatorname{AKO}}^{+}[k]\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q,\operatorname{AKO}}^{+}[k]|\vee 1}\Biggr{]}

\displaystyle\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q,\operatorname{AKO}}^{+}[k]\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q,\operatorname{AKO}}^{+}[k]|\vee 1}\Biggr{]}

\displaystyle=\mathbb{E}\Biggl{[}\frac{|\cup_{i=1}^{k}(\widehat{\mathcal{S}}_{q_{i}}^{+}\backslash\mathcal{S}^{*})|}{|\cup_{i=1}^{k}\widehat{\mathcal{S}}_{q_{i}}^{+}|\vee 1}\Biggr{]}

\displaystyle\leq\mathbb{E}\Biggl{[}\sum_{i=1}^{k}\frac{|\widehat{\mathcal{S}}_{q_{i}}^{+}\backslash\mathcal{S}^{*}|}{|\cup_{i=1}^{k}\widehat{\mathcal{S}}_{q_{i}}^{+}|\vee 1}\Biggr{]}

\displaystyle\leq\mathbb{E}\Biggl{[}\sum_{i=1}^{k}\frac{|\widehat{\mathcal{S}}_{q_{i}}^{+}\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q_{i}}^{+}|\vee 1}\Biggr{]}

\displaystyle=\sum_{i=1}^{k}\mathbb{E}\Biggl{[}\frac{|\widehat{\mathcal{S}}_{q_{i}}^{+}\backslash\mathcal{S}^{*}|}{|\widehat{\mathcal{S}}_{q_{i}}^{+}|\vee 1}\Biggr{]}

\leq i = 1 \sum k q_{i} = q,

y = X β + u .

y = X β + u .

Pr (y_{i} = 1) = \frac{exp ( x _{i}^{⊤} β )}{1 + exp ( x _{i}^{⊤} β )} for i \in {1, \dots, n} .

Pr (y_{i} = 1) = \frac{exp ( x _{i}^{⊤} β )}{1 + exp ( x _{i}^{⊤} β )} for i \in {1, \dots, n} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Code & Models

Repositories

LedererLab/aggregated_knockoffs
none

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Abstract

Recent discoveries suggest that our gut microbiome plays an important role in our health and wellbeing. However, the gut microbiome data are intricate; for example, the microbial diversity in the gut makes the data high-dimensional. While there are dedicated high-dimensional methods, such as the lasso estimator, they always come with the risk of false discoveries. Knockoffs are a recent approach to control the number of false discoveries. In this paper, we show that knockoffs can be aggregated to increase power while retaining sharp control over the false discoveries. We support our method both in theory and simulations, and we show that it can lead to new discoveries on microbiome data from the American Gut Project. In particular, our results indicate that several phyla that have been overlooked so far are associated with obesity.

keywords:

false discovery rate control; knockoffs; variable selection; gut microbiome

\pubvolume

23 \issuenum2 \articlenumber230

\externaleditorAcademic Editor: Jürgen Pilz \datereceived20 January 2021 \dateaccepted11 February 2021 \datepublished16 February 2021 \hreflinkhttps://doi.org/10.3390/e23020230

\TitleAggregating Knockoffs for False Discovery Rate Control with an Application to Gut Microbiome Data

\TitleCitationAggregating Knockoffs for False Discovery Rate Control with an Application to Gut Microbiome Data \AuthorFang Xie *\orcidA and Johannes Lederer \AuthorNamesFang Xie and Johannes Lederer \AuthorCitationXie, F.; Lederer, J.

\corres Correspondence: [email protected]

1 Introduction

Research on associations between the microbiome in the human gut and health and disease has surged in recent years (Evans et al., 2013; Huttenhower et al., 2012; Koliada et al., 2017; Ley et al., 2006). Data on the microbiome are abundant in view of citizen science endeavors such as the American Gut Project (Knight Lab, ), but these non-standard ways of data collection limit data quality.

Another difficulty in the analysis of microbiome data is the high dimensionality, which means that the number of parameters is large. The high-dimensionality is due to the diversity of the microbiome: at the phylum level, there are typically several dozen types of microbes; at the genus level, there are even several hundred types of microbes. Finding the microbes that are connected to a trait, therefore, requires the use of variable selection techniques. Consequently, in view of the low data quality and the high dimensionality, research on microbiome data is in particular need for controlling false discovery rates.

As a specific example where false discovery control can be important, we are interested in finding phyla and genera of microbes that are related to obesity. Successful detection of such groups of microbes could eventually lead to new means for weight control. We model the task as a variable selection problem in logistic regression with the obesity as the dependent variable and the (log-transformed) counts of the microbial relative abundancies as the variables. The number of parameters, that is, the number of potential phyla/genera, is large, while the number of actually relevant phyla/genera is assumed to be small; hence, we speak of sparse, high-dimensional logistic regression. A number of corresponding variable selection methods have been established, including especially $\ell_{1}$ -penalized logistic regression (Ng, 2004), which can be equipped with knockoffs to do FDR control (Barber and Candès, 2015). However, in our application, these standard pipelines perform insufficiently: for example, phyla that are known to be associated with obesity are missed, and few phyla are selected in the first place.

To increase the variable selection performance, we propose a simple aggregation scheme. It consists of two steps: the knockoff method is run $k$ times with a decreasing FDR target level, and the selections are then combined. We show that this aggregation scheme keeps FDR control intact while improving power in practice.

Our three key contributions are:

We introduce an aggregation scheme that provably retains the original methods’ guarantees—see Theorem 2.2. 2. 2.

We show numerically that the aggregation can increase the original methods’ power—see Sections 3.1 and 3.2. 3. 3.

We show that the resulting pipelines for FDR control can be readily applied to empirical data and lead to new discoveries—see Section 3.3.

The remainder of the paper is organized as follows. In Section 2, we introduce our pipeline and establish its theory. In Sections 3.1 and 3.2, we verify the accuracy of the pipeline in simulations. In Section 3.3, we then show the usefulness of the pipeline in selecting phyla and genera connected with obesity. In Section 4, we finally conclude with a brief discussion.

2 Methods and Theory

In this section, we introduce and study our aggregation scheme for knockoffs. We presume $n$ independent data pairs $(\boldsymbol{x}_{1},y_{1}),\dots,(\boldsymbol{x}_{n},y_{n})$ , where each $\boldsymbol{x}_{i}\in\mathbb{R}^{p}$ is a vector of variables and $y_{i}\in\mathbb{R}$ an outcome. We keep the relationship between the outcomes and variables completely general at this point— the relationship could be linear, logistic, or anything else—but we assume that this relationship is captured by a parameter $\boldsymbol{\beta}\in\mathbb{R}^{p}$ . Our target for inference is then the active set $\mathcal{S}^{*}:=\{j:\beta_{j}\neq 0\}$ .

Two important quality measures for an estimate $\widehat{\mathcal{S}}_{q}[\boldsymbol{X},\boldsymbol{y}]$ of $\mathcal{S}^{*}$ are the FDR level

[TABLE]

and the power

[TABLE]

where $|\cdot|$ denotes the cardinality of a set. Our focus is on estimators that provide FDR control at level $q\in[0,1]$ , that is,

[TABLE]

while having large power. In the following, we recall that the knockoff filter provides FDR control. We then equip the knockoff filter with an aggregation step to improve its power.

2.1 A Brief Introduction to the Knockoff Filter

The main ingredient of the knockoff method (Barber and Candès, 2015) is a "knockoff version" $\widetilde{\boldsymbol{X}}\in\mathbb{R}^{n\times p}$ of the design matrix $\boldsymbol{X}$ . This new matrix $\widetilde{\boldsymbol{X}}$ is essentially $\boldsymbol{X}$ with additional noise such that (1) the estimator can distinguish predictors from $\boldsymbol{X}$ and $\widetilde{\boldsymbol{X}}$ but (2) both design matrices still have a similar correlation structure. The idea is then to compare the selections of predictors in $\widetilde{\boldsymbol{X}}$ and in $\boldsymbol{X}$ when estimating on the extended design $[\widetilde{\boldsymbol{X}}\,\boldsymbol{X}]$ .

Denoting $\widetilde{\boldsymbol{X}}=(\widetilde{\boldsymbol{x}}_{1},\ldots,\widetilde{\boldsymbol{x}}_{n})^{\top}$ and $\boldsymbol{X}=(\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})^{\top}$ , the underlying assumption is that $\boldsymbol{x}_{i}\sim\mathcal{N}(\boldsymbol{0}_{p},\boldsymbol{\Sigma})$ for a positive definite matrix $\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$ . The knockoffs are then generated according to (Barber et al., 2018)

[TABLE]

with

[TABLE]

where $\boldsymbol{a}\in\mathbb{R}^{p}$ is such that $\boldsymbol{V}$ is positive definite.

The variable selection is then based on the estimator

[TABLE]

where $l:\mathbb{R}^{p}\mapsto\mathbb{R}$ is a loss function, $[\boldsymbol{X}~{}\widetilde{\boldsymbol{X}}]\in\mathbb{R}^{n\times 2p}$ the extended design matrix, and $h_{\tau}:\mathbb{R}\mapsto[0,\infty)$ a penalty function with tuning parameter $\tau>0$ . Examples include the lasso (where $l\bigl{[}\boldsymbol{\alpha}|[\boldsymbol{X}~{}\widetilde{\boldsymbol{X}}],\boldsymbol{y}\bigr{]}:=|\!|{\boldsymbol{y}-[\boldsymbol{X}~{}\widetilde{\boldsymbol{X}}]\boldsymbol{\alpha}}|\!|_{2}^{2}$ and $h_{\tau}[\alpha_{j}]:=\tau|\alpha_{j}|$ ) and the penalized logistic regression (where $l\bigl{[}\boldsymbol{\alpha}|[\boldsymbol{X}~{}\widetilde{\boldsymbol{X}}],\boldsymbol{y}\bigr{]}=-\sum_{i=1}^{n}(y_{i}\boldsymbol{x}_{i}^{\top}\boldsymbol{\alpha}-\ln(1+\exp(\boldsymbol{x}_{i}^{\top}\boldsymbol{\alpha})))$ and $h_{\tau}[\alpha_{j}]:=\tau|\alpha_{j}|$ ).

The basic test statistics of the knockoff approach then capture the maximal tuning parameter of each variable entering the model:

[TABLE]

for $j\in\{1,\ldots,p\}$ . We surpress the dependence on $[\boldsymbol{X},\boldsymbol{y}]$ in the following for notational ease. The final test statistic $W:=(W_{1},\ldots,W_{p})^{\top}$ then combines the basic test statistics into

[TABLE]

This statistic compares how much earlier the original predictor enters the model as compared to the fake predictor—or the other way around. The threshold of the standard knockoff procedure for a given FDR level $q\in[0,1]$ is then

[TABLE]

where $\mathcal{W}:=\{|W_{j}|:j\in\{1,\ldots,p\}\}$ .

The knockoff procedure has another variant, called knockoff+, which differs in the threshold:

[TABLE]

These definitions finally yield the active sets

[TABLE]

The active sets $\widehat{\mathcal{S}}_{q}^{+}$ indeed provide FDR control at level $q$ , that is, satisfy inequality (1)—see Barber and Candès (2015), [Theorem 2]; the active sets $\widehat{\mathcal{S}}_{q}$ provide an approximate version of it—see Barber and Candès (2015), [Theorem 1].

2.2 Aggregating Knockoffs

We now introduce the aggregation scheme and its theory. The aggregation scheme applies the knockoff method $k$ times and combines the results:

Step 1: Given a target FDR $q\in[0,1]$ , choose a sequence $q_{1},\ldots,q_{k}\in[0,1]$ such that $q=\sum_{i=1}^{k}q_{i}$ . Apply the standard knockoff (or knockoff+) procedure $k$ times, once for each target FDR $q_{i}$ , and denote the corresponding $k$ -estimated active sets by $\widehat{\mathcal{S}}_{q_{1}},\ldots,\widehat{\mathcal{S}}_{q_{k}}$ (or $\widehat{\mathcal{S}}_{q_{1}}^{+},\ldots,\widehat{\mathcal{S}}_{q_{k}}^{+}$ ).

Step 2: Combine these $k$ -estimated active sets by taking the union:

[TABLE]

The intuition behind this scheme is that increasing the number of knockoffs stabilizes the outcome and improves the power. While applied here to the knockoff method, we emphasize that the aggregation scheme can be applied to every model and estimator as long as the there is a method for FDR control to start with. Hence, rather than the standard knockoffs as used below, we could equally well use model-X knockoffs (Candès et al., 2018), deep knockoffs (Romano et al., 2019), or KnockoffGAN (Jordon et al., 2019).

On the other hand, the aggregation scheme retains the knockoffs’ theoretical guarantees: {Theorem} Given a target FDR level $q\in[0,1]$ , the set $\widehat{\mathcal{S}}_{q,\operatorname{AKO}}^{+}[k]$ of the aggregation scheme for knockoff+ provides FDR control at level $q$ :

[TABLE]

Similarly, the scheme retains the approximate FDR control of standard knockoffs—we will demonstrate this in our simulations.

Proof of Theorem 2.2.

The proof is—maybe suprisingly—simple. By Theorem 2 in (Barber and Candès, 2015), we have for all $q_{i}$ , $i\in\{1,\ldots,k\}$ ,

[TABLE]

Hence,

[TABLE]

as desired. ∎

For $k=1$ , our method equals the standard knockoff (or knockoff+) procedure; in practice, we recommend $k\approx 5$ as a trade-off between computational effort and statistical effect. We also recommend $q_{i}=q/2^{i-1}$ , which, strictly speaking, does not meet our assumptions on $q_{1},\ldots,q_{k}$ exactly, but it works well empirically—see the simulations.

2.3 Other Approaches

While preparing this manuscript, two other ways of aggregating knockoffs were proposed (Holden and Hellton, 2018; Gimenez and Zou, 2019). While we use $k$ knockoffs in $k$ processes individually and aggregate at the end, they use multiple knockoffs simultaneously in one process. They can also show that their schemes provide valid FDR control, but we can argue that our approach is considerably simpler. We illustrate in the Appendix A.1 that we typically get more power than the multiple knockoffs method proposed by (Holden and Hellton, 2018); we have not gotten the scheme of (Gimenez and Zou, 2019) to run yet.

3 Simulations and a Real Data Analysis

In this section, we test our method empirically both on synthetic and on real data. We incorporate the specifics of microbiome data: First, in line with recent proposals for generalized linear models in this context (Lu et al., 2019), we study linear regression as well as logistic regression, and we transform count data with the standard log-transformation (Aitchison, 1982). Second, since gut microbiome data tends to be correlated (Naqvi et al., 2010), we ensure that the synthetic data are also correlated. Third, since gut microbime data also tend to be zero-inflated, we replace zero values by 0.5 times the observed minimum abundance, which is standard approach in microbiome analysis (Aitchison, 2003; Kurtz et al., 2015). We compare our modified aggregating knockoff pipeline applied to $\ell_{1}$ -penalized regression (called AKO henceforth) with the standard knockoff pipeline (KO).

Throughout, we set $k=5$ and $q_{i}=q/2^{i-1}$ . We also show the results for other choice of $q_{1},\dots,q_{k}$ in the Appendix A.2.

3.1 Simulation 1: Linear Regression

We first consider linear regression. The dimensions of the synthetic data are $(n,p)\in\{(200,100),(400,200)\}$ . The rows $\boldsymbol{x}_{i}$ of the design matrix are sampled i.i.d. from $\mathcal{N}(\boldsymbol{0}_{p},\boldsymbol{\Sigma})$ with the elements in $\boldsymbol{\Sigma}$ satisfying $\Sigma_{ij}=\rho^{|i-j|}$ with correlation factor $\rho=0.5$ . The noise is drawn from $\boldsymbol{u}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\boldsymbol{I}_{n})$ . The true parameter $\boldsymbol{\beta}$ has $20$ nonzero coefficients, which are selected uniformly at random from $\{1,\ldots,p\}$ and set to $1$ before the entire vector is rescaled to have signal-to-noise ratio ${|\!|{\boldsymbol{X}\boldsymbol{\beta}}|\!|_{2}^{2}}/n\sigma^{2}=5$ with $\sigma^{2}=1$ . The outcome $\boldsymbol{y}$ is then generated by the linear model

[TABLE]

The test statistic $W$ is based on the lasso method as described in the preceeding section.

We compare the empirical FDR and power averaged over $r=100$ repetitions of the simulations for each method. The ideal result would be an average FDR of at most $q$ and a power equal to 1.

The results in Figure 1a,b show that our pipeline retains the KO’s FDR control while increasing the power.

{paracol}

2 \switchcolumn

3.2 Simulation 2: Logistic Regression

We now consider logistic regression (the paper (Klose and Lederer, 2020) was the first one to apply our scheme to logistic regression). The above settings remain the same except for the outcomes $y_{i}$ being generated by

[TABLE]

The results in Figure 2a,b show again that our pipeline retains the KO’s FDR control while increasing the power. For more simulations in various settings, please refer to Appendix A.3.

{paracol}

2 \switchcolumn

3.3 Influence of the Gut Microbiome on Obesity

A well-functioning gut microbiome is essential for health (Huttenhower et al., 2012); for example, there is strong evidence that the composition of the microbiome is connected to obesity (Escobar et al., 2014; Gérard, 2016; Turnbaugh and Gordon, 2009). Existing research about this connection has focused only on phyla of bacteria that are abundant in most guts, such as Actinobacteria, Bacteroidetes, Firmicutes, and Verrucomicrobia (Bai et al., 2019; Clarke et al., 2012; Depommier et al., 2019; Gao et al., 2018; Koliada et al., 2017). In the following analysis, in contrast, we include the microbiome in its entirety. Our findings suggest that also phyla beyond those ones mentioned in the literature are connected to obesity.

The data for our analysis are from the American Gut Project (Knight Lab, ). The scope is American adults with age between 20 and 69 and BMI between $15\,\mathrm{kg/m^{2}}$ and $60\,\mathrm{kg/m^{2}}$ in the collection up to January, 2018. This includes $n=8404$ subjects, of which $278$ are underweight (uw; BMI below $18.5\,\mathrm{kg/m^{2}}$ ), $4972$ have normal weight (nor; 18.5–25 $\,\mathrm{kg/m^{2}}$ ), $2253$ are overweight (ow; 25–30 $\,\mathrm{kg/m^{2}}$ ), and $901$ are obese (ob; above $30\,\mathrm{kg/m^{2}}$ ). We transform the data and deal with the zero-counts as decribed earlier. The total number of phyla in our scope is $p=55$ .

The underlying model for these data is $\ell_{1}$ -penalized logistic regression model with outcomes $y_{i}=1$ (ob) when BMI $\geq 30$ and $y_{i}=0$ (non-ob) otherwise. The target FDR level is $q=0.1$ . Seven different groupings are considered to get the most out of the data: (i) all four groups (uw+nor+ow+ob), (ii) uw+ob, (iii) nor+ob, (iv) ow+ob, (v) uw+nor+ob, (vi) uw+ow+ob and (vii) nor+ow+ob. The AKO is applied with $\{q_{i}=q/2^{i-1}:i\in\{1,\ldots,k\}\}$ .

The results in Table 3.3 show that our pipeline selects more phyla in general. Since the theory and the above similations suggest that both methods have similar FDR, the results indicate that our pipeline has more power. In particilar, phyla that are known to be connected with obesity, such as Actinobacteria, Bacteroidetes, Firmicutes, and Verrucomicrobia (Bai et al., 2019; Clarke et al., 2012; Depommier et al., 2019; Gao et al., 2018; Koliada et al., 2017), are selected by AKO more often across the seven groupings.

The results also highlight Spirochaetes, which have not been associated with obesity in the literature, yet. The standard pipeline does not seem to have the power to select it, and similarly, two other known FDR control methods for microbiome data—the standard BH procedure (Benjamini and Hochberg, 1995) and the TreeFDR (Xiao et al., 2017)—select Spirochaetes only for large FDR levels ( $q\gtrsim 0.8$ ). (Futher comparisons with the BH procedure, TreeFDR, and the plain KO can be found in the Appendix A.4.) In contrast, our method selects Spirochaetes even at very small FDR levels (such as $q=0.01$ ), which strongly suggests a connection between Spirochaetes and obesity.

Our pipeline can, of course, also be applied to more detailed taxonomic ranks. As an illustration, we report the results of an application to genera for ALL—cf. (i) in Table 3.3—in Table 3.3. The data sampling is same with the phyla analysis. The only difference is that the total number of genera is $p=969$ . The rest of results for other six different groupings are given in the Appendix B. We find correspondingly that AKO selects more genera than the standard KO.

{specialtable}

[H] \tablesize \widetableSelected bacterial phyla by our pipeline (AKO) and the original pipeline (KO) at FDR level $q=0.1$ for seven groupings. AKO consistently selects more phyla than KO.

(i) all (ii) uw + ob

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslash

\PreserveBackslashBacteroidetes

\PreserveBackslash

\PreserveBackslashCyanobacteria

\PreserveBackslash

\PreserveBackslashCyanobacteria

\PreserveBackslash

\PreserveBackslashFirmicutes

\PreserveBackslashProteobacteria

\PreserveBackslash

\PreserveBackslashSpirochaetes

\PreserveBackslash

\PreserveBackslashSynergistetes

\PreserveBackslash

\PreserveBackslashSynergistetes

\PreserveBackslashTenericutes

\PreserveBackslash

\PreserveBackslashVerrucomicrobia

\PreserveBackslash

(iii) nor + ob (iv) ow + ob

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslash

\PreserveBackslashActinobacteria

\PreserveBackslashBacteroidetes

\PreserveBackslash

\PreserveBackslashCyanobacteria

\PreserveBackslash

\PreserveBackslashFirmicutes

\PreserveBackslash

\PreserveBackslashLentisphaerae

\PreserveBackslash

\PreserveBackslashProteobacteria

\PreserveBackslash

\PreserveBackslashProteobacteria

\PreserveBackslash

\PreserveBackslashSpirochaetes

\PreserveBackslash

\PreserveBackslashSpirochaetes

\PreserveBackslashSynergistetes

\PreserveBackslash

\PreserveBackslashSynergistetes

\PreserveBackslash

\PreserveBackslashTM7

\PreserveBackslashTenericutes

\PreserveBackslash

\PreserveBackslashVerrucomicrobia

\PreserveBackslash

\PreserveBackslashThermi

\PreserveBackslash

(v) uw + nor + ob (vi) uw + ow + ob

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslash

\PreserveBackslashActinobacteria

\PreserveBackslashBacteroidetes

\PreserveBackslash

\PreserveBackslashCyanobacteria

\PreserveBackslash

\PreserveBackslashFirmicutes

\PreserveBackslash

\PreserveBackslashLentisphaerae

\PreserveBackslash

\PreserveBackslashProteobacteria

\PreserveBackslash

\PreserveBackslashProteobacteria

\PreserveBackslash

\PreserveBackslashSpirochaetes

\PreserveBackslash

\PreserveBackslashSpirochaetes

\PreserveBackslashSynergistetes

\PreserveBackslash

\PreserveBackslashSynergistetes

\PreserveBackslash

\PreserveBackslashTM7

\PreserveBackslashTenericutes

(vii) nor+ow+ob

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslash

\PreserveBackslashActinobacteria

\PreserveBackslash

\PreserveBackslashBacteroidetes

\PreserveBackslash

\PreserveBackslashCyanobacteria

\PreserveBackslash

\PreserveBackslashProteobacteria

\PreserveBackslash

\PreserveBackslashSpirochaetes

\PreserveBackslash

\PreserveBackslashSynergistetes

\PreserveBackslash

\PreserveBackslashTenericutes

\PreserveBackslash

\PreserveBackslashVerrucomicrobia

\PreserveBackslash

{paracol}2 \switchcolumn

{specialtable}

[H] \tablesize Selected bacterial genera by our pipeline (AKO) and the original pipeline (KO) at FDR level $q=0.1$ for ALL—cf. (i) in Table 3.3. AKO selects more genera than the original KO.

\PreserveBackslashPhylum

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslashActinomyces

\PreserveBackslash

\PreserveBackslashCollinsella

\PreserveBackslash

\PreserveBackslashEggerthella

\PreserveBackslashCyanobacteria

\PreserveBackslashYS2

\PreserveBackslash

\PreserveBackslashStreptophyta

\PreserveBackslashFirmicutes

\PreserveBackslashBacillus

\PreserveBackslash

\PreserveBackslashLactobacillus

\PreserveBackslash

\PreserveBackslashLactococcus

\PreserveBackslash

\PreserveBackslashClostridium

\PreserveBackslash

\PreserveBackslashLachnospira

\PreserveBackslash

\PreserveBackslashRuminococcus

\PreserveBackslash

\PreserveBackslashPeptostreptococcaceae

\PreserveBackslash

\PreserveBackslashAcidaminococcus

\PreserveBackslash

\PreserveBackslashMegasphaera

\PreserveBackslash

\PreserveBackslashMogibacteriaceae

\PreserveBackslash

\PreserveBackslashErysipelotrichaceae

\PreserveBackslash

\PreserveBackslashCatenibacterium

\PreserveBackslashProteobacteria

\PreserveBackslashRF32

\PreserveBackslash

\PreserveBackslashHaemophilus

\PreserveBackslashTenericutes

\PreserveBackslash

\PreserveBackslashRF39

4 Discussion

Our aggregation scheme for knockoffs is supported by theory (Section 2) and simulations (Sections 3.1 and 3.2) and may lead to new discoveries in microbiomics (Section 3.3).

While we focus on a specific pipeline, our concept applies very generally. For example, it does not depend on the underlying statistical model or estimator but only on the availability of FDR control. In particular, the FDR control can be established via standard knockoffs or any other scheme. This flexibility is particularly interesting in practice: while the standard knockoffs rely on normality, other knockoff procedures, such as model-X knockoffs (Candès et al., 2018), deep knockoffs (Romano et al., 2019), and KnockoffGAN (Jordon et al., 2019), allow for much more general designs. Hence, the standard knockoffs can be readily swapped for those procedures in our pipeline (without any changes to the methodology or theory) when indicated in an application.

Considerably later than our paper has been put online, two other papers on the topic have also been put online (Srinivasan et al., 2020; Nguyen et al., 2020)—apparently without being aware of our manuscript. A comparison to those results would also be of interest.

\authorcontributions

Conceptualization, F.X. and J.L.; methodology, F.X. and J.L.; software, F.X.; validation, F.X. and J.L.; writing—original draft preparation, F.X.; writing—review and editing, J.L.. All authors have read and agreed to the published version of the manuscript.

\funding

This research received no external funding.

\institutionalreview

Not applicable.

\informedconsent

Not applicable.

\dataavailability

All the software, simulations, and data analyses are provided at github.com/lederlab. The data for our analysis are downloaded from the website of the American Gut Project, which is public.

Acknowledgements.

We thank Aaditya Ramdas and Sophie-Charlotte Klose for their valuable input. \conflictsofinterestThe authors declare no conflict of interest.\appendixtitlesyes \appendixstart

Appendix A Additional Explanations

A.1 Further Simulations for Comparison to Multiple Knockoffs (MKO) (Holden and Hellton, 2018)

We compare the three methods—KO, AKO, and MKO—in the logistic cases. The results shown in Figure 3 indicate that the MKO is more conservative than both our pipeline and the KO, and our pipeline always has higher power than KO and MKO.

A.2 Choice of $q_{1},\dots,q_{k}$

The theory applies to every choice of $q_{1},\dots,q_{k}$ that satisfies the condition in step 1 on page 2.2, but indeed, each choice has slightly different characteristics in practice. To give an example, we have added the result for $q_{i}=q/k$ (AKO.ave, Figure 4) to Figure 1a in the main text. We observe in general that the results do look different for each choice of the sequence but that the differences are usually only small.

A.3 Various Settings for the Simulation Part

Our conclusions hold over a wide spectrum of settings, including different dimensionalities, sparsity levels, correlations, and so forth. To illustrate this, we vary the settings of the linear case (see Figure 1a,b) once more in Figure 5 below.

A.4 Better than other Competitors (under the AGP Data)

Knockoff methods seem to be better suited for the microbiome data than their competitors. For illustration, we have included BH (Benjamini and Hochberg, 1995) and TreeFDR (Xiao et al., 2017) in the application that corresponds to Table 3.3 (i) in the paper. We find (in Table A.4) on those data more generally that the knockoff methods select many more predictors than other methods.

{specialtable}

[H] Selected bacterial phyla by four methods—BH, TreeFDR, KO, and AKO (correponds to Table 3.3 (i)).

(i) ALL

\PreserveBackslashBH

\PreserveBackslashTreeFDR

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslash

\PreserveBackslashActinobacteria

\PreserveBackslash

\PreserveBackslashBacteroidetes

\PreserveBackslash

\PreserveBackslashCyanobacteria

\PreserveBackslash

\PreserveBackslashProteobacteria

\PreserveBackslash

\PreserveBackslashSpirochaetes

\PreserveBackslash

\PreserveBackslashSynergistetes

\PreserveBackslash

\PreserveBackslashTenericutes

\PreserveBackslashVerrucomicrobia

\PreserveBackslash

\PreserveBackslashVerrucomicrobia

Appendix B Additional Results on the Genera Rank

We present in this section the results on the genera rank of different groupings; see Tables B–B.

{specialtable}

[H] Analysis at the genus level rank for the grouping (ii) uw+ob. AKO selects more genera than the original KO.

\PreserveBackslashPhylum

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslashCollinsella

\PreserveBackslashFirmicutes

\PreserveBackslash

\PreserveBackslashLachnospira

\PreserveBackslash

\PreserveBackslashAcidaminococcus

\PreserveBackslash

\PreserveBackslashCatenibacterium

\PreserveBackslashTenericutes

\PreserveBackslashRF39

{specialtable}

[H] Analysis at the genus level for the grouping (iii) nor + ob. AKO selects more genera than the original KO.

\PreserveBackslashPhylum

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslashActinomyces

\PreserveBackslash

\PreserveBackslashCollinsella

\PreserveBackslashCyanobacteria

\PreserveBackslashYS2

\PreserveBackslashFirmicutes

\PreserveBackslashBacillus

\PreserveBackslash

\PreserveBackslashLactococcus

\PreserveBackslash

\PreserveBackslashLachnospira

\PreserveBackslash

\PreserveBackslashRuminococcus

\PreserveBackslash

\PreserveBackslashAcidaminococcus

\PreserveBackslash

\PreserveBackslashMegasphaera

\PreserveBackslash

\PreserveBackslashMogibacteriaceae

\PreserveBackslash

\PreserveBackslashErysipelotrichaceae

\PreserveBackslash

\PreserveBackslashCatenibacterium

\PreserveBackslashProteobacteria

\PreserveBackslashRF32

\PreserveBackslash

\PreserveBackslashHaemophilus

\PreserveBackslashTenericutes

\PreserveBackslashRF39

\PreserveBackslash

\PreserveBackslashML615J-28

{specialtable}

[H] Analysis at the genus level for the grouping (iv) ow + ob. AKO selects more genera than the original KO.

\PreserveBackslashPhylum

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslashEggerthella

\PreserveBackslashCyanobacteria

\PreserveBackslashYS2

\PreserveBackslash

\PreserveBackslashStreptophyta

\PreserveBackslashFirmicutes

\PreserveBackslash

\PreserveBackslashBacillus

\PreserveBackslash

\PreserveBackslashClostridium

\PreserveBackslash

\PreserveBackslashLachnospira

\PreserveBackslash

\PreserveBackslashAcidaminococcus

\PreserveBackslash

\PreserveBackslash1-68

\PreserveBackslash

\PreserveBackslashErysipelotrichaceae

\PreserveBackslash

\PreserveBackslashCatenibacterium

\PreserveBackslashProteobacteria

\PreserveBackslashHaemophilus

{specialtable}

[H] Analysis at the genus level for the grouping (v) uw + nor + ob. AKO selects more genera than the original KO.

\PreserveBackslashPhylum

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslashActinomyces

\PreserveBackslash

\PreserveBackslashCollinsella

\PreserveBackslashCyanobacteria

\PreserveBackslashYS2

\PreserveBackslashFirmicutes

\PreserveBackslashBacillus

\PreserveBackslash

\PreserveBackslashLactococcus

\PreserveBackslash

\PreserveBackslashLachnospira

\PreserveBackslash

\PreserveBackslashRuminococcus

\PreserveBackslash

\PreserveBackslashAcidaminococcus

\PreserveBackslash

\PreserveBackslashMegasphaera

\PreserveBackslash

\PreserveBackslashMogibacteriaceae

\PreserveBackslash

\PreserveBackslashSHA-98

\PreserveBackslash

\PreserveBackslashErysipelotrichaceae

\PreserveBackslash

\PreserveBackslashCatenibacterium

\PreserveBackslashProteobacteria

\PreserveBackslashRF32

\PreserveBackslash

\PreserveBackslashHaemophilus

\PreserveBackslashTenericutes

\PreserveBackslashRF39

\PreserveBackslash

\PreserveBackslashML615J-28

{specialtable}

[H] Analysis at the genus level for the grouping (vi) uw + ow + ob. AKO selects more genera than the original KO.

\PreserveBackslashPhylum

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslashEggerthella

\PreserveBackslashCyanobacteria

\PreserveBackslashYS2

\PreserveBackslash

\PreserveBackslashStreptophyta

\PreserveBackslashFirmicutes

\PreserveBackslashBacillus

\PreserveBackslash

\PreserveBackslashLactobacillus

\PreserveBackslash

\PreserveBackslashClostridium

\PreserveBackslash

\PreserveBackslashLachnospira

\PreserveBackslash

\PreserveBackslashVeillonellaceaes

\PreserveBackslash

\PreserveBackslashAcidaminococcus

\PreserveBackslash

\PreserveBackslash1-68

\PreserveBackslash

\PreserveBackslashErysipelotrichaceae

\PreserveBackslash

\PreserveBackslashCatenibacterium

\PreserveBackslash

\PreserveBackslashEubacterium

\PreserveBackslashProteobacteria

\PreserveBackslash

\PreserveBackslashRF32

\PreserveBackslash

\PreserveBackslashHaemophilus

{specialtable}

[H] Analysis at the genus level for the grouping (vii) nor+ow+ob. AKO selects more genera than the original KO.

\PreserveBackslashPhylum

\PreserveBackslashKO

\PreserveBackslashAKO

\PreserveBackslashActinobacteria

\PreserveBackslashActinomyces

\PreserveBackslash

\PreserveBackslashCollinsella

\PreserveBackslash

\PreserveBackslashEggerthella

\PreserveBackslashCyanobacteria

\PreserveBackslashYS2

\PreserveBackslashFirmicutes

\PreserveBackslashBacillus

\PreserveBackslash

\PreserveBackslashLachnospira

\PreserveBackslash

\PreserveBackslashRuminococcus

\PreserveBackslash

\PreserveBackslashAcidaminococcus

\PreserveBackslash

\PreserveBackslashMegasphaera

\PreserveBackslash

\PreserveBackslashErysipelotrichaceae

\PreserveBackslash

\PreserveBackslashCatenibacterium

\PreserveBackslashProteobacteria

\PreserveBackslashRF32

\PreserveBackslash

\PreserveBackslashHaemophilus

\PreserveBackslashTenericutes

\PreserveBackslash

\PreserveBackslashRF39

\reftitle

References

Bibliography30

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Evans et al. (2013) Evans, J.M.; Morris, L.S.; Marchesi, J.R. The gut microbiome: the role of a virtual organ in the endocrinology of the host. J. Endocrinol. 2013 , 218 , R 37–R 47.
2Huttenhower et al. (2012) Huttenhower, C.; Gevers, D.; Knight, R.; Abubucker, S.; Badger, J.H.; Chinwalla, A.T.; Creasy, H.H.; Earl, A.M.; Fitz Gerald, M.G.; Fulton, R.S.; others. The Human Microbiome Project Consortium: Structure, function and diversity of the healthy human microbiome. Nature 2012 , 486 , 207–214.
3Koliada et al. (2017) Koliada, A.; Syzenko, G.; Moseiko, V.; Budovska, L.; Puchkov, K.; Perederiy, V.; Gavalko, Y.; Dorofeyev, A.; Romanenko, M.; Tkach, S. Association between body mass index and Firmicutes/Bacteroidetes ratio in an adult Ukrainian population. BMC Microbiol. 2017 , 17 , 120.
4Ley et al. (2006) Ley, R.E.; Turnbaugh, P.J.; Klein, S.; Gordon, J.I. Microbial ecology: human gut microbes associated with obesity. Nature 2006 , 444 , 1022.
5(5) Knight Lab. American Gut Project. Available online: http://americangut.org (accessed on 11 June 2019).
6Ng (2004) Ng, A.Y. Feature selection, L 1 vs. L 2 regularization, and rotational invariance. In Proceedings of the 21st International Conference on Machine Learning, Banff, Alberta, Canada 2004; p. 78.
7Barber and Candès (2015) Barber, R.F.; Candès, E.J. Controlling the false discovery rate via knockoffs. Ann. Stat. 2015 , 43 , 2055–2085.
8Barber et al. (2018) Barber, R.F.; Candès, E.J.; Samworth, R.J. Robust inference with knockoffs. ar Xiv 2018 , ar Xiv:1801.03896.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Code & Models

Videos

Abstract

keywords:

1 Introduction

2 Methods and Theory

2.1 A Brief Introduction to the Knockoff Filter

2.2 Aggregating Knockoffs

Proof of Theorem 2.2.

2.3 Other Approaches

3 Simulations and a Real Data Analysis

3.1 Simulation 1: Linear Regression

3.2 Simulation 2: Logistic Regression

3.3 Influence of the Gut Microbiome on Obesity

4 Discussion

Acknowledgements.

Appendix A Additional Explanations

A.1 Further Simulations for Comparison to Multiple Knockoffs (MKO) (Holden and Hellton, 2018)

A.2 Choice of q1,…,qkq_{1},\dots,q_{k}q1​,…,qk​

A.3 Various Settings for the Simulation Part

A.4 Better than other Competitors (under the AGP Data)

Appendix B Additional Results on the Genera Rank

A.2 Choice of $q_{1},\dots,q_{k}$