Feature Selection in High-Dimensional Models via EBIC with Energy Distance Correlation

Abstract

In this paper, the LASSO method with extended Bayesian information criteria (EBIC) for feature selection in high-dimensional models is studied. We propose the use of the energy distance correlation in place of the ordinary correlation coefficient to measure the dependence of two variables. The energy distance correlation detects linear and non-linear association between two variables, unlike the ordinary correlation coefficient, which detects only linear association. EBIC is adopted as the stopping criterion. It is shown that the new method is more powerful than Luo and Chen’s method for feature selection. This is demonstrated by simulation studies and illustrated by a real-life example. It is also proved that the new algorithm is selection-consistent.

1. Introduction

Advancements in technology have led to the production of sophisticated machines which are able to measure many details about every observational or experimental unit in a system. This results in data with more features (p or predictors) than the number of observational or experimental units (sample size n), referred to as high-dimensional data. Most of these data come from genetic research, e-commerce, biomedical imaging, and functional magnetic resonance imaging, among many others.

Since there are many more features (p) than the sample size (n), they are analyzed using a sparse high-dimensional regression (SHR) model:

$$y_i={\beta }_0+\sum _{j=1}^p{x}_{ij}{\beta }_j+{ϵ}_i,\mathrm{for}i=1,\dots ,n.$$

(1)

It is assumed that there is only a relatively small number of the nonzero ${\beta }_j$’s. The main goal in their analysis is feature selection. As stated by [1], feature selection typically has two goals. The first is for model building using desirable prediction properties. The second is for identifying the features with nonzero coefficients. For convenience, such features are referred to as relevant features in this paper.

One approach to the SHR model is to estimate the ${\beta }_j$’s by a regularization method, which is done by simultaneously minimizing the penalized least squares below:

$$\sum _{i=1}^n{\left(y_i-{\beta }_0-\sum _{j=1}^p{\beta }_j{x}_{ij}\right)}^2+\sum _{j=1}^p{p}_{\lambda }\left(\right|{\beta }_j\left|\right),$$

where $\lambda$ is the regulating parameter and $p_{\lambda }$ is a penalty function. When $p_{\lambda }$ is based on the $L_1$ norm, thus ${∥\beta ∥}_1={\sum }_{j=1}^p\left|{\beta }_j\right|$, it is referred to as the LASSO [2]. The $L_1$ norm penalty is able to shrink the coefficients of redundant predictors to zero. Thus, the LASSO usually results in sparse models that are easier to interpret. Other penalty functions such as the SCAD [3] and adaptive LASSO (ALasso) [4] have also been reported in the literature. SCAD smoothly clips the $L_1$ penalty (for small $|{\beta }_j|$), and assigns a constant penalty (for large $\mid {\beta }_j\mid$s). On the other hand, adaptive LASSO utilizes the minimax concave penalty [5], $p_{\lambda }(\mid {\beta }_j\mid )=\lambda w_j\mid {\beta }_j\mid$, where $w_j$ represents the weights. The regulating parameter, $\lambda$, is usually chosen using cross-validation (CV).

Another approach to analyzing the SHR model, large-p-small-n problem, is sequential variable selection, which is designed to reduce the dimension of the data such that $d<p$. There are two forms of sequential variable selection. The first uses the sure screening property, which selects from the many features a subset which contains the relevant features (predictors). This is usually followed by a regularization method such as SCAD or ALasso to identify and estimate the relevant predictors from the reduced feature space. The other form is to sequentially select the relevant features through a repetitive process which terminates when a stopping criterion is met.

A recent addition to sequential feature selection is the sequential LASSO cum EBIC in ultra high-dimensional feature space (SLasso) by [1], which sequentially solves partially penalized least squares problems and uses EBIC as the stopping criteria. The EBIC proposed by [6] are suitable for model selection in large model spaces. It has the ordinary BIC as a special case. For large model spaces, the ordinary BIC tends to select a model with many spurious variables. Let k and $k+1$ be the number of predictors in two models, respectively. Using EBIC as the selection or stopping criteria, the model with k predictors is selected if the EBIC(k + 1) > EBIC(k).

In SLasso, sequentially solving the partially penalized least squares reduces to selecting the feature(s) which maximize the ordinary correlation coefficient between the features and the response variable at each step. It is well-known that the Pearson correlation coefficient is used for measuring the strength of linear associations. Thus, maximizing the Pearson correlation coefficient might not work well for data structures where the relationship between at least one feature and the response variable is nonlinear.

In this article, we propose the use of the energy distance correlation instead of the ordinary correlation coefficient to identify and maximize both the linear and nonlinear relationships that might exist between each feature and the response. Energy distance is a metric that measures the distance between the distributions of random vectors. The name ‘energy’ is motivated by analogy to the potential energy between objects in a gravitational space. The potential energy is zero if and only if the locations (the gravitational centers) of the two objects coincide, and increases as their distance in space increases [7]. The energy distance correlation has an explicit relationship with the product-moment correlation, but unlike the classical definition of correlation, energy distance correlation is zero only if the random vectors are independent. The empirical energy distance correlation is based on Euclidean distances between sample elements rather than sample moments.

The remainder of the article is arranged as follows: in Section 2, we discuss the derivation of the energy distance correlation, extended Bayesian information criteria and our proposed method (energy distance correlation with EBIC (Edc + EBIC)). In Section 3, we report simulation studies comparing Edc + EBIC with various other methods and provide an analysis of real data. In Section 4, we conclude the article with a discussion of the results.

2. EBIC with Energy Distance Correlation

2.1. Energy Distance Correlation

The authors in [8] proposed the energy distance correlation between two random variables. Suppose that $W\in {\mathbb{R}}^p$ and $Z\in {\mathbb{R}}^q$ are two random vectors with $\mathbb{E}\parallel W\parallel <\infty$, and $\mathbb{E}\parallel Z\parallel <\infty$, where $\parallel .\parallel$ is the euclidean norm and $\mathbb{E}$ is the expected value. Let F and G be the cumulative distribution function (CDF) of W and Z, respectively. Further, let $W^{{}^{\prime }}$ denote an independent and identically distributed (iid) copy of W; that is, W and $W^{{}^{\prime }}$ are iid. Similarly, Z and $Z^{{}^{\prime }}$ are iid.

The squared energy distance can be defined in terms of expected distances between the random vectors

$$D^2(F,G):=2\mathbb{E}\parallel W-Z\parallel -\mathbb{E}\parallel W-W^{\prime }\parallel -\mathbb{E}\parallel Z-Z^{\prime }\parallel \ge 0,$$

and the energy distance between distributions F and G is defined as the square root of $D^2(F,G)$.

The energy distance correlation between random vectors W and Z with finite first moments is the nonnegative number $\mathcal{R}(W,Z)$ defined by

$$\mathcal{R}(W,Z)=\left\{\begin{array}{cc}\frac{{\nu }^2(W,Z)}{\sqrt{{\nu }^2\left(W\right){\nu }^2\left(Z\right)}},\hfill & {\nu }^2\left(W\right){\nu }^2\left(Z\right)>0\hfill \\ 0,\hfill & {\nu }^2\left(W\right){\nu }^2\left(Z\right)=0\hfill \end{array}\right.$$

where ${\nu }^2(W,Z)$ is the energy distance covariance between W and Z, ${\nu }^2\left(W\right)$ and ${\nu }^2\left(Z\right)$ are the energy distance variance of W and Z respectively.

For a statistical sample $(w,z)=\{(w_k,z_k),k=1,2,\dots ,n\}$ from a pair of real-valued or vector-valued random variables $(W,Z)$, the sample energy distance correlation, ${\mathcal{R}}_n(W,Z)$, is calculated by first computing the n by n distance matrices $\left(a_{j,k}\right)$ and $\left(b_{j,k}\right)$ containing all pairwise distances $\left(a_{j,k}\right)=\parallel W_j-W_k\parallel ,$ $j,k=1,2,\dots ,n$ and $\left(b_{j,k}\right)=\parallel Z_j-Z_k\parallel ,$ $j,k=1,2,\dots ,n$ where $\parallel .\parallel$ denotes euclidean norm. Secondly, calculate all doubly centered distances $A_{j,k}=a_{j,k}-{\overline{a}}_{j.}-{\overline{a}}_{.k}+{\overline{a}}_{..},B_{j,k}=b_{j,k}-{\overline{b}}_{j.}-{\overline{b}}_{.k}+{\overline{b}}_{..}$ where ${\overline{a}}_{j.}$ is the $j^{th}$ row mean, ${\overline{a}}_{.k}$ is the $k^{th}$ column mean, and ${\overline{a}}_{..}$ is the grand mean of the distance matrix of the w sample. The notation is similar for the b values.

The squared sample distance covariance (a scalar) is the arithmetic average of the products $A_{j,k}{B}_{j,k}.$

$${\nu }_n^2(w,z)={\displaystyle \frac{1}{n^2}}\sum _{j=1}^n\sum _{k=1}^n{A}_{j,k}{B}_{j,k}.$$

The sample energy distance variance for sample w

$${\nu }_n^2\left(w\right)={\displaystyle \frac{1}{n^2}}\sum _{j,k=1}^n{A}_{j,k}^2.$$

The sample energy distance variance for sample z

$${\nu }_n^2\left(z\right)={\displaystyle \frac{1}{n^2}}\sum _{j,k=1}^n{B}_{j,k}^2.$$

The sample energy distance correlation is

$${\mathcal{R}}_n(W,Z)=\left\{\begin{array}{cc}\frac{{\nu }_n^2(W,Z)}{\sqrt{{\nu }_n^2\left(W\right){\nu }_n^2\left(Z\right)}},\hfill & {\nu }_n^2\left(W\right){\nu }_n^2\left(Z\right)>0\hfill \\ 0,\hfill & {\nu }_n^2\left(W\right){\nu }_n^2\left(Z\right)=0.\hfill \end{array}\right.$$

Some basic properties of the distance correlation are as follows:

(i)

$0\le R_n(W,Z)\le 1;$

(ii)

If $E\left(\right|W{|}_p+\left|Z{|}_q\right)<\infty ,$ then $R_n(W,Z)=0$ if and only if W and Z are independent;

(iii)

Suppose that $R_n(W,Z)=1$. Then, there exist a vector a, a nonzero real number b and an orthogonal matrix C such that $Z=a+bWC.$

For further details on the energy distance correlation, see [7].

2.2. EBIC

Ref. [6] derived the EBIC which have special cases as AIC and BIC. Let $\{(y_i,x_i):i=1,2,\dots ,n\}$ be independent observations. Suppose that the conditional density function of $y_i$ given $x_i$ is $f\left(y_i\right|x_i,\beta ),$ where $\beta \in \mathsf{\Theta }\subset R^{p_n},p_n$ being a positive integer. The likelihood function of $\beta$ is given by

$$L_n\left(\beta \right)=f(x;\beta )=\prod _{i=1}^{n}f\left(y_i\right|x_i,\beta ).$$

Denote $y={(y_1,y_2,\dots ,y_n)}^{\tau }$. Let $s\subset \{1,2,\dots ,p_n\}$ and $\beta \left(s\right)$ be the parameter vector $\beta$ with those components outside s set to 0. Let S be the underlying model space, i.e., $S=\{s:s\subseteq \{1,2,\dots ,p_n\}\},$ let $p\left(s\right)$ be a prior for model s. Assume that, given s, the prior density of $\beta \left(s\right)$ is $\pi \left(\beta \right(s\left)\right).$ The posterior is

$$p\left(s\right|y)=\frac{m\left(y\right|s\left)p\right(s)}{{\sum }_{s\in S}m\left(y\right|s)p\left(s\right)},$$

where $m\left(y\right|s)$ is the likelihood in model s, i.e.,

$$m(Y|s)=\int f(y;\beta (s))\pi (\beta (s))d\beta (s).$$

Suppose S is partitioned into ${\cup }_{j=1}^p{S}_j,$ such that models within each $S_j$ have an equal dimension. Let $\tau \left(S_j\right)$ be the size of $S_j$. Assign the prior distribution $P\left(S_j\right)$ proportional to ${\tau }^{\eta }\left(S_j\right)$ for some $\eta$ between 0 and 1. For each $s\in S_j$, assign equal probability, $p\left(s\right|S_j)=1/\tau \left(S_j\right)$; this is equivalent to $P\left(s\right)$ for $s\in S_J$ proportional to ${\tau }^{-\gamma }\left(S_j\right)$, where $\gamma =1-\eta .$ Then, the extended BIC family is given by

$$EBI{C}_{\gamma }\left(s\right)=-2log{L}_n\left\{\widehat{\beta }\left(s\right)\right\}+\left|s\right|log\left(n\right)+2\gamma ln\left(\tau \left(S_{\left|s\right|}\right)\right),0\le \gamma \le 1,$$

where $\widehat{\beta }\left(s\right)$ is the maximum likelihood estimator of $\beta \left(s\right)$ and $\left|s\right|$ is the number of components in s.

2.3. Energy Distance Correlation with EBIC (Edc + EBIC) Algorithm

We propose a sequential model selection method which we call energy distance correlation with EBIC, and for convenience abbreviate it as Edc + EBIC. Let $y_i,i=1,\dots ,n$ be a continuous response variable and $x_j,j=1,\dots ,p$ be an $n\times p$ data matrix. Let S be the index set of all predictors. Let $s_0=\{j:{\beta }_j\ne 0,j=1,\dots ,p\}$. For $s\subset S$, let $s^{-}=s^c\cap s_0$. If $s\subset s_0$, then $s^{-}$ is the complement of s in $s_0$. Let $p_0=\left|s_0\right|$ be the number of elements in the set $s_0$.

At the initial stage we standardize all the variables. Next, we find the energy distance correlation between the response variable and each of the predictor variables—$\{\mathcal{R}(x_j,y)j=1,\dots ,p.\}$. We then select the predictor (feature) which has the highest distance correlation with the response and store it in the active set $s_{*1}$.

Let $\mathcal{L}\left(s\right)$ be the linear space spanned by the columns of $X\left(s\right)$ and $H\left(s\right)$ its corresponding projection matrix, i.e., $H\left(s\right)=X\left(s\right){\left[X^{\tau }\left(s\right)X\left(s\right)\right]}^{-1}{X}^{\tau }\left(s\right)$. Next, we compute $I-H\left(s_{*1}\right)$, $\mathrm{EBIC}\left(s_{*1}\right)$, $\tilde{y}=[I-H\left(s_{*k}\right)]y\mathrm{and}{\tilde{x}}_j=[I-H\left(s_{*k}\right)]x_j.$ The variable $\tilde{y}$ is the unexplained part of y by $X\left(s_{*1}\right)$. This gives $X\left(s_{*1}\right)$ close to a zero chance of being selected in the subsequent steps.

For the general step where $k>1$, we calculate $\{\mathcal{R}({\tilde{x}}_j,\tilde{y})j=1,\dots ,p.\}$ and update the active set to $s_{*k+1}$, which is the union of all the previous selected variables and the current one. We then compute $\mathrm{EBIC}\left(s_{*k+1}\right)$ and compare it with $\mathrm{EBIC}\left(s_{*k}\right)$. The procedure stops if $\mathrm{EBIC}\left(s_{*k+1}\right)>\mathrm{EBIC}\left(s_{*k}\right)$. The selected variables which we call the relevant variables will be $X\left(s_{*k}\right)$. We can then fit a linear regression model between the response y and the relevant variables.

We wish to note that care must be taken in fitting this model because some of the predictors might be non-linearly related to y, and thus some of the predictors may have to enter into the model in their quadratic or cubic form, etc. Alternatively, a Box–Cox transformation can be performed on the data before fitting the model.

The algorithm details are given in the following.

• Initial Step: With $y,x_j,j=1,\dots ,p$ satisfying $y^{\tau }1=0,x_j^{\tau }1=0$ and $y^{\tau }y=n,x_j^{\tau }{x}_j=n$, compute $\mathcal{R}(x_j,y)$ for $j\in S$. Let

  $$s_{TEMP}=\{j:\mathcal{R}(x_j,y)=\underset{j^{\prime }\in S}{max}\mathcal{R}(x_{j^{\prime }},y)\}.$$
  Let $s_{*1}=s_{TEMP}$ be the active set. Compute $I-H\left(s_{*1}\right)$ and $\mathrm{EBIC}\left(s_{*1}\right),$ where $H\left(s\right)=X\left(s\right){\left[X^{\tau }\left(s\right)X\left(s\right)\right]}^{-1}{X}^{\tau }\left(s\right).$
• General Step: In the selection step k, compute $\mathcal{R}({\tilde{x}}_j,\tilde{y})$ for $j\in s_{*k}^c,$ where $\tilde{y}=[I-H\left(s_{*k}\right)]y,{\tilde{x}}_j=[I-H\left(s_{*k}\right)]x_j.$ Let

  $$s_{TEMP}=\{j:\mathcal{R}({\tilde{x}}_j,\tilde{y})=\underset{j^{\prime }\in S_{*k}^c}{max}\mathcal{R}({\tilde{x}}_j,\tilde{y})\}.$$
  Let $s_{*k+1}=s_{*k}\cup s_{TEMP}$. Compute $\mathrm{EBIC}\left(s_{*k+1}\right)$. If $\mathrm{EBIC}\left(s_{*k+1}\right)>\mathrm{EBIC}\left(s_{*k}\right)$, stop; otherwise, continue computing $I-H\left(s_{*k+1}\right)$.
• When the process terminates, return the least-squares estimates for parameters in the selected model.

2.4. Selection Consistency of Edc + EBIC

We attempt to establish the large sample property for the Edc + EBIC. We will show that under regular conditions, the Edc + EBIC is selection-consistent. The proof essentially follows the approach in [9]. We proceed with the regularity conditions.

Assumption 1. 

Random vectors X and Y possess the subexponential tail probabilities, uniformly in p, specified as follows. There is a constant $a_0>0$, such that for any $0<a\le 2a_0,$${sup}_p{max}_{1\le k\le p}E\{exp\left(a{\parallel X_k\parallel }_1^2\right)\}<\infty$ and $E\{exp\left(a{\parallel Y\parallel }_q^2\right)\}<\infty$.

Assumption 2. 

The minimum distance correlation of predictors on which y functionally depends satisfies ${min}_{j\in s_0}\mathcal{R}({\tilde{X}}_j,\tilde{Y})\ge 2c{n}^{-d},$ for some constants $0<c<1$ and $0\le d<1/2$.

Assumption 3. 

For the index set S of all predictors, let $s_0=\{j:{\beta }_j\ne 0,j=1,\dots ,p\}$ and $p_0=\left|s_0\right|$ ($p_0$ is the number of elements in the set $s_0$). For $s\subset S$ let $s^{-}=s^c\cap s_0$. If $s\subset s_0$ then $s^{-}$ is the complement of s in $s_0$. For $s\subset s_0$, ${max}_{j\in s_0^c}\mathcal{R}({\tilde{X}}_j,\tilde{Y})<q{max}_{j\in s^{-}}\mathcal{R}({\tilde{X}}_j,\tilde{Y})$ for some $0<q<1$, where $\tilde{Y}=[I-H\left(s_{*k}\right)]Y,{\tilde{X}}_j=[I-H\left(s_{*k}\right)]X_j.$ For $k=0$, $s_{*0}$ is defined as the empty set ∅.

Details for requiring Assumption 1 and 2 are stated in [9]. Intuitively, Assumption 1 is required to make it easy to establish a relationship between the energy distance correlation and the squared Pearson correlation to aid with the derivations in the proof. Assumption 2 requires that the energy distance correlation for the relevant predictors cannot be too small. Assumption 3 requires that the maximum energy distance correlation between the selected features and the residual response $\tilde{Y}$ is smaller than the maximum energy distance correlation between the remaining features and the residual response in the sequential step of the algorithm.

Theorem 1. 

Suppose that Assumptions 1–3 hold. The proposed Edc + EBIC with the energy distance correlation is consistent, i.e.,

$$\underset{n\to \infty }{lim}P(s_{*k^{*}}=s_{0n})=1,$$

where $s_{*k^{*}}$ is the set of features selected at the $k^{*th}$ step of Edc + EBIC such that $|s_{*k^{*}}|=p_{0n,}$$s_{0n}$ is the set of relevant features and $p_{0n}=\left|s_{0n}\right|$.

Proof. 

Suppose that $X\in {\mathbb{R}}^p$ and $Y\in {\mathbb{R}}^q$ with cumulative distribution function (CDF) F and G, respectively, where $\mathbb{E}\parallel X\parallel <\infty$, and $\mathbb{E}\parallel Y\parallel <\infty$. The energy distance correlation $\mathcal{R}(X,Y)$ is the square root of the standardized coefficient:

$$\mathcal{R}(X,Y)=\left\{\begin{array}{cc}\frac{{\nu }^2(X,Y)}{\sqrt{{\nu }^2\left(X\right){\nu }^2\left(Y\right)}},\hfill & {\nu }^2\left(X\right){\nu }^2\left(Y\right)>0\hfill \\ 0,\hfill & {\nu }^2\left(X\right){\nu }^2\left(Y\right)=0\hfill \end{array}\right.$$

where $0\le \mathcal{R}(X,Y)\le 1$. In the numerator is the distance covariance defined by [8], as

$$dco{v}^2(x,y)=S_1+S_2-2S_3,$$

where $S_j,j=1,2,\mathrm{and}3$ are defined as:

$$\begin{array}{cc}\hfill S_1=& \mathbb{E}\parallel X-X^{\prime }\parallel \parallel Y-Y^{\prime }\parallel \hfill \\ \hfill S_2=& \mathbb{E}\parallel X-X^{\prime }\parallel \mathbb{E}\parallel Y-Y^{\prime }\parallel \hfill \\ \hfill S_3=& \mathbb{E}\parallel X-X^{\prime }\parallel \parallel Y-Y^{″}\parallel \hfill \end{array}$$

(2)

where $(X,Y),(X^{\prime },Y^{\prime })$, and $(X^{″},Y^{″})$ are independently and identically distributed.

For a random sample $\{(x_i,y_i),i=1,\dots ,n\}$ from $(x,y)$, [8] estimated $S_1,S_2,S_3$ as:

$$\begin{array}{cc}\hfill \widehat{S_1}=& \frac{1}{n^2}\sum _{k,l=1}^n|x_k-x_l{|}_p{|y_k-y_l|}_q\hfill \\ \hfill \widehat{S_2}=& \frac{1}{n^2}\sum _{k,l=1}^n|x_k-x_l{|}_p\frac{1}{n^2}\sum _{k,l=1}^n{|y_k-y_l|}_q\hfill \\ \hfill \widehat{S_3}=& \frac{1}{n^3}\sum _{k=1}^n\sum _{l,m=1}^n|x_k-x_l{|}_p{|y_k-y_l|}_q\hfill \end{array}$$

so the sample distance covariance is ${\widehat{dcov}}^2=\widehat{S_1}+\widehat{S_2}-2\widehat{S_3}$.

The remaining part of the proof is to show that the energy distance correlation is uniformly consistent and has the sure screening property. The numerator and denominator of the energy distance correlation are similar, so to show the uniform consistency of the energy distance correlation it suffices to show that both the numerator and the denominator are uniformly consistent.

The uniform consistency of the numerator, ${\widehat{dcov}}^2=\widehat{S_1}+\widehat{S_2}-2\widehat{S_3}$, of the energy distance correlation between the random vectors $(x,y)$ is shown by [9]. However, in the general step of the sequential algorithm for Edc + EBIC, the energy distance correlation is calculated between the residuals $\tilde{y}=[I-H\left(s_{*k}\right)]y,\mathrm{and}{\tilde{x}}_j=[I-H\left(s_{*k}\right)]x_j$ at each step of the algorithm. Thus, to show the uniform consistency of Edc+EBIC it is equivalent to follow the proof by [9].

Additionally, in [9] they showed that the energy distance correlation has the sure screening property. They showed that the energy distance is able to select a subset of the features which contains the relevant features. Their argument applies here because we used the energy distance correlation as well, thus the Edc + EBIC has the sure screening property.

Therefore, the Edc + EBIC is selection-consistent since it is uniformly consistent and has the sure screening property. The proof is complete. □

3. Simulation Studies and Data Analysis

3.1. Sure Independence Screening Using Energy Distance Correlation

We establish the need for Ebc + EBIC by firstly examining the performance of a sure independence screening method introduced by [9] called Distance Correlation Sure Independence Screening (DC-SIS). This is similar to the Sure Independence Screening (SIS) introduced by [10].

In SIS, they perform a componentwise regression between each predictor and the response and select the first $n-1\mathrm{or}[n/log(n\left)\right]$ predictors with the largest estimates. Performing a componentwise regression is equivalent to finding the ordinary correlation between the response and each predictor when the two variables are standardized. Hence, in DC-SIS, they replaced the ordinary correlation with the energy distance correlation.

We examine the performance of DC-SIS through a simulation study. We are interested in observing, on average, the model size selected by SCAD or ALasso if we screened the data first using DC-SIS. We present two simulation set-ups. For each simulation we generated two hundred datasets, and for each dataset we ran SCAD, ALasso, DC-SIS + SCAD, DC-SIS + ALasso and found the average model size and the standard deviation.

In [10], details of two simulation setups we adapted for this subsection are discussed, namely independent features setup and dependent features setup. In

Table 1. Comparing model size selected with or without screening for $n=200,s=8,p=1000$.

Methods	MSize (SD)
SCAD	12.87 (7.292)
DC-SIS + SCAD	10.7 (3.1575)
ALasso	25.24 (9.0365)
DC-SIS + ALasso	11.74 (4.419)

and

Table 2. Comparing model size selected with or without screening for $n=800,s=14,p=3000$.

Methods	MSize (SD)
SCAD	16.62 (2.78807)
DC-SIS + SCAD	16.69 (3.5525)
ALasso	14.78 (0.7860)
DC-SIS + ALasso	14.78 (3.8522)

we present results under the independent features setup and in

Table 3. Comparing model size selected with or without screening for $n=200,p=1000,s=5$.

Methods	MSize (SD)
SCAD	12.215 (12.2560)
DC-SIS + SCAD	7.335 (2.6109)
ALasso	44.485 (13.8041)
DC-SIS + ALasso	8.21 (2.5844)

,

Table 4. Comparing model size selected with or without screening for $n=200,p=1000,s=8$.

Methods	MSize (SD)
SCAD	14.625 (10.3324)
DC-SIS + SCAD	10.905 (2.3158)
ALasso	19.95 (3.2789)
DC-SIS + ALasso	12.36 (3.9875)

and

Table 5. Comparing model size selected with or without screening for $n=800,p=3000,s=14$.

Methods	MSize (SD)
SCAD	19.185 (5.0146)
DC-SIS + SCAD	17.675 (4.6038)
ALasso	38.125 (5.6372)
DC-SIS + ALasso	31.845 (13.1011)

we present results under the dependent features setup. In each simulation, n is the sample size, p is the number of features and s is the true model size. For the screening using the energy distance correlation we chose $d=[n/logn]$ features and applied SCAD or ALasso.

In Table 1 and Table 2, we report the average selected model size and their standard deviations. We observe that applying the sure screening by distance correlation before either SCAD or ALasso in all cases did not lead to significant differences in the average model size when SCAD and ALasso were applied directly to the data. This suggests that either applying distance correlation before SCAD or ALasso did not yield the intended result, and thus needs some improvement.

In Table 3, Table 4 and Table 5 we report the selected model size and the standard deviation. We observe that applying DC-SIS followed by either SCAD or ALasso did not yield any significant difference in the average model size, as was also observed in the independent features setup.

3.2. Simulation Studies to Compare Edc + EBIC with Other Feature Selection Methods

In this simulation study we adopted two simulation setups from [1], which they call group A and group B, respectively. Under their group A we considered four settings of the covariance structure for the design matrix X, namely GA1, GA2, GA3, and GA5. In their group B setup we considered all three settings of the covariance structure for the design matrix X, namely GB1, GB2, and GB3. We compared the performance of adaptive LASSO (ALasso) [11], SCAD [12], SIS+SCAD [10], SLasso [1], and the energy distance correlation with EBIC (Edc + EBIC) based on the model size (MSize), positive discovery rate (PDR), $\mathrm{PDR}=\frac{|s_{*k^{*}}\cap s_o|}{|s_o|}$, and false discovery rate, $\mathrm{FDR}=\frac{|s_{*k^{*}}\cap s_0^c|}{|s_{*k^{*}}|}$ averaged over 200 and 500 simulations, respectively.

We considered the diverging pattern $(n,p,p_o)=(n,\left[5e^{n^{0.3}}\right],\left[4n^{0.16}\right])$, meaning that as the sample size increased, the number of predictors increased and the number of relevant predictors also increased. The coefficients were generated as independent random variables distributed as ${(-1)}^u(4n^{-0.15}+\left|z\right|),$ where $u\sim PBernoulli\left(0.4\right)$ and z is a normal random variable with mean 0 and satisfies $P\left(\right|z|\ge 0.1)=0.25.$ The variance of the error term in the linear model was determined by

$$h=\frac{{\beta }^{\tau }\Sigma \beta }{{\beta }^{\tau }\Sigma \beta +{\sigma }^2}=0.8$$

where $\Sigma$ is the variance-covariance matrix of relevant features. The response variable is simulated from the sparse high-dimensional regression (SHR) model

$$y_i={\beta }_0+\sum _{j=1}^p{\beta }_j{x}_{ij}+{ϵ}_i,i=1,...,n$$

3.2.1. Group A Simulations and Results

We used two sample sizes, $n=100$ and $n=200$. By the diverging pattern considered for the simulation, for a sample size of 100 we have $(n,p,p_0)=(100,268,8)$ and for a sample size of 200 we have $(n,p,p_0)=(200,672,9)$. Under the two sample sizes we had eight and nine relevant features, respectively, and we expected a well-performed selection model to select the right number of relevant predictors (model size) on average. The details about the covariance structure for GA1, GA2, GA3, and GA5 are in [1].

In

Table 6. We compared the methods using PDR, FDR, and model size (MSize) averaged over 200 simulation replications. The relevant predictors were 8 and the sample size was 100. The standard deviations are in parentheses.

Setting	Methods	MSize	PDR	FDR
GA1	ALasso	34.105 (13.95)	1.000 (0.000)	0.721 (0.120)
	SCAD	25.735 (5.020)	1.000 (0.000)	0.676 (0.065)
	SIS + SCAD	8.100 (1.790)	0.866 (0.239)	0.157 (0.167)
	SLasso	8.565 (0.848)	1.000 (0.000)	0.058 (0.081)
	Edc + EBIC	8.365 (1.375)	0.978 (0.125)	0.056 (0.085)
GA2	ALasso	34.455 (11.095)	1.000 (0.000)	0.754 (0.108)
	SCAD	25.650 (6.720)	0.876 (0.141)	0.709 (0.075)
	SIS + SCAD	7.335 (1.740)	0.813 (0.182)	0.103 (0.107)
	SLasso	6.055 (1.725)	0.688 (0.185)	0.080 (0.104)
	Edc + EBIC	6.075 (1.713)	0.717 (0.195)	0.050 (0.082)
GA3	ALasso	14.710 (3.847)	1.000 (0.000)	0.423 (0.131)
	SCAD	26.27 (5.244)	1.000 (0.000)	0.680 (0.070)
	SIS + SCAD	8.165 (1.160)	0.951 (0.113)	0.062 (0.078)
	SLasso	8.625 (1.005)	1.000 (0.000)	0.062 (0.089)
	Edc + EBIC	8.265 (1.373)	0.976 (0.132)	0.048 (0.074)
GA5	ALasso	23.845 (7.005)	0.964 (0.057)	0.652 (0.092)
	SCAD	24.070 (6.147)	0.997 (0.020)	0.642 (0.102)
	SIS + SCAD	7.605 (2.020)	0.832 (0.245)	0.127 (0.141)
	SLasso	7.650 (2.182)	0.856 (0.217)	0.089 (0.113)
	Edc + EBIC	7.180 (2.453)	0.842 (0.270)	0.050 (0.081)

and

Table 7. We compared the methods using PDR, FDR, and model size (MSize) averaged over 200 simulation replications. The relevant predictors were 9 and the sample size was 200. The standard deviations are in parentheses.

Setting	Methods	MSize	PDR	FDR
GA1	ALasso	27.670 (12.996)	1.000 (0.000)	0.638 (0.180)
	SCAD	17.035 (7.746)	1.000 (0.000)	0.454 (0.168)
	SIS + SCAD	9.215 (1.507)	1.000 (0.000)	0.112 (0.123)
	SLasso	8.710 (0.944)	1.000 (0.000)	0.072 (0.088)
	Edc + EBIC	8.42 (0.712)	1.000 (0.000)	0.045 (0.071)
GA2	ALasso	27.92 (9.686)	1.000 (0.000)	0.675 (0.124)
	SCAD	15.11 (6.241)	1.000 (0.000)	0.397 (0.171)
	SIS + SCAD	9.16 (1.509)	1.000 (0.000)	0.108 (0.171)
	SLasso	8.72 (0.947)	1.000 (0.000)	0.073 (0.089)
	Edc + EBIC	8.49 (0.763)	1.000 (0.000)	0.051 (0.076)
GA3	ALasso	27.115 (12.867)	1.000 (0.000)	0.632 (0.177)
	SCAD	16.245 (7.770)	1.000 (0.000)	0.434 (0.162)
	SIS + SCAD	9.22 (1.617)	1.000 (0.000)	0.110 (0.128)
	SLasso	8.70 (0.857)	1.000 (0.000)	0.072 (0.084)
	Edc + EBIC	8.47 (0.694)	1.000 (0.000)	0.050 (0.071)
GA5	ALasso	38.95 (8.308)	0.939 (0.075)	0.797 (0.054)
	SCAD	19.075 (7.427)	1.000 (0.000)	0.520 (0.159)
	SIS + SCAD	9.975 (1.858)	1.000 (0.000)	0.174 (0.132)
	SLasso	8.765 (1.125)	0.989 (0.061)	0.087 (0.094)
	Edc + EBIC	8.44 (0.768)	0.998 (0.025)	0.048 (0.073)

, we report the simulation results under the conditions for GA1, GA2, GA3, and GA5 using sample sizes 100 and 200, respectively. We observed that under all setups, Edc + EBIC improved in terms of the average model size, PDR, and FDR, as we increased the sample size. This demonstrates the selection consistency of Edc + EBIC.

3.2.2. Group B Simulations and Results

We considered three different covariance structures named GB1, GB2, and GB3 for the features (predictors), as used in [1]. We also increased the signal-to-noise ratio by increasing the value of the expected predictors.

GB1. All the features had constant pairwise correlation $p_{ij}=0.5.(n,p,p_0)=(100,200,15).$$\sigma =1.5.$ The coefficients of the relevant features were specified as $|{\beta }_j|=2.5$ for $1\le j\le 5,1.5$ for $6\le j\le 10,0.5$ for $11\le j\le 15.$ The signs of the coefficients were determined as ${(-1)}^{u_i}$, where the $u_i$s were iid Bernoulli random variables with probability of success $p=0.5.$

GB2. This structure was the same as in GB1, that is, $(n,p,p_0)=(100,200,15)$ and $\sigma =1.5.$ The covariance structure of the features was specified such that the partially orthogonality condition [11] was satisfied. Specifically, while $s_0$ was taken as $\{1,\dots 5,11,\dots ,15,21,\dots ,25\}$, the correlations were specified as ${\rho }_{ij}=0.5^{|i-j|}$ for $1\le i\le 215$ and $1\le j\le 215$. The coefficients were specified as $\left|\beta \right|=2.5$ for $1\le j\le 5,1.5$ for $10\le j\le 15,0.5$ for $21\le j\le 25.$ The signs of the coefficients were determined in the same way as in GB1.

GB3. $(n,p,p_0)=(100,1000,10)$ and $\sigma =1$. The relevant features were generated as iid standard normal variables with coefficients (3, 3.75, 4.5, 5.25, 6, 6.75, 7.5, 8.25, 9, 9.75). The irrelevant features were generated as

$$x_j=0.25Z_j+\sqrt{0.75}\sum _{k\in s_0}{X}_k,j\notin s_0,$$

where $Z_{j}s$ are iid standard normal and independent from the relevant features.

In

Table 8. We compared the methods using PDR, FDR, and model size (MSize) averaged over 500 simulation replications. The standard deviations are in parentheses.

Setting	Methods	MSize	PDR	FDR
GB1	ALasso	23.32 (3.018)	0.766 (0.062)	0.501 (0.066)
	SCAD	14.08 (1.644)	0.853 (0.065)	0.085 (0.068)
	SIS + SCAD	10.656 (1.688)	0.694 (0.112)	0.025 (0.067)
	SLasso	14.916 (2.194)	0.893 (0.081)	0.092 (0.089)
	Edc + EBIC	14.094 (2.035)	0.869 (0.088)	0.067 (0.076)
GB2	ALasso	40.474 (11.7331)	0.447 (0.0858)	0.710 (0.0605)
	SCAD	20.966 (7.6121)	0.517 (0.0614)	0.315 (0.1896)
	SIS + SCAD	10.314 (1.0797)	0.403 (0.0427)	0.042 (0.0712)
	SLasso	13.65 (2.038)	0.499 (0.052)	0.077 (0.081)
	Edc + EBIC	14.006 (1.657)	0.67 (0.014)	0.0273 (0.0785)
GB3	ALasso	22.464 (2.4414)	1.000 (0.000)	0.5495 (0.0498)
	SCAD	11.000 (0.000)	1.000 (0.000)	0.091 (0.000)
	SIS + SCAD	9.964 (0.6897)	0.992 (0.0764)	0.107 (0.050)
	SLasso	10.182 (0.475)	0.667 (0.006)	0.015 (0.039)
	Edc + EBIC	10.158 (0.440)	1.000 (0.000)	0.0139 (0.038)

, we report the simulation results under the conditions for GB1, GB2, and GB3. We observed that SLasso and Edc + EBIC performed better. SLasso had the highest PDR while Edc + EBIC had the lowest FDR. Thus, when there was some correlation among the features, Edc + EBIC still performed well.

3.2.3. Real Data Example

The new method was applied to the gene expression data used in [1]. For the data and details of data collection and variable definitions, see [1].

This study aimed to find the probes among the remaining 18,975 probes most closely related to TRIM32. The response variable was the expression level of probe 1389163_at. The features were the expression levels of the remaining 18,975 probes. Of the 18,975 probes, the top 3000 probes with the largest variances were considered. The expression levels were standardized to have mean 0 and standard deviation 1.

In our analysis of the data, for each of the 100 replications we selected a random sample of size 100 from 120 rats and applied Edc + EBIC to the sample. From these 100 replications, Edc + EBIC selected the distinct probes $1367705\_at$ and $1367728\_at$. In [1], the probes selected by five (5) variable selection methods are reported. The probes selected by Edc + EBIC did not intersect with any of the probes these methods selected. Among the probes selected by these five (5) methods, some intersected, but this is not a surprise because these methods essentially maximized only the linear relationship between TRIM32 and each of the probes. Since Edc + EBIC is capable of detecting and maximizing both the linear and nonlinear relationships that might exist between TRIM32 and each of the probes, as evidenced in the simulation studies by the high PDR and low FDR, we are convinced that the two probes selected by our method are the most associated with TRIM32.

4. Conclusions and Discussion

From the simulation results in Table 6 and Table 7 we observed that, as the sample size (n) increased, Edc + EBIC selected on average the expected number of predictors and did so with decreasing standard deviations, meaning that through the simulation runs more and more of the selected predictors were close to the expected number of relevant predictors. We also observed the positive discovery rate of $100\%$, indicating that on average for each simulation run, out of the selected features all of the relevant features were selected. Of greater importance was the small false discovery rates recorded as the sample size increased.

We found through simulation studies that when we applied the Energy Distance Correlation Sure Independence Screening proposed by [9] for variable screening followed by a regularization method such as SCAD and ALasso, the average model size selected was higher than expected and with high standard deviations.