Ensemble-based Ultrahigh-dimensional Variable Screening

1 Introduction

The evolution of data acquisition technologies and computing power has allowed researchers nowadays to collect and store data with high dimensionality and complex structure much more efficiently. Examples can be found in gene expression microarray data, single nucleotide polymorphism (SNP) data, magnetic resonance imaging (MRI) data, highfrequency financial data, and others. One common task is to extract useful variables from a high dimensional feature space to explain or predict a response variable. Traditional variable selection methods such as forward selection, backward elimination and best subset selection become computationally expensive or even infeasible at these conditions. To address these problems, a family of penalized least squares based methods has been developed. Examples include Lasso and Adaptive Lasso ([Tibshirani, 1996; Zou, 2006]), SCAD [Fan and Li, 2001], elastic net [Zou and Hastie, 2005], and MCP [Zhang, 2010]. However, when the dimensionality p is much larger than the sample size n or even grows exponentially with n, the aforementioned penalization methods can perform poorly or even become infeasible due to the simultaneous challenges of computational expediency, statistical accuracy and algorithm stability [Fan et al., 2009]. For example, in MRI studies, images with dimension 1024×1024×200 can be acquired for each subject, and due to the high cost of the MRI scanning, studies might only contain less than 100 subjects. If we treat the signal from each voxel as a feature, the dimension of feature space p is much higher than the sample size n.

∗These authors contributed equally to this work.

†Correspondence Author.

A natural idea to address these challenges is to reduce the dimensionality p from a large scale to a relatively large scale d using a fast screening algorithm, and then the ultrahighdimensional problem can be greatly simplified into a moderately high-dimensional one. Subsequently, standard penalized variable selection methods can be applied to the remaining variables. Fan and Lv [2008] first introduced the sure independence screening (SIS) by ranking the marginal correlation of each covariate and the response. The good numerical performance and novel theoretical properties have made SIS popular in ultrahigh dimensional reduction. As a result, SIS and its extensions have been generalized to many important settings including generalized linear model [Fan et al., 2010], multi-index semi-parametric models [Zhu et al., 2011], nonparametric regression [Fan et al., 2011], quantile regression [He et al., 2013]. Other marginal screening methods based on different measure of association between predictors and response have also been studied, such as Kendall's τ [Li et al., 2012a], distance correlation [Li et al., 2012b]. We refer to [Liu et al., 2015] for a more comprehensive list of references.

Most feature screening methods are designed for specific models, and they all enjoy good theoretical and numerical performances under a set of conditions on the data generating process. However, in practice, we rarely know the actual relationship between features and the response, and different covariates might have different relationships with the response.

The heterogeneity presented here motivates us to consider various screening methods simultaneously in practice, so the natural question here is how to aggregate the results from these screening methods to achieve more robust and accurate performances. To the authors' best knowledge, there is no existing work on addressing these problems. In this paper, we propose an ensemble-based framework for combining results from multiple variable screening methods.

It is well known that ensemble methods is usually more accurate than a single learner, and they have already achieved great success in many real-world tasks [Zhou, 2012]. In the area of variable selection, the idea of ensemble has been studied before ([Bach, 2008; Bol´on-Canedo and AlonsoBetanzos, 2019]). In section 2, we introduce the details of ensemble-based variable screening. In section 3, we establish the theoretical properties of proposed method. In section 4, we conduct extensive simulation studies under different models. In section 5, we conduct a real data analysis using functional magnetic resonance imaging (fMRI) data to predict attention deficit hyperactivity disorder (ADHD) status. In summary, we have the following contributions and findings:

1. We propose an ensemble-based variable screening framework to efficiently combine results from different variable screening methods. The proposed framework is very flexible, and can be paralleled.
2. We prove that the proposed ensemble-based variable screening method inherits the nice theoretical properties (e.g. consistency, sure screening property) of the base candidates.
3. We conduct extensive simulation studies and real data analysis to illustrate the numerical performance of proposed framework. The proposed ensemble-based method more robust again model specification than using single variable screening method, that is to say, even though the proposed ensemble based screener may not outperform all base screener in a specific model setting, but it has the most consistent and robust performance across different model settings.

2 Ensemble-based Variable Screening

2.1 AGeneral Framework

We consider the problem of variable screening in ultrahighdimensional feature space, where we observe response variable Y and the associated covariate vector X = (X1, . . . ,Xp)T . Consider the conditional distribution function of Y given x, denoted by F(y|X) = P(Y < y|X). Define two sets of variables:

If j ∈ A, Xj is referred to as an active feature, otherwise an inactive feature. The goal is to reduce dimensionality p from a large scale to a moderate scale by a fast and efficient method, while including all the variables in set A.

The key idea of the marginal screening procedure is to rank all predictors by using a utility measure between the response and each predictor and then to retain the top variables for further investigation. It usually involves three steps:

1. Calculate the screening statistic vector F ⊤ = [f1, . . . , fp]. For example, fj = cor(Xj, Y ) is used in SIS.
2. Obtain the rank vector R⊤ = [r1, . . . , rp] by ranking F from largest to smallest (usually a large fj indicates a stronger relationship between Xj and Y).
3. Choose the ⌊γp⌋ top ranked features as the active set ̂ A, and γ ∈ (0, 1) is predetermined. To be specific, the vote vector V ⊤ = [v1, . . . , vp] with vj = I(rj ≤ ⌊γp⌋).

Different screening methods usually result in different ranking, and in the next section we introduce the concept of variable screening ensembles to combine the results from different methods.

2.2 Constructing Variable Screening Ensembles

Generally, an ensemble is constructed in two steps: generating the base learners and combining them. To get a good ensemble, it is generally believed that the base learners should be as accurate as possible, and as diverse as possible. Consider K base variable screening candidates, each one of them returns screening utilities of all p variables. Now we introduce the variable screening ensemble matrix:

where W⊤ i = [wi1, . . . ,wip], i = 1, . . . ,K denotes the screening utilities obtained using screening candidate i. Wi can be the vector of the screening statistics vector Fi of all p features, for example, the correlation between X and Y in SIS.

If a single base variable screening candidate has been used K times through bootstrapping or parameter tuning, we obtain homogeneous ensembles. On the other hand, if multiple screening candidates have been used, this leads to heterogeneous ensembles. Since the screening statistics produced by different base candidates can have different scale and/or range, Wi andWj may not be comparable. In this case, we can choose the rank of all p features Ri for Wi obtained using base candidate i. Another choice is to use the binary vote Vi introduced in the last section. Intuitively, using the rank R will be more robust compared to directly using screening statistic F , and it will be more efficient compared to using binary vote V since V can be sensitive to the choice of threshold parameter γ. In section 4, we further use numerical studies to confirm this intuition. After obtaining the ensemble matrix, we discuss how to efficiently combine the columns of the matrix.

2.3 Aggregation Functions

After constructing the ensemble matrix, the next step would be finding appropriate combination functions to aggregate the results. A general combination function is a multivariate function that projects the j-th column of the ensemble matrix to a real number, that is,

There are many combination functions available in the literature [Zhou, 2012]. In our ensemble screening framework, we consider two major types of combination functions.

Mean Combination: Taking means of prediction results is a commonly adopted approach in machine learning literature [Nguyen et al., 2018]. In our approach, mean combination function combines the output of the screening algorithms by taking the mean of each column in the ensemble matrix. Provided the screening utilities associated with the j-th column as w1j , . . . , wKj , the mean ensemble function is defined as

Median Combination: Instead of choosing the mean of w1j , . . . , wKj , we use the median, which is more robust when those screening utilities are skewed or with outliers.

As an alternative to mean, median provides a more robust result when the distribution of wij is skewed. For example, we use the ranking vector R⊤ i obtained from base candidate i as each row of the ensemble matrix, W⊤ i . For active variable Xj, if most base candidates return a small rank except for one or two "outliers", then the mean-aggregated result will be largely influenced by these "outliers", while the medianaggregated result will still be robust. Formally we have

2.4 Variables Selected

Denote the combined information of the ensemble matrix by applying mean or median aggregation function as {E1, . . . ,Ep}. For any given threshold γ ∈ (0, 1), the ⌊γp⌋ top ranked predictors are selected as the active set:

3 Technical Results

As variable screening only serves as the first step in high dimensional data analysis, the most important property as far as practical application is concerned the sure screening property [Fan and Lv, 2008]. That is, with probability approaching 1, the screening algorithm keeps all of the true active variables.

We require that the sure screening property holds for each screening algorithms. For the quantile based variable screening algorithms, e.g. QaSIS , we require the sure screening property holds at each quantile level τk which means the selected variables contain the true active set Ak with a probability tending to one.

Regarding the screening utilities, we require all the screening algorithms in the ensemble have consistency of the screening utilities, that is for 1 ≤ i ≤ K,

where wij is the corresponding screening utility for screening algorithm i and variable j, ̂ wij is an estimator for wij, and δn is some threshold number that is usually related to n. This requirement is reasonable as it is shown in most of the variable screening literatures [Fan and Lv, 2008; Fan et al., 2011; He et al., 2013; Zhu et al., 2011].

Lemma 3.1 (Consistency of aggregated screening utilities) Given number of K screening algorithms which are based on screening utilities wij, i = 1, . . . ,K and j = 1, . . . , p. Denote f as the combination function. Assume the following:

in which δn is some threshold constant related to n. We have the consistency of the aggregated screening utility which is:

Theorem 3.2 (Sure screening property) Denote wj as the jth aggregated screening utility, and w∗ j as the sample estimate. Denote A as the active variable set. Assume min j∈Awj > 2δn and our screening method select the variables with w∗ j > δn. Given the assumptions of the previous lemma, we have the sure screening property for our ensemble screening approach:

in which δn is some threshold number related to n.

Remark: To guarantee the sure screening property holds for our ensemble procedure, we just need δn to be the smallest among all the thresholds we combined. As long as all of the candidate methods show sure screening property, given the assumptions in Lemma 3.1 hold, our ensemble approach also enjoys sure screening property.

Lemma 3.3 (Lower bound for simultaneous screening probability) Denote ̂ T1 and ̂ T2 as two selected sets by applying two different screening algorithms. Define ̂ T simu = ̂ T1 ∩ ̂ T2 as the variable set selected by both of the screening algorithms. Define ̂Π1K , ̂Π2K , ̂Πsimu K as the probabilities of selected sets of screening algorithm 1, screening algorithm 2 and the simultaneous set containing a variable set K, K ⊆ {1, . . . , p}. Then we have:̂Πsimu K ≥ 2min { ̂Π1K , ̂Π2K } - 1.

Lemma 3.4 Let K ⊂ {1, . . . , p} be a set of variables and ̂ Ti be the set of selected variables by applying a variable screening algorithm i. If

Theorem 3.5 (Error control) Denote S = |T | and N = |F| as the number of underlying true important and unimportant variables. Correspondingly, denote ̂ S = |T ∩ ̂ T simu |

and ̂ N = |F ∩ ̂ F simu | as the number of estimated important and unimportant variables. In addition denote V = E(|F ∩ ̂ T simu |) as the expected number of falsely selected variables in ̂ T simu . Assume exchangeability which is P(k ∈ ̂ T ) = E( ̂ N)/N, where k ∈ (1, . . . , p). Also assume that the candidate variable screening process is not worse than random guessing. Given the screening threshold Tn and the threshold of selection probability πthr, we have:

Remark: If a variable is important, the underlying true ranking of the variable is higher than the other unimportant variables. With a threshold Tn, the true important variables are supposed to rank within Tn. If Tn is a relative small number say 30, and the probability threshold πthr is decently greater than 50%, given p = 1000 the expectation number of falsely discovered variables could be controlled within a small number. Intuitively, when each candidate screening method is good enough and the threshold Tn is small, the number of falsely selected variables could be controlled at a very low level.

4 Numerical Studies

4.1 Simulation Settings

In this section, we demonstrate the numerical performance of our proposed approach. For the base screening candidates, we consider three popular methods, the SIS in [Fan and Lv, 2008], the sure independence ranking and screening (SIRS) in [Zhu et al., 2011], and the quantile adaptive sure independence screening (QaSIS) in [He et al., 2013]. The SIS is based on linear model, and both SIRS and QaSIS are model free. These three methods all have their advantages and disadvantages, and none of them uniformly performs better than the others in all models. For QaSIS, different quantile levels can be used for the screening, we first consider the ensemble of QaSIS with different quantile levels, and there will be a more homogeneous scenario. Secondly we consider a more heterogeneous case by combining all these three different methods.

We consider two evaluation metrics: the first one is R, that is the smallest number of features that we need to include to ensure that all the active variables are selected; the second one (denoted by S) is the proportion of the active predictors selected when the threshold Tn = ⌊n/ log(n)⌋ is adopted. An effective variable screening procedure is expected to have the value of R reasonably small comparing to the number of active variables and the value of S close to one. All the experiments have been repeated 100 times, and the median and interquartile range (IQR) have been reported for both R and S.

Depending on the choice of Wi in the ensemble matrix and the aggregation function, different variable screening ensembles can be obtained. We use VSrE (Variable Screening Ensemble), VSrE-R and VSrE-V to denote the ensembles obtained by choosing W as the original screening statistic F, the rank R and the binary vote V respectively (defined in section 2.1). For the aggregation function, mean and median are used. For example, mean-based ensemble is denoted as VSrE(mean). Note that the median aggregation function for the vote based VSrE may not be a good choice since median is not a good summary statistics for binary data.

4.2 Models Settings

We consider the following four models, from linear, nonlinear to heteroscedastic cases. The random error term ε follows a standard normal distribution in all four models and is independent of X.

Model 1(n = 200, p = 2000) [Fan and Lv, 2008]: This is a model of the form Y = Xβ + ε, in which all the columns of X are generated from a standard normal distribution. There are in total eight true active predictors which are generated as follows: set a = 4 log(n)/n 12 and the coefficients corresponding to the active predictors are derived by (-1)u (a+ |z|), where u follows a Bernoulli distribution with p = 0.4 and z is drawn from a standard normal distribution.

Model 2(n = 200, p = 2000) [Zhu et al., 2011]: The random data is generated from Y = 2(X1 + 0.8X2 + 0.6X3 + 0.4X4 + 0.2X5) + exp(X20 + X21 + X22) · ε, where X = (X1, X2, . . . , X2000) follows a multivariate normal distribution with mean 0 and covariance Σ = (σij)p×p with σij = ρ|i-j| . We consider ρ = 0.8 here. The number of active variables here is 8.

Model 3(n = 400, p = 1000) [Fan et al., 2011]: First, define the following functions g1(x) = x; g2(x) = (2x - 1)2 ; g3(x) = sin(2πx)/(2 - sin(2πx)); g4(x) = 0.1 sin(2πx) + 0.2 cos(2πx) + 0.3 sin(2πx)2 + 0.4 cos(2πx)3 + 0.5 sin(2πx)3 . The random data are generated from: Y = 5g1(X1)+3g2(X2)+4g3(X3)+6g4(X4)+ √1.74ε, where X = (X1,X2, . . . ,X1000) has the similar structure as in Model 2 with ρ = 0.4.

Model 4(n = 400, p = 5000) [He et al., 2013]: Y = 2(X2 1 +X2 2 ) + [10-1 exp(X1 + X2 + X18 + X19 + . . . + X30)] · ε., where X = (X1,X2, . . . ,X5000) has the similar structure as Model 2 with ρ = 0.8. The number of active variables here is 15.

4.3 Combining Quantile Levels

We first consider VSrE constructed by applying QaSIS separately on different quantile levels. In order to show the effect of different choices of quantile sets, we show the results of applying four different quantile sets which have 5, 10, 19 and 37 equally spaced quantile levels respectively as follows: T1 = {0.1, 0.3, . . . , 0.9}, T2 = {0.05, 0.15, . . . , 0.95}, T3 = {0.05, 0.1, . . . , 0.95}, T4 = {0.05, 0.075, . . . , 0.95}. As an illustration of the choice of VSrEs and combination functions, all five types of methods are applied. Model 1 is used for this simulation. The line chart in Figure 1 shows the comparison of single screening method (QaSIS(0.5) and QaSIS(0.75)) and the proposed ensemble-based methods using 4 different sets of quantile levels, while the boxplots in Figure 2 focus on the case when quantile level set T2 with K = 10 is used. From Figure 1 and Figure 2, we have the following observations:

1. All proposed ensemble-based methods outperform the single base screening method, with regard to both R and S.
2. The performance of ensemble-based method improves with the increase of the number of base screening methods K. But the improvement becomes substantial as K reaches a certain level, while the computation time increases linearly with K.
3. When the screening statistic is comparable across ensemble candidates, both mean- and median- based VSrE

and VSrE-R provides reasonably good results.

4.4 Combining Different Variable Screening Methods

Instead of constructing a VSrE by applying QaSIS on different quantile levels, we could apply different screening methods to the sample and combine the results. In the following, we show the power of simply combining three commonly used variable screening methods which are SIS, SIRS and QaSIS. Quantile level 0.75 is used for QaSIS. As the screening utilities of different methods may have very different magnitudes, it is more appropriate to use VSrE-R in order to avoid the situation that some screening methods with relative large screening utilities might dominate the combined result.

Figure 3 compares the ensemble-based VSrE-R(mean) and VSrE-R(median) with the single screening candidates under Model 1 to Model 4. Figure 4 provides more detailed results on Model 2 about both R and S. Based on Figure 3 and Figure 4, we have the following observations:

1. The proposed VSrE-R(median) has the most robust performance across all four models.
2. Each ensemble candidate has its own advantage: SIS performs the best under Model 1, SIRS under Model 3 and QaSIS under Model 4. The ensemble-based methods inherit the good performance of the base screening candidates, while still robust to model misspecification as long as some of the base learners perform reasonably well.
3. For the aggregation function, median is a more robust choice.

5 Real Data Analysis

Attention deficit hyperactivity disorder (ADHD) is a brain disorder marked by an ongoing pattern of inattention and hyperactivity-impulsivity that interferes with functioning or development. Recent development of medical imaging such as functional magnetic resonance imaging (fMRI) and diffusion tensor imaging (DTI) shows promising potential in predicting patients outcomes and understanding the underlying pathophysiology of diseases [Greicius et al., 2007].

We use the ADHD-200 Consortium data which is a publicly available resting-state fMRI (rs-fMRI) data [Milham et al., 2012] in this study. fMRI measures brain activity by detecting changes associated with blood flow [Huettel et al., 2004]. This data set contains 120 subjects (n = 120) from the NYU site (New York University Child Study Center) of the ADHD-200 Consortium, 42 are typically developing children showing ADHD negative and 78 are diagnosed as ADHD. The Anatomical Automatic Labeling (AAL) atlas [TzourioMazoyer et al., 2002] was used for the parcellation.

In this study, our goal is to find which brain connectivity pair of ROIs (region of interest) is contributing to the ADHD express levels hence each brain connectivity is considered to be a variable. For each subject, there are 172 time points and the AAL has 116 ROIs. We obtain the mean time series for each of the 116 regions by averaging the fMRI time series over all voxels in the region, hence initially we have p = (116 × 116 - 116)/2 = 6670 predictors. Because of the large p small n scenario (n = 120, p = 6670), a variable screening procedure is necessary to remove some noise brain connectivities in order to apply lower dimensional variable selection approaches. We apply VSrE-R(median) by combining SIS, SIRS and QaSIS with quantile level 0.75 and use QaSIS at quantile level 0.75 as a comparison. In the variable selection step, we employ both LASSO and SCAD. The tuning parameters of LASSO and SCAD are selected. In the next step, we choose those ROI connectivities that are simultaneously selected by both LASSO and SCAD. For the classification, we adopt the support vector machine classifier (SVM) [Hearst et al., 1998] with linear kernel. All tuning parameter are selected by 10 folds cross validation.

Unlike the simulation settings, we do not know the true active variables, so we can not use metrics like R and S. Instead we measure the usefulness of the selected features by their predictive power. Figure 5 shows the prediction accuracy using 100 bootstrapped samples. The top panel is the histogram of accuracy using ensemble-based method while the bottom panel is the one using QaSIS with quantile level 0.75. We can see a substantial improvement using the proposed method.

Figure 5: Prediction accuracy of ADHD status using single screening candidate (QaSIS (0.75)) and proposed ensemble-based method

Also, some of the partial correlations screened out by the proposed ensemble-based method include regions associated with a network called the default-mode network, which is a large and robustly replicable network of brain regions that is associated with task-irrelevant mental processes and mindwandering. Similar findings are reported in [Uddin et al., 2008; Tian et al., 2006].

6 Discussion

In this paper, we introduce a general ensemble-based framework to combine results from different variable screening methods. The simulation studies confirm our intuition that the ensemble-based methods indeed inherit the good performance of the candidates, that is to say, as long as some candidates have decent performance, the ensemble-based method will work reasonably well. This is very important in practice, since when we do not have much information about the actual data generating process, the ensemble-based methods will be a more robust and safer choice than relying on a single method. Also, we find that the median rank-based aggregation method has the most robust performance. It is also worthy to point out that the proposed framework can be easily implemented in a parallel way so it will not add too much computation burden.