ECKO: Ensemble of Clustered Knockoffs for multivariate inference on fMRI data

TL;DR

This paper introduces ECKO, a new method for multivariate inference in high-dimensional fMRI data that clusters voxels and ensembles results to improve variable selection while controlling false discoveries.

Contribution

ECKO is a novel algorithm that applies knockoff inference on voxel clusters and uses ensembling to stabilize results in high-dimensional settings.

Findings

ECKO controls FDR at the nominal level in brain imaging data.

ECKO shows increased sensitivity compared to existing methods.

Empirical results demonstrate robustness and improved variable detection.

Abstract

Continuous improvement in medical imaging techniques allows the acquisition of higher-resolution images. When these are used in a predictive setting, a greater number of explanatory variables are potentially related to the dependent variable (the response). Meanwhile, the number of acquisitions per experiment remains limited. In such high dimension/small sample size setting, it is desirable to find the explanatory variables that are truly related to the response while controlling the rate of false discoveries. To achieve this goal, novel multivariate inference procedures, such as knockoff inference, have been proposed recently. However, they require the feature covariance to be well-defined, which is impossible in high-dimensional settings. In this paper, we propose a new algorithm, called Ensemble of Clustered Knockoffs, that allows to select explanatory variables while controlling the false discovery rate (FDR), up to a prescribed spatial tolerance. The core idea is that knockoff-based inference can be applied on groups (clusters) of voxels, which drastically reduces the problem's dimension; an ensembling step then removes the dependence on a fixed clustering and stabilizes the results. We benchmark this algorithm and other FDR-controlling methods on brain imaging datasets and observe empirical gains in sensitivity, while the false discovery rate is controlled at the nominal level.

Equations

$$\mathbf{y}=f(\mathbf{X}{\mathbf{w}}^{\ast })+\sigma \mathbf{ϵ},$$

$$\text{SNR}=\frac{{\Vert \mathbf{X}{\mathbf{w}}^{\ast }\Vert }_2^2}{{\sigma }^2{\Vert \mathbf{ϵ}\Vert }_2^2}.$$

$$\text{FDP}=\frac{\#\{k\in \hat{S}:k\notin S\}}{\#\{k\in \hat{S}\}}$$

$$\delta \text{-FDP}=\frac{\#\{k\in \hat{S}:\forall j\in S,d(j,k)>\delta \}}{\#\{k\in \hat{S}\}}$$

$$\begin{array}{rl}\text{FDR}& =\mathbb{E}[\text{FDP}],\\ \delta \text{-FDR}& =\mathbb{E}[\delta \text{-FDP}].\end{array}$$

$$\underset{\mathbf{w}\in {\mathbb{R}}^{2p}}{\min }\frac{1}{2}\Vert \mathbf{y}-[\mathbf{X},\tilde{\mathbf{X}}]\mathbf{w}{\Vert }_2^2+\lambda \Vert \mathbf{w}{\Vert }_1,$$

$$p_j=\frac{\#\{k:z_k\le -z_j\}}{p}.$$

$$\forall j\in [p],p_j=\text{quantile-aggregation}(\{p_j^{(b)}:b\in [B]\})$$

Full text

1 Introduction

Medical images are increasingly used in predictive settings, in which one wants to classify patients into disease categories or predict some outcomes of interest. Besides predictive accuracy, a fundamental question is that of opening the black box, *i.e. *understanding the combinations of observations that explains the outcome. A particular relevant question is that of the importance of image features in the prediction of an outcome of interest, conditioned on other features. Such conditional analysis is a fundamental step to allow causal inference on the implications of the signals from image regions in this outcome; see e.g. [12] for the case of brain imaging. However, the typical setting in medical imaging is that of high-dimensional small-sample problems, in which the number of samples $n$ is much smaller than the number of covariates $p$. This is further aggravated by the steady improvements in data resolution. In such cases, classical inference tools fail, both theoretically and practically. One solution to this problem is to reduce the massive number of covariates by utilizing dimension reduction, such as clustering-based image compression, to reduce the number of features to a value close to n; see e.g. [4]. This approach can be viewed as the bias/variance trade-off: some loss in the localization of the predictive features —bias— is tolerated as it comes with less variance —hence higher power— in the statistical model. This is particularly relevant in medical imaging, where localizing predictive features at the voxel level is rarely important: one is typically more interested in the enclosing region.

However, such a method suffers from the arbitrariness of the clustering step and the ensuing high-variance in inference results with different clustering runs, as shown empirically in [6]. [6] also introduced an algorithm called Ensemble of Clustered Desparsified Lasso (ECDL), based on the inference technique developed in [13], that provides p-values for each feature, and controls the Family Wise Error Rate (FWER), i.e. the probability of making one or more false discoveries. In applications, it is however more relevant to control the False Discovery Rate (FDR) [3], which indicates the expected fraction of false discoveries among all discoveries, since it allows to detect a greater number of variables. In univariate settings, the FDR is easily controlled by the Benjamini-Hochberg procedure [3], valid under independence or positive correlation between features. It is unclear whether this can be applied to multivariate statistical settings. A promising method which controls the FDR in multivariate settings is the so-called knockoff inference [2, 5], which has been successfully applied in settings where $n\approx p$. However, the method relies on randomly constructed knockoff variables, therefore it also suffers from instability. Our contribution is a new algorithm, called Ensemble of Clustered Knockoffs (ECKO), that i) stabilizes knockoff inference through an aggregation approach; ii) adapts knockoffs to $n\ll p$ settings. This is achieved by running the knockoff inference on the reduced data and ensembling the ensuing results.

The remainder of our paper is organized as follows: Section 2 establishes a rigorous theoretical framework for the ECKO algorithm; Section 3 describes the setup of our experiments with both synthetic and brain imaging data predictive problems, to illustrate the performance of ECKO, followed by details of the experimental results in Section 4; specifically, we benchmark this approach against the procedure proposed in [7], that does not require the clustering step, yet only provides asymptotic ($n\rightarrow\infty$) guarantees. We show the benefit of the ECKO approach in terms of both statistical control and statistical power.

2 Theory

2.1 Generalized Linear Models and High Dimensional Setting

Given a design matrix $\mathbf{X}\in\mathbb{R}^{n\times p}$ and a response vector $\mathbf{y}\in\mathbb{R}^{n}$, we consider that the true underlying model is of the following form:

[TABLE]

where $\mathbf{w}^{*}\in\mathbb{R}^{p}$ is the true parameter vector, $\sigma\in\mathbb{R}^{+}$ the (unknown) noise magnitude, $\bm{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{n})$ the noise vector and $f$ is a function that depends on the experimental setting (e.g. $f=Id$ for the regression problem or e.g. $f=sign$ for the classification problem). The columns of $\mathbf{X}$ refer to the explanatory variables also called features, while the rows of $\mathbf{X}$ represent the coordinates of different samples in the feature space. We focus on experimental settings in which the number of features $p$ is much greater than the number of samples $n$ i.e. $p\gg n$. Additionally, the (true) support denoted by $S$ is given by $S=\{k\in[p]:\mathbf{w}_{k}^{*}\neq 0\}$. Let $\hat{S}$ denotes an estimate of the support given a particular inference procedure. We also define the signal-to-noise ratio (SNR) which allows to assess the noise regime of a given experiment:

[TABLE]

A high SNR means the signal magnitude is strong compared to the noise, hence it refers to an easier inference problem.

2.2 Structured Data

In medical imaging and many other experimental settings, the data stored in the design matrix $\mathbf{X}$ relate to structured signals. More precisely, the features have a peculiar dependence structure that is related to an underlying spatial organization, for instance the spatial neighborhood in 3D images. Then, the features are generated from a random process acting on this underlying metric space. In our paper, the distance between the $j$-th and the $k$-th features is denoted by $d(j,k)$.

2.3 FDR control

In this section, we introduce the false discovery rate (FDR) and a spatial generalization of the FDR that we called $\delta$-FDR. This quantity is important since a desirable property of an inference procedure is to control the FDR or the $\delta$-FDR. In the following, we assume that the true model is the one defined in Section 2.1.

Definition 1

False discovery proportion (FDP).* Given an estimate of the support $\hat{S}$ obtained from a particular inference procedure, the false discovery proportion is the ratio of the number selected features that do not belong to the support (false discoveries) divided by the number of selected features (discoveries):*

[TABLE]

Definition 2

$\delta$-FDP.* Given an estimate of the support $\hat{S}$ obtained from a particular inference procedure, the false discovery proportion with parameter $\delta>0$, denoted $\delta$-FDP is the ratio of the number selected features that are at a distance more than $\delta$ from any feature of the support, divided by the number of selected features:*

[TABLE]

One can notice that for $\delta=0$, the FDP and the $\delta$-FDP refer to same quantity i.e. $0\text{\emph{-FDP}}=FDP$.

Definition 3

False Discovery Rate (FDR) and $\delta$-FDR.* The false discovery rate and the false discovery rate with parameter $\delta>0$ which is denoted by $\delta$-FDR are respectively the expectations of the FDP and the $\delta$-FDP:*

[TABLE]

2.4 Knockoff Inference

Initially introduced by [2] to identify variables in genomics, the knockoff filter is an FDP control approach for multivariate models. This method has been improved to work with mildly high-dimensional settings in [5], leading to the so-called model-X knockoffs:

Definition 4

Model-X knockoffs [5].* The model-X knockoffs for the family of random variables $\mathbf{X}=(X_{1},\dots X_{p})$ are a new family of random variables $\tilde{\mathbf{X}}=(\tilde{X}_{1},\dots,\tilde{X}_{p})$ constructed to satisfy the two properties:*

For any subset $\mathcal{K}\subset\left\{1,\dots,p\right\}$: $(\mathbf{X},\tilde{\mathbf{X}})_{\text{swap}(\mathcal{K})}\stackrel{{\scriptstyle d}}{{=}}(\mathbf{X},\tilde{\mathbf{X}})$,

where the vector $(\mathbf{X},\tilde{\mathbf{X}})_{\text{swap}(\mathcal{K})}$ denotes the swap of entries $X_{j}$ and $\tilde{X}_{j}$, $\forall j\in\mathcal{K}$ 2. 2.

\tilde{\mathbf{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}}}\mathbf{y}\mid\mathbf{X}* where $\mathbf{y}$ is the response vector.*

In a nutshell, knockoff procedure first creates extra null variables that have a correlation structure similar to that of the original variables. A test statistic vector is then calculated to measure the strength of the original versus its knockoff counterpart. An example of such statistic is the lasso-coefficient difference (LCD) that we use in this paper:

Definition 5

Knockoff procedure with intermediate p-values [2, 5].**

Construct knockoff variables, produce matrix concatenation: $[\mathbf{X},\mathbf{\tilde{X}}]\in\mathbb{R}^{n\times 2p}$. 2. 2.

Calculate LCD by solving

[TABLE]

and then, for all $j\in[p]$, take the difference $z_{j}=\left\lvert\hat{\mathbf{w}}_{j}(\lambda)\right\rvert-\left\lvert\hat{\mathbf{w}}_{j+p}(\lambda)\right\rvert$. 3. 3.

Compute the p-values $p_{j}$, for $j\in[p]$:

[TABLE] 4. 4.

Derive q-values by Benjamini-Hochberg procedure: $(q_{j})_{j\in[p]}=\text{BHq}\left((p_{j})_{j\in[p]}\right)$ 5. 5.

Given a desired FDR level $\alpha\in(0,1)$: $\hat{S}=\{j:q_{j}\leq\alpha\}$.

Remark 1

The above formulation is distinct from that of [2, 5], but it is equivalent. We use it to introduce the intermediate variables $p_{j}$ for all ${j\in[p]}$.

Our first contribution is to extend this procedure computing $q_{j}$ by aggregating different draws of knockoffs before applying the Benjamini-Hochberg (BHq) procedure. More precisely, we first compute $B$ draws of knockoff variables and, using Equation 6, we derive the corresponding p-values $p_{j}^{(b)}$, for all $j\in[p]$ and $b\in[B]$. Then, we aggregate them for each $j$ in parallel, using the quantile aggregation procedure introduced in [11]:

[TABLE]

We then proceed with the fourth and fifth steps of the knockoff procedure described in Def. 5.

2.5 Dimension reduction

Knockoff (KO) inference is intractable in high-dimensional settings, as knockoff generation requires the estimation and inversion of covariance matrices of size $(2p\times 2p)$. Hence we leverage data structure by introducing a clustering step that reduces data dimension before applying KO inference. As in [6], assuming the features’ signals are spatially smooth, it is relevant to consider a spatially-constrained clustering algorithm. By averaging the features with each clustering, we reduce the number of parameters from $p$ to $q$, the number of clusters, where $q\ll p$. KO inference on cluster-based signal averages will be referred to as Clustered Knockoffs (CKO). However, it is preferable not to fully rely on a particular clustering, as a small perturbation on the input data has a dramatic impact on the clustering solution. We followed the approach used in [9] that aggregates solutions across random clusterings. More precisely, they build $C$ different clusterings from $C$ different random subsamples of size $\lfloor 0.7n\rfloor$ from the full sample $\mathbf{X}$, but always using the same clustering algorithm.

2.6 The Ensemble of Clustered Knockoff Algorithm

The problem is to aggregate the q-values obtained across CKO runs on different clustering solutions. To do so, we transfer the q-values from clusters (group of voxels) to features (voxels): given a clustering solution $c\in[C]$, we assign to each voxel the q-value of its corresponding cluster. More formally, if, considering the $c$-th clustering solution, the $k$-th voxel belongs to the $j$-th cluster denoted by $G^{(c)}_{j}$ then the q-value $\tilde{q}^{(c)}_{k}$ assigned to this voxel is: $\tilde{q}^{(c)}_{k}=q_{j}^{(c)}$ if $k\in G^{(c)}_{j}$. This procedure hinges on the observation that the FDR is a resolution-invariant concept —it controls the proportion of false discoveries. In the worst case, this results in a spatial inaccuracy of $\delta$ in the location of significant activity, $\delta$ being the diameter of the clusters. Finally, the aggregated q-value $\tilde{q}_{k}$ of the $k$-th voxel is the average of the q-values $\tilde{q}^{(c)}_{k}$, $c\in[C]$, received across C different clusterings: given the FDR definition Equation 5, FDPs can naturally be averaged. The algorithm is summarized in Alg. 1 and represented graphically in Figure 1.

2.7 Theoretical Results

Ensemble of Clustered Knockoffs (ECKO).

Theorem 1

$\delta$-FDR control by the Ensemble of Clustered Knockoffs procedure at the voxel level.* Assuming that the true model is the one defined in Equation 1, using the q-values $\tilde{q}_{k}$ defined in Section 2.6, the estimated support $\hat{S}=\{k:\tilde{q}_{k}\leq\alpha\}$ ensure that the $\delta$-FDR is lower than $\alpha$.*

Sketch of the proof

(details are omitted for the sake of space). We first establish that the aggregation procedure yields q-values $q_{j}^{c}$ that control the FDR. This follows simply from the argument given in the proof of Theorem 3.1 in [11]. Second, we show that broadcasting the values from clusters ($q$) to voxels ($\tilde{q}$) still controls the FDR, yet with a possible inaccuracy of $\delta$, where $\delta$ is the supremum of clusters diameters: the $\delta$-FDR is controlled. This comes from the resolution invariance of FDR and the definition of $\delta$-FDR. Third, averaging-based aggregation of the q-values at the voxel level, controls the $\delta$-FDR. This stems from the definition of the FDR as an expected value.

2.8 Alternative approaches

In the present work, we use two alternatives to the proposed CKO/ECKO approach: the ensemble of clustered desparsified lasso (ECDL) [6] and the APT framework from [7]. As we already noted, ECDL is structured as ECKO. The main differences are that it relies on desparsified lasso rather than knockoff inference and returns p-values instead of q-values. The APT approach was proposed to return feature-level p-values for binary classification problems (though the generalization to regression is straightforward). It directly works at the voxel level, yet with two caveats:

• •

Statistical control is granted only in the $n\rightarrow\infty$ limit

• •

Unlike ECDL and ECKO, it is unclear whether the returned score represents marginal or conditional association of the input features with the output.

For both ECDL and APT, the returned p-values are converted to q-values using the standard BHq procedure. The resulting q-values are questionable, given that BHq is not valid under negative dependence between the input q-values [3]; on the other hand, practitioners rarely check the hypothesis underlying statistical models. We thus use the procedure in a black-box mode and check its validity a posteriori.

3 Experiments

3.0.1 Synthetic data.

To demonstrate the improvement of the proposed algorithm, we first benchmark the method on 3D synthetic data set that resembles a medical image with compact regions of interest that display some predictive information. The size of the weight vector $\mathbf{w}$ is $50\times 50\times 50$, with 5 regions of interest (ROIs) of size $6\times 6\times 6$. A design matrix $\mathbf{X}$ that represents random brain signal is then sampled according to a multivariate Gaussian distribution. Finally, the response vector $\mathbf{y}$ is calculated following linear model assumption with Gaussian noise, which is configured to have $SNR\approx 3.6$, similar to real data settings. An average precision-recall curve of 30 simulations is calculated to show the relative merits of single cluster Knockoffs inference versus ECKO and ECDL and APT. Furthermore, we also vary the Signal-to-Noise Ratio (SNR) of the simulation to investigate the accuracy of FDR control of ECKO with different levels of difficulty in detecting the signal.

3.0.2 Real MRI dataset.

We compare single-clustered Knockoffs (CKO), ECKO and ECDL on different MRI datasets downloaded from the Nilearn library [1]. In particular, the following datasets are used:

• •

Haxby [8]. In this functional-MRI (fMRI) dataset, subjects are presented with images of different objects. For the benchmark in our study, we only use the brain signal and responses for images related to faces and houses of subject 2 ($n=216,p=24083$).

• •

Oasis [10]. The original collection include data of gray and white matter density probability maps for 416 subjects aged 18 to 96, 100 of which have been clinically diagnosed with very mild to moderate Alzheimer’s disease. The purpose for our inference task is to find regions that predict the age of a subject ($n=400,p=153809$).

We chose $q=500$ in all experiments for the algorithms that require clustering step (KO, ECKO and ECDL). In the two cases, we start with a qualitative comparison of the returned results. The brain maps are ternary: all regions outside $\hat{S}$ are zeroed, while regions in $\hat{S}$ get a value of $+1$ or $-1$, depending on whether the contribution to the prediction is positive or negative. For ECKO, a vote is performed to decide whether a voxel is more frequently in a cluster with positive or negative weight.

4 Results

4.0.1 Synthetic data.

A strong demonstration of how ECKO makes an improvement in stabilizing the single-clustering Knockoffs (CKO) is shown in Fig. 2. There is a clear distinction between selection of the orange area at lower right and the blue area at upper right in the CKO result, compared to the ground truth. Moreover, CKO falsely discovers some regions in the middle of the cube. By Contrast, ECKO’s selection is more similar to the true 3D weight cube. While it returns a wider selection than ECKO, ECDL also claims more false discoveries, most visibly in the blue area on upper-left corner. At the same time, APT returns adequate results, but is more conservative than ECKO.

Fig. 3(a) is the result of averaging 30 simulations for the 3D brain synthetic data. ECKO and ECDL obtain almost identical precision-recall curve: for a precision of at least 90%, both methods have recall rate of around 50%. Meanwhile, CKO falls behind, and in fact it cannot reach a precision of over 40% across all recall rates. APT yields the best precision-recall compromise, slightly above ECKO and ECDL.

When varying SNR (from $2^{-1}$ to $2^{5}$) and investigating the average proportion of false discoveries ($\delta$-FDR) made over the average of 30 simulations (Fig. 3(b)), we observe that CKO fails to control $\delta$-FDR at nominal level 10% in general. Note that accurate $\delta$-FDR control would be obtained with larger $\delta$ values, but this makes the whole procedure less useful. The ECDL controls $\delta$-FDR at low SNR level. However, when the signal is strong, ECDL might select more false positives. ECKO, on the other hand, is always reliable —albeit conservative— keeping FDR below the nominal level even when SNR increases to larger magnitude.

Oasis & Haxby dataset. When decoding the brain signal on subject 2 of the Haxby dataset using response vector label for watching ’Face vs. House’, there is a clear resemblance of selection results between ECKO and ECDL.

Using an FDR threshold of 10%, both algorithms select the same area (with a difference in size), namely a face responsive region in the ventral visual cortex, and agree on the sign of the effect. However, on Oasis dataset, thresholding to control the FDR at 0.1 yields empty selection with ECDL and APT, while ECKO still selects some voxels. This potentially means that ECKO is statistically more powerful than ECDL and APT.

5 Conclusion

In this work, we proposed an algorithm that makes False Discovery Rate (FDR) control possible in high-dimensional statistical inference. The algorithm is an integration of clustering algorithm for dimension reduction and aggregation technique to tackle the instability of the original knockoff procedure. Evaluating the algorithm on both synthetic and brain imaging datasets shows a consistent gain of ECKO with respect to CKO in both FDR control and sensitivity. Furthermore, empirical results also suggest that the procedure achieves non-asymptotic statistical guarantees, yet requires the $\delta$-relaxation for FDR.

The number of clusters represents a bias-variance trade-off: increasing it can reduce the bias (in fact, the value of $\delta$), while reducing it improves the conditioning for statistical inference, hence the sensitivity of the knockoff control. We set it to 500 in our experiments. Learning it from the data is an interesting research direction.

We note that an assumption of independence between hypothesis tests is required for the algorithm to work, which is often not the case in realistic scenarios. Note that this is actually the case for all FDR-controlling procedures that rely on the BHq algorithm. As a result, making the algorithm work with relaxed assumption is a potential direction for our future study. Furthermore, the double-aggregation procedure makes the algorithm more expensive, although it results in embarrassingly parallel loops. An interesting challenge is to reduce the computation cost of this procedure. Another avenue to explore for the future is novel generative schemes for knockoff, based e.g. on deep adversarial approaches.

Acknowledgement This research is supported by French ANR (project FASTBIG ANR-17-CE23-0011) and Labex DigiCosme (project ANR-11-LABEX-0045-DIGICOSME). The authors would like to thank Sylvain Arlot and Matthieu Lerasle for fruitful discussions and helpful comments.