Importance of feature selection stability in the classifier evaluation on high-dimensional genetic data

Machine learning models with embedded feature selection

We focus only on these models, due to the following reasons:

• using any classifier on data such as gene expression, where the number of objects is much lower than the number of features will likely conduct feature selection anyway;

• we have been working on a classifier that minimizes convex and piecewise linear criterion function (CPL) with L₁ regularization which we would like to test;

• computational efficiency is a significant constraint, as employing wrapper methods such as SVM-RFE (Guyon et al., 2002; Guyon & Elisseeff, 2003) would render our extensive computations, which spanned weeks, unfeasible within a practical timeframe;

• and finally using other feature selection methods such as wrapper or filter would require the use of a classifier to measure accuracy anyway.

We can distinguish two types of machine learning models with embedded feature selection—one is based on L₁ regularization, while the other one is based on decision trees. One of the widely used anti-overfitting techniques is the regularization of the model parameters values by L₁ or L₂ norm. It can be used by many classifiers such as logistic regression, support vector machines, or neural networks. When both norms reduce overfitting, only L₁ norm tends to set weights for the features to zero, which works like feature selection. Let us introduce L₁ norm with our CPL classifier.

We assume that data is represented as m feature vectors x_j = [x_j1, …, x_jn]^T $\left(j=1,\dots ,m\right)$ of the same dimensionality n. The components x_ji of the vectors x_j are called features. Each feature vector x_j is labelled with y_j (y_j ∈ { − 1, 1}) its membership to the decision class. For example y_j =  − 1 describe patients suffering from a certain disease and the y_j = 1 describe patients without the disease. Then we define the convex and piecewise linear criterion function as: (1)${\displaystyle {\Phi }_C\left(\mathbf{w},\theta \right)=\sum _{j=1}^{m}max\left(0;1-y_j\left({\mathbf{w}}^T{\mathbf{x}}_j+\theta \right)\right)+\left(1-C\right)\sum _{i=1}^n\left|w_i\right|,}$ where w = [w₁, …, w_n]^T ∈ Rⁿ is the weight vector, and θ ∈ R¹ is the threshold of the hyperplane H(w, θ) and 0 ≥ C ≥ 1 is a hyperparameter that controls the strength of the regularization penalty. A lower value of C increases the strength of the regularization and leads to more weights being pushed towards zero, resulting in a sparser model with fewer features. Conversely, a higher value of C reduces the impact of regularization, allowing the model to fit the training data more closely.

The basis exchange algorithm allows us to find the minimum efficiently, even in the case of large multidimensional datasets (Bobrowski, 1991). In our experiments we used our own python implementation.

Similar to our criterion function Φ_C (Eq. (1)), L₁ regularization can be introduced to LR and SVM algorithms. In our experiments we used the scikit-learn (Pedregosa et al., 2011) implementation.

Another embedded feature selection method is one based on decision trees. In the process of building the tree on every internal node the feature needs to be selected for test that is based on gini or entropy criterion (Tangirala, 2020). For datasets with much more features than objects n >  > m only a small number of features will be used even for a fully grown tree without pruning. In practice many ensemble techniques are used with decision trees like bagging or boosting. We decided to use the RF algorithm in our experiments also from the scikit-learn package (Pedregosa et al., 2011).

Feature selection stability measures

The issue of feature selection method sensitivity to minor changes in training data has been explored by several authors. If the selection of features changes significantly when a different sample from the same training data is used, the method is deemed unstable. On the other hand, if the feature subset remains largely unchanged despite alterations in the data, the method is considered stable. Although this notion is straightforward, there is currently no universally accepted metric for quantifying stability, and various suggestions have been put forth in the literature.

One of the first measures used to assess stability were those based on the intersection of feature sets, such as the Tanimoto distance (Kalousis, Prados & Hilario, 2007) and the Jaccard Index (Saeys, Abeel & Van de Peer, 2008). These measures allow for the evaluation of the similarity between feature sets of different cardinalities but do not account for the total number of features, and therefore, do not provide a correction for chance. However, in our case, where the total number of features is much higher than the number of selected features, this correction is minimal. (2)${\displaystyle S_J\left(F^{c_i},F^{c_j}\right)=\frac{|F^{c_i}\cap F^{c_j}|}{|F^{c_i}\cup F^{c_j}|}=\frac{r}{c_i+c_j-r}}$ where F^c_i and F^c_j are two subsets of features with cardinalities c_i and c_j, respectively, and r is the number of shared features. This measure ranges from 0 to 1, taking the value of 0 when there are no common features and the value of 1 when both subsets are identical.

In our previous research, we used a measure proposed by Lustgarten, Gopalakrishnan & Visweswaran (2009) as an extension of the Kuncheva similarity measure (Kuncheva, 2007). The Kuncheva metric introduces a correction for chance but is applicable only to two datasets with the same cardinality. Lustgarten’s extension allows for its use with subsets F^c_i and F^c_j of different cardinalities c_i and c_j, which is the usual case in embedded feature selection methods. It is calculated as follows: (3)${\displaystyle S_L\left(F^{c_i},F^{c_j}\right)=\frac{r-\frac{c_i{c}_j}{n}}{min\left(c_i,c_j\right)-max\left(0,c_i+c_j-n\right)},}$ where n is the total number of features. This measure ranges from (−1, 1] and provides a correction for chance. Negative values correspond to situations where common features appear less frequently than the expected overlap under random selection. Feature stability is then defined as the average similarity of all pairs of feature sets: (4)${\displaystyle ASM=\frac{2}{k\left(k-1\right)}\sum _{i=1}^{k-1}\sum _{j=i+1}^k{S}_{J/L}\left(F^{c_i},F^{c_j}\right),}$ where k is the number of obtained features sets.

Lately Nogueira, Sechidis & Brown (2018) proposed a new feature selection stability measure based on the frequency of selected features. (5)${\displaystyle Z=\left[\begin{array}{ccccc}\hfill z_{11}\hfill & \hfill z_{12}\hfill & \hfill z_{13}\hfill & \hfill \dots \hfill & \hfill z_{1n}\hfill \\ \hfill z_{21}\hfill & \hfill z_{22}\hfill & \hfill z_{23}\hfill & \hfill \dots \hfill & \hfill z_{2n}\hfill \\ \hfill ....\hfill \\ \hfill z_{k1}\hfill & \hfill z_{k2}\hfill & \hfill z_{k3}\hfill & \hfill \dots \hfill & \hfill z_{kn}\hfill \end{array}\right]}$ Each row Z_k of Z matrix refers to one feature selection subset, columns represents each feature, then observed frequency of selection of f feature is $p_f=\frac{1}{k}{\sum }_{i=1}^k{z}_{if}$. While it is still a generalization of Kuncheva measure in the sense that for subsets with the same cardinality it has the same values, it is computationally effective and holds the so called Maximum Stability (Deterministic Selection) property. The Nogueira features stability measure ϕ(Z) reaches its maximum, if-and-only-if all feature sets in Z are identical. Nogueira, Sechidis & Brown (2018) show that their measure holds this property while the Lustgraten measures violate the backward implication of it. (6)${\displaystyle \varphi \left(Z\right)=1-\frac{\frac{1}{n}\sum _{f=1}^n{s}_f^2}{\frac{\bar{n}}{n}\left(1-\frac{\bar{n}}{n}\right)}}$ where $s_f^2=\frac{k}{k-1}{p}_f\left(1-p_f\right)$ is the unbiased sample variance of the selection of the f feature and $\bar{n}=\frac{1}{k}{\sum }_{i=1}^k{\sum }_{f=1}^n{z}_{if}$.

It it worth noticing that both measures hold the other four desired properties mentioned by Nogueira, Sechidis & Brown (2018):

1. Fully defined. The stability estimator ϕ(Z) (Eq. (6)) should be defined for any collection Z of feature sets, thus allowing for the total number of features selected to vary.

2. Strict monotonicity. The stability estimator ϕ(Z) (Eq. (6)) should be a strictly decreasing function of the sample variances $s_f^2$ of the variables Z_f (column in the matrix Z).

3. Bounds. The stability ϕ(Z) (Eq. (6)) should be upper/lower bounded by constants not dependent on the overall number of features or the number of features selected.

4. Correction for chance. Subsets of features can share some of the features even if selected randomly.

trains-p-diff procedure for generating dataset splits

Cross-validation is a technique used in machine learning to evaluate the performance of a model by splitting the dataset into training and testing sets. The basic idea is to train the model on a portion of the data and then use the remaining data to evaluate the model’s performance. Typically, we divide a set of objects into a fixed number of k roughly equal portions, and then create a training set from k − 1 portions and, using the remaining one portion, assess the quality of the model. The procedure is repeated k times, and each time a different portion of the data acts as the test set.

In our research, the aim is to evaluate the discriminant model in terms of the stability of the features it uses and the impact of the composition of the training set on the features selected by the model. With this aim in mind, it is particularly important to construct the training sets correctly. In each cross-validation step, the number of objects in the training set should be the same and each training set should differ from any other training set used in the other cross-validation steps by a precisely specified number p (p ∈ {1, 2, …}) objects. In order to meet these conditions, we have developed our own method for generating training and test sets for cross-validation.

The parameters of the trains-p-diff procedure are the size of the training set t in each split generated in successive cross-validation steps and the number of objects p by which the training sets of any two splits differ. The dataset D on which the splits are to be made should have a count of at least t + p (|D| ≥ t + p). The procedure will result in train/test splits with counts of t and |D| − t respectively. Each selected pair of training parts from any two splits will have exactly t − p common objects and p objects occurring in only one of the sets from the pair. An example showing how the procedure works is shown in Table 1.

Using the trains-p-diff procedure described above, we can obtain reproducible training and test sets that allow the results of different discriminant models to be compared. We can also accurately control the size and degree of differentiation in the composition of the training sets, allowing us to assess the stability of the model.

Datasets

We selected seven gene expression datasets named: Breast, Colorectal, Leukemia, Liver, Prostate, Renal and Throat. They were obtained from the Curated Microarray Database (CuMiDa) repository provided at https://sbcb.inf.ufrgs.br/cumida (Feltes et al., 2019). All datasets were characterized by a significant number of features, in this case describing the expression of specific genes, and the presence of two decision classes. Table 2 shows the most important parameters of the datasets.

The CuMiDa database is a collection of several dozen carefully selected microarray data sets for Homo sapiens. These data sets were chosen from over 30,000 microarray experiments available in the Gene Expression Omnibus (GEO) (Edgar, Domrachev & Lash, 2002), following strict filtering and evaluation criteria. In addition to technical aspects, proper data curation of GEO datasets also supports ethical standards. Thorough documentation of experimental methods, participant consent, and compliance with ethical guidelines help protect the rights and privacy of individuals contributing to scientific research.

Experimental setup

For every experiment we used the trains-p-diff procedure described earlier for p ∈ P (P = {1, 2, 3..10, 12, 14, …, 20}), and for every p there was usually 10 of datasets splits to train and test, sometimes less if the dataset was small and the p high. This was repeated three times for different random seeds. For every split a classifier was trained and metrics were obtained. Accuracy was calculated on a test set, number of selected features and feature selection stability metrics were calculated based on a fitted classifier model (see Fig. 1).

We performed two kinds of experiments - in the first one, called the single classifier, we ran one classifier on all datasets with different values for C. From this experiment we can see how different values of C influence the results. This experiment did not really allow us to compare different classifiers, which is why we performed the second experiment where for each dataset we ran all classifiers with the C parameter set in such a way, that all of them selected a similar number of features (see Fig. 2). In this design we can compare classifiers and for example say which of them is more stable in terms of selecting the same features.

As mentioned earlier, implementations of classifiers LR, SVM and RF offered by the scikit-learn (Pedregosa et al., 2011) library, as well as our own implementation of the CPL classifier, were used in the experiments. Most of the parameters of the classifiers were left at default values apart from those listed below:

• LogisticRegression(C, penalty =‘l1’, solver =‘liblinear’, max_iter =100000)

• LinearSVC(C, penalty =‘l1’, dual =False, max_iter =100000)

• RandomForestClassifier(n_estimators =100, random_state =0, min_samples_leaf =2)

• GenetClassifier(C)

In the case of the LR, SVM and CPL classifiers, the influence on feature selection, in particular on the number of selected features, was obtained by manipulating the values of the C parameter. In the RF classifier, on the other hand, features were selected on the basis of a threshold value set for the collection of feature_importances_ (the Gini importance) created on the basis of the trained model.