{'fraction_non_alphanumeric': 0.045625335480407946, 'fraction_numerical': 0.016835016835016835, 'mean_word_length': 4.5434135785772245, 'pattern_counts': {'":': 1, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
FDR-HS: An Empirical Bayesian Identification of Heterogenous Features in Neuroimage Analysis

Recent studies found that in voxel-based neuroimage analysis, detecting and differentiating"procedural bias"that are introduced during the preprocessing steps from lesion features, not only can help boost accuracy but also can improve interpretability. To the best of our knowledge, GSplit LBI is the first model proposed in the literature to simultaneously capture both procedural bias and lesion features. Despite the fact that it can improve prediction power by leveraging the procedural bias, it may select spurious features due to the multicollinearity in high dimensional space. Moreover, it does not take into account the heterogeneity of these two types of features. In fact, the procedural bias and lesion features differ in terms of volumetric change and spatial correlation pattern. To address these issues, we propose a"two-groups"Empirical-Bayes method called"FDR-HS"(False-Discovery-Rate Heterogenous Smoothing). Such method is able to not only avoid multicollinearity, but also exploit the heterogenous spatial patterns of features. In addition, it enjoys the simplicity in implementation by introducing hidden variables, which turns the problem into a convex optimization scheme and can be solved efficiently by the expectation-maximum (EM) algorithm. Empirical experiments have been evaluated on the Alzheimer's Disease Neuroimage Initiative (ADNI) database. The advantage of the proposed model is verified by improved interpretability and prediction power using selected features by FDR-HS.


Introduction
In recent years, the issue of model interpretability attracts an increasing attention in voxel-based neuroimage analysis of disease prediction, e.g. [9,5]. Examples include, but not limited to, the preprocessed features on structural Magnetic Resonance Imaging (sMRI) images that usually contain the following voxel-wise features: (1) lesion features that are contributed to the disease (2) procedural bias introduced during the preprocessing steps and shown to be helpful in classification [12,3] (3) irrelevant or null features which are uncorrelated with disease label. Our goal is to stably select non-null features, i.e. lesion features and procedural bias with high power/recall and low false discovery rate (FDR).
The lesion features have been the main focus in disease prediction. In dementia disease such as Alzheimer's Disease (AD), such features are thought to be geometrically clustered in atrophied regions (hippocampus and medial temporal lobe etc.), as shown by the red voxels in Fig. 1 (A). To explore such spatial patterns, multivariate models with Total Variation [10] regularization can be applied by enforcing smoothness on the voxels in neighbor, e.g. the n 2 GFL [16] can stably identify the early damaged regions in AD by harnessing the lesions.
Recently, another type of features called procedural bias, which are introduced during the preprocessing steps, are found to be helpful for disease prediction [12]. Again, taking AD as an example, the procedural bias refer to the mistakenly enlarged Gray Matter (GM) voxels surrounding locations with cerebral spinal fluid (CSF) spaces enlarged, e.g. lateral ventricle, as shown in Fig. 1 (A). This type of features has been ignored in the literature until recently, when the GSplit LBI [12] was targeted on capturing both types of features via a split of tasks of TV regularization (for lesions) and disease prediction with general linear model (with procedural bias). By leveraging such bias, it can outperform models which only focus on lesions in terms of prediction power and interpretability.
However, GSplit LBI may suffer from inaccurate feature selection due to the following limitations in high dimensional feature space: 7 : (1) multicollinearity: high correlation among features in multivariate models [14]; (2) "heterogenous features": the procedural bias and lesion features differ in terms of volumetric change (enlarged v.s. atrophied) and particularly spatial pattern (surroundingly distributed v.s. spatially cohesive). Specifically, the multicollinearity could select spurious null features which are inter-correlated with non-nulls. Moreover, GSplit LBI fails to take into account the heterogeneity since it enforces correlation on features without differentiation. Such problems altogether may result in inaccurate selection of non-nulls, especially procedural bias. As shown in Fig. 1 (B) and Table 2, the procedural bias selected by GSplit LBI are unstably scattered on regions that are less informative than ventricle. Moreover, the collinearity among features tends to select a subset of features among correlated ones, as discussed in [17]. Such a limitation leads to the ignorance of many meaningful regions (such as medial temporal lobe, thalamus etc.) of GSplit LBI in selecting lesion features, as identified by the purple frames of FDR-HS in Fig. 1 (B). 7 Please refer supplementary material for detailed and theoretical discussion Moreover, the two problems above may get worse as dimensionality grows. In our experiments with a fine resolution (4 × 4 × 4 of 20,091 features), the prediction accuracy of GSplit LBI deteriorates to 89.77% (as shown in Table 3), lower than 90.91% reported in [12] with a coarse resolution (8 × 8 × 8 of 2,527 features). To resolve the problems above, we propose a "two-groups" empirical Bayes method to identify heterogenous features, called FDR-HS standing for "FDR Heterogenous smoothing" in this paper. As a univariate FDR control method, it avoids the collinearity problem by proceeding voxel-by-voxel, as discussed in [7]. Moreover, it can deal with heterogeneity by regularizing on features with different levels of spatial coherence in different feature groups, which remedies the problem of losing spatial patterns that most conventional mass-univariate models suffer from, such as two sample T-test, BH q [4] and LocalFDR [7]. By introducing a binary latent variable, our problem turns into a convex optimization and can be solved efficiently via EM algorithm like [13]. The method is applied to a voxelbased sMRI analysis for AD with a fine resolution (4×4×4 of 20,091 features). As a result, our proposed method exhibits a much stabler feature extraction than GSplit LBI, and achieves much better classification accuracy at 91.48%.

Method
Our dataset consists of p voxels and N samples {x i , y i } N 1 where x ij denotes the intensity value of the j th voxel of the i th sample and y i = {±1} indicates the disease status (−1 denotes AD). The FDR-HS method is proposed to select non-null features. Such method is the combination of "two-groups" model and heterogenous regularization, which is illustrated in Fig. 2 and discussed below. Model Formulation. Assuming for each voxel i ∈ {1, ..., p}, the statistic z i is sampled from the following mixture: where s i is a latent variable indicating if the voxel i belongs to the group of null features (s i = 0) or the group of non-null ones (s i = 1), c i = p(s i = 1) = sigmoid(β i ) = e βi / 1 + e βi and z i = Φ −1 (F N −2 (t i )) with t i computed by twosample t-test. Correspondingly, f 0 (·) is density function of nulls, i.e. uncorrelated with AD and f 1 (·) is that of non-nulls, i.e. procedural bias and lesions. The loss function can thus be defined as negative log-likelihood of z i : which can be viewed as logistic regression (when f 0 and f 1 are replaced with binaries, as (2.6)) with identity design matrix since (2.1) proceeds voxel-by-voxel. Hence, it does not have the problem of multicollinearity. Selecting Features. To select features, we compute the posterior distribution of s i conditioned on z i and β i (estimated β i ) and features with are selected. The γ ∈ (0, 1) is pre-setting threshold parameter. Heterogenous Spatial Smoothing. However, (2.1) may lose spatial structure of non-nulls, especially lesion features. Besides, note that the procedural bias and lesion features are heterogenous in terms of volumetric change and level of spatial coherence. Hence, to capture the spatial structure of heterogenous features, we split the graph of voxels which denotes as G 8 into three subgraphs, i.e. G = G 1 ∪ G 2 ∪ G 3 with: where G 1 denotes the subgraph restricted on enlarged voxels (procedural bias since -1 denotes AD); G 2 denotes the subgraph restricted on degenerate voxels (lesion features); G 3 denotes the bipartite graph with the edges connecting enlarged and degenerate voxels. The optimization function can be redefined as: where D G k β = (i,j)∈E k β i − β j for k ∈ {1, 2, 3} denote graph difference operator on G k=1,2,3 . By setting the group of regularization hyper-parameters {λ pro , λ les , λ pro-les } with different values, we can enforce spatial smoothness on three subgraphs at different level in a contrast to the traditional homogeneous regularization in [13]. The choice of each hyper-parameter, similar to [13], it is a trade-off between over-fitting and over-smoothing. Too small value tends to select features more than needed, while too large value will oversmooth hence the features are less clustered. Note that lesion features are more spatially coherent than procedural bias and they are located in different regions, the reasonable choice of regularization hyper-parameters tend to have λ les ≤ λ pro ≤ λ pro-les .
Optimization. Note that the function (2.5) is not convex. Hence we adopted the same idea in [13] that introduced the latent variables s i and = 1 if . The (β) and g(β) are modified as: To solve (2.7), we can implement Expectation-Maximization (EM) algorithm to alternatively solve β and s. Suppose currently we are in the (k + 1) th iteration.
In the E-step, we can estimate s i by expectation value conditional on (β k , z i ): and expand (β|s k ) using a second-order Taylor approximation at the β k . Then the M-step turns into a generalized lasso problem with square loss: Note that X and ‹ D G are sparse matrices, hence (2.8) can be efficiently solved by Alternating Direction Method of Multipliers (ADMM) [6] which has a complexity of O(p log p). Estimation of f 0 and f 1 . Before the iteration, we need to estimate f 0 (z) and f 1 (z). The marginal distribution of z can be regarded as mixture models with p components: z ∼ 1 which is equivalent to LocalFDR [7]. We can therefore implement the CM (Central Matching) [7] method to estimate {f 0 (z),c} and kernel density to estimate f (z). The f 1 (z) can thus be given as

Experimental Results
In this section, we evaluate the proposed method by applying it on the ADNI database http://adni.loni.ucla.edu. The database is split into 1.5T and 3.0T (namely 15 and 30) MRI scanner magnetic field strength datasets. The 15 dataset contains 64 AD, 110 MCI (Mild Cognitive Impairment) and 90 NC, while the 30 dataset contains 66 AD and 110 NC. After applying DARTEL VBM [2] preprocessing pipeline on the data with scale of 4×4×4 mm 3 voxel size, there are in total 20,091 voxels with average values in GM population on template greater than 0.1 and they are served as input features. We designed experiments on 1.5T AD/NC, 1.5T MCI/NC and 3.0T AD/NC tasks, namely 15ADNC, 15MCINC and 30ADNC, respectively.

Prediction Results
To test the efficacy of selected features by FDR-HS and compare it with other univariate models (as listed in Table 1), we feed them into elastic net classifier, which has been one of the state-of-the-arts in the prediction of neuroimage data [11]. The hyper-parameters are determined by grid-search. In details, the threshold hyper-parameter of p-value in T-test and q-value in BH q are optimized through {0.001, 0.01, 0.02, 0.05, 0.1}; the threshold hyper-parameter for choosing non-nulls, i.e. γ for FDR-HS (2.3) and the counterpart of LocalFDR [7], are chosen from {0.1, 0.2, ..., 0.5}. Besides, the regularization parameters λ pro , λ les and λ pro-les of FDR-HS are ranged in {0.1, 0.2, ..., 2}. For elastic net, the regularization parameter is chosen from {0.1, 0.2, ..., 2, 5, 10}; the mixture parameter α is from {0, 0.01, ..., 1}. Moreover, we compare our model to GSplit LBI and elastic net, adopting the same optimized strategy for hyper-parameters in [12] (the top 300 negative voxels are identified as procedural bias [12]) and those of elastic net following after the univariate models, as mentioned above.
A 10-fold cross-validation strategy is applied and the classification results for all tasks are summarized in Table 1. As shown, our method yields better results than others in all cases, that includes: (1) FDR-HS can select features with more prediction power than other univariate models due to the ability to capture heterogenous spatial patterns; (2) FDR-HS can achieve better classification results than multivariate methods in high dimensional settings, in which the non-nulls may be represented by other nulls that are highly correlated with them.

Feature Selection Analysis
We used 2-d images of 30ADNC to visualize the features of all methods under the hyper-parameters that give the best accuracy. As shown in Fig. 3, the lesion features selected by FDR-HS are located clustered in early damaged regions; while procedural bias are surrounding around lateral ventricle. Besides, such a result is given by λ les < λ pro < λ pro-les , which agrees with that the larger value results in features with lower level of spatial coherence. In contrast, the lesions selected by T-test and BH q are scattered and redundant; some procedural bias around lateral ventricle are missed by BH q and LocalFDR. Moreover, GSplit LBI selected procedural bias on regions with CSF space less enlarged than lateral ventricle; besides, it ignored lesions located in medial temporal lobe, Thalamus and Fusiform etc., which are believed to be the early damaged regions [1,8].
Besides, we also evaluated the stability of selected features using multi-set Dice Coefficient (mDC) measurement defined in [16]. Larger mDC implies more stable feature selection. As shown in Table 2, our model can obtain more stable results than GSplit LBI which suffer the "collinearity" problem.

Conclusions
In this paper, a "two-groups" Empirical-Bayes model is proposed to stably and efficiently select interpretable heterogenous features in voxel-based neuroimage analysis. By modeling prior probability voxel-by-voxel and using a heterogenous regularization, the model can avoid multicollinearity and exploit spatial patterns of features. With experiments on ADNI database, the features selected by our models have better interpretability and prediction power than others.

Supplementary Information A Multicollinearity problem
One can implement sparse multivariate models for feature selection and classification by minimizing the penalized optimization function. Particularly, n 2 GFL [16] was proposed to stably capture the degenerate voxels by harnessing the sparsity and geometrically clustering properties under the β ≥ 0 constraint and regularization function Ω(β) = Dβ 1 with D = [I; ρD G ]. 9 On the basis of this, GSplit LBI was proposed to capture additional procedural bias via a variable splitting scheme. However, such multivariate models suffer from the multicollinearity problem, i.e. high correlation among features, including the following aspects in details: (1) the collinearity between non-nulls and nulls can make the irrepresentable condition of genlasso which ensures the successfully recovery of the true support set [15] hard to satisfy. Such problem can violate the non-nulls to be selected; (2) the collinearity among non-nulls of sparsity model such as lasso/genlasso tends to only select one feature/region among correlated ones. (3) the procedural bias may be represented by other variables which are highly correlated with them due to the "collinearity" between non-nulls and nulls during minimization of General Linear Model (GLM). Specifically, the general penalized optimization function of the multivariate model is: ä −y i ·(β 0 +x T i β) with β 0 denoting the bias parameter and Ω(β) is the regularization function. Different choice of Ω(β) leads to different model. Since GSplit LBI and genlasso enforce same sparsity regularizations on lesion features, we firstly discuss the problem of genlasso under high dimensional space: (1) the irrepresentable condition which ensures the model selection consistency (i.e. successfully recover the true support set) is not easy to satisfy (2) only select one feature/cluster among correlated ones. To understand (1), note in [15] that under linear model the necessary condition for genlasso to satisfy model selection consistency is that IC 1 < 1. When D = I, the slightly stronger version of this condition is X T S c X S X T S X S −1 1 < 1 where S denotes the true support set (lesion voxels). Such irrepresentable condition implies the "decorrelated" property of S and S c . However, such condition is hard to satisfy under high dimensional space since covariates are more easier to be correlated.
To see (2), note in [17] that the lasso only selects one feature among a group of correlated ones. We claim that genlasso also suffers from this limitation. Specifically, it can be shown in Lemma 1 that the Total Variation regularization of genlasso tends to select only single region among a group of correlated ones. Lemma 1. Let G S1 = (V S1 , E S2 ) and G S2 = (V S2 , E S1 ) denotes two subgraphs (regions) that are not connected with other nodes in V , i.e.
Proof. It can be easily verified from the definition of A.1.
Now we discuss the problem of the procedural bias in GSplit LBI, i.e. the procedural bias may be represented by other variables which are highly correlated with them due to the "collinearity". Note that to select procedural bias, GSplit LBI adopts an variable splitting scheme Dβ −γ 2 2 by introducing an augmented variable γ, the loss function is redefined as: and implement it by the following iterative algorithm: where S k := supp(γ k ). Note that in A.3f, the β les , which is the projection of β onto the support set of γ, is the estimator to capture lesion features. The elements with comparably large magnitude among the remainder of such projection are regarded as procedural bias. However, under high dimensional data, such definition may suffer multicollinearity problem that the procedural bias are "submerged" or represented by other null variables that have high correlation with them. In detail, note that the procedural bias are less clustered than lesion features and also that in A.3b, when κ → ∞, α → 0, we have β k+1 → arg min β (β k 0 , β, γ k ), i.e. minimizing GLM model, it can then be shown in the following lemma that the algorithm may choose null variables while the procedural bias may not be successfully recovered.

B Heterogenous Features
We mentioned that the procedural bias and lesion features are heterogenous in terms of volumetric change (enlarged v.s. atrophied) and spatial patterns (surroundingly distributed v.s. spatially cohesive). The heterogeneity in terms of volumetric change is easy to understand, since the procedural bias commonly refer to enlarged GM voxels in voxel-based dementia analysis and lesion features refer to atrophied ones. To illustrate the heterogeneous levels of spatial coherence, we evaluate the edge density in 3D coordinate system by introducing 3D edge density (3dED) measurement for both selected lesions and procedural bias. In detail, 3dED is defined as: where S ± k denote the support set of lesion features and procedural bias in the k-th fold, respectively. For each fold k, we need to compute: For the graphs with p voxels which are embedded in 3-d coordinate space, the maximum number of edges is equal to: where R + = {x ∈ R, x ≥ 0}. Besides, the c 1 , c 2 , c 3 can be taken as length, width and height of a cube, r 1 r 2 + l are the remainder of p for a cube (c 1 , c 2 , c 3 ). The r 1 , r 2 can be taken as the length and width of the rectangle which is located on one of the surfaces of (c 1 , c 2 , c 3 ). The l is the residual of p − c 1 c 2 c 3 for rectangle (r 1 , r 2 ). According to B.1, the procedural bias are turned to be selected much less clustered than lesion features by all univariate models, e.g. the BH q yields 0.4677 for lesion features and 0.1621 for procedural bias; while localFDR has 0.4662 and 0.1699; FDR-HS has 0.5365 and 0.2535, which validates the heterogenous assumption in terms of the level of spatial coherence.