RENT -- Repeated Elastic Net Technique for Feature Selection

Feature selection is an essential step in data science pipelines to reduce the complexity associated with large datasets. While much research on this topic focuses on optimizing predictive performance, few studies investigate stability in the context of the feature selection process. In this study, we present the Repeated Elastic Net Technique (RENT) for Feature Selection. RENT uses an ensemble of generalized linear models with elastic net regularization, each trained on distinct subsets of the training data. The feature selection is based on three criteria evaluating the weight distributions of features across all elementary models. This fact leads to the selection of features with high stability that improve the robustness of the final model. Furthermore, unlike established feature selectors, RENT provides valuable information for model interpretation concerning the identification of objects in the data that are difficult to predict during training. In our experiments, we benchmark RENT against six established feature selectors on eight multivariate datasets for binary classification and regression. In the experimental comparison, RENT shows a well-balanced trade-off between predictive performance and stability. Finally, we underline the additional interpretational value of RENT with an exploratory post-hoc analysis of a healthcare dataset.


Introduction
A predictive task involves a dataset consisting of N -dimensional row vectors X = (x T 1 , . . . , x T I ) ∈ R I×N and an associated vector of target values y = (y 1 , . . . , y I ) ∈ T I , where the target space T may represent a set of classes (classification task) or a subset of the real numbers (regression task). In this study, our focus lies on generalized linear models (GLMs), which model the target as a linear combination of the inputs with weights β ∈ R N , followed by a transformation. The columns of the data matrix describe object characteristics, denoted as features. Since data arXiv:2009.12780v3 [cs.
LG] 22 Nov 2021 acquisition techniques evolve steadily, situations where the number of features N exceeds the number of objects I often occur. In such setups, mathematical obstacles, like spurious correlations and multicollinearity issues causing model overfitting, trigger the necessity to reduce the number of features by using some feature selection approach [1]. These issues are characteristic of various domains, including healthcare [2,3], biomedicine [4], text mining [5] and botany [6]. A successful feature selection approach will decrease the model complexity, improve the model stability and provide more useful model interpretations.
A feature selector θ F decomposes the data space into a direct sum of selected features (V 1 ) and non-selected features (V 2 ) according to the given feature set F ⊂ {1, . . . , N }, and projects all objects from R N to the subspace V 1 , i.e.
The goal of good feature selection is to determine the feature set F , which enables a predictive model to obtain the most accurate prediction. Predictive quality is measured using a metric q (ŷ, y), such as F1 score, whereŷ, y ∈ T Itest denote the vectors containing predicted target valuesŷ after feature selection and ground truth target values y, both referring to a set of test data X test of size |X test | = I test . An optimal feature set F is characterized by F = arg max F ⊂{1,...,N } q(ŷ F , y).
A taxonomy of feature selection techniques distinguishes between filter, wrapper, and embedded approaches. Filter approaches rank features by an importance criterion, such as mutual information or correlation coefficients between features and target variables. Baseline filters include the Fisher score [7] and the Laplacian score [8] as well as algorithms from the relief family [9]. Approaches like mRMR [10] or the stratified feature weight method [11] aim to resolve the issue that correlated and redundant features are not well handled by classical filters [12]. A combination of different filter approaches is suggested in [13]. Wrapper approaches select features concerning their prediction performance. By training supervised models on different subsets of the entire feature set, the subset delivering the most accurate predictions on a test set is chosen. This strategy often causes overfitting issues and high computational costs [14]. Prominent wrapper approaches are forward/backward selection, such as recursive feature selection [15], and heuristic searches like simulated annealing or genetic algorithms [16].
The third category of feature selection methods, embedded feature selection, integrates the selection step directly into the learning algorithm. A class of embedded methods, which is particularly important in this work, comprises regularization for GLMs: During parameter estimation, regularization terms are added as penalties to the target function. While the well-established LASSO [17] uses an L1 term λ 1 (β) = |β| for this purpose and delivers a sparse parameter vector, L2-regularization λ 2 (β) = β 2 handles multicollinearities by pulling the L2-norm of the parameter vector β towards zero. The effects of both regularization terms are combined in the elastic net λ enet (β) [18], defined as with parameters α ∈ [0, 1] and γ to weight the regularization terms and to define the regularization strength, respectively. Other representatives of embedded feature selection models are tree-based models, such as decision trees or regression trees. Ensembles of tree-based architectures are called random forests [19]. Graph-based approaches together with elastic net regularization further play a key role in recent works [20][21][22], where the authors demonstrate a graph-based structurally interacting elastic net method incorporating pairwise relationships between objects via a feature graph, or [23], proposing a solution for 2,0 -norm regularized feature selection via linear discriminant analysis.
Most feature selection approaches suffer from the phenomenon that minor changes in the random initialization or train-test-split of the model lead to major variations in the selected feature set-this issue is referred to as lack of stability and is investigated in [24] and [25]. In agreement with [26] and [27], the authors argue that L1 regularisation on GLMs is generally unstable. They claim that the issue can be resolved by investigating ensemble feature selection, where θ F is derived from a set of independently trained (elementary) feature selectors θ F1 , . . . , θ F K , such that θ F = φ (θ F1 , . . . , θ F K ). The operator φ acts as a meta-model based on information from the elementary models θ F k , k = 1, . . . , K. A basic approach is to build such a meta-model by counting the frequency of selection for each feature across all feature sets F k , expressed by where t 1 ∈ [0, 1] is a scalar representing a minimum selection frequency threshold. This approach assumes that each elementary feature set F k consists of a subset of important features (with correspondingly higher probabilities for being selected), and a small, random subset of unimportant features. The final, selected feature set F is less likely to contain unimportant features than each of the elementary feature sets F k . Hence, model stability is increased as shown by Meinshausen and Bühlmann [28], who propose such a feature selector named stability selection. Even though the stability selection framework is intuitive and reasonable, the corresponding feature weights may be small-not significantly different from zero-or have alternating signs across the elementary models. Thus, features might be selected although resulting in ambiguous or contradictory information and hence, deteriorating interpretability and predictive performance. This means that further insights into the predictive power of features have to be gained from the distribution of weights, which is not considered by Meinshausen and Bühlmann [28].
The present work suggests the novel repeated elastic net technique (RENT) for feature selection. RENT is based on the idea of model ensembles discussed in [28]. Besides merely calculating the frequency of each feature, we also focus on the empirical distribution of the feature weights resulting from elastic net regularized models. Thereby, we extend the model ensemble framework to combine three rigid selection criteria: 1) how often is a feature selected?; 2) to which degree do the feature weights alternate between positive and negative values?; 3) are feature weights significantly different from 0? The final feature selection of RENT consists of the features that satisfy all three selection criteria. When required, the RENT framework can be extended with additional custom criteria to refine the feature selection process according to the user's a priori insights and requirements. By taking elastic net regularization into account, RENT aims at optimizing predictive performance and model stability simultaneously. In contrast, the concept of stability selection focuses on model stability as the primary target. We suggest a hyperparameter selection procedure based on the Bayesian information criterion (BIC) to balance the number of features and the predictive performance.
In the experiments section, we explore and evaluate RENT extensively using real-world datasets for both classification and regression problems. In addition, we use the information provided by the ensemble models for an exploratory post-hoc analysis with statistical tools, including principal component analysis. Our implementation (in Python code) is publicly available and published in the Journal of Open Source Software [29].

Repeated Elastic Net Technique for Feature Selection
In this section, we present the methodological concept of RENT which relies on regularized logistic regression for binary classification problems and regularized linear regression for regression problems. We introduce the idea of elastic net regularization combined with repeated training of machine learning models on unique subsets of the training data to investigate feature selection stability. Finally, we define three quality metrics that influence the feature selection.

Ensemble Training and Selection Criteria
Given a set of training data X train = {x i : i = 1, . . . , I train } where x i denotes an object from the N -dimensional feature space, our concept builds on sampling K unique i.i.d. (independent and identically distributed) subsets X (k) train ⊂ X train of size I   Table 3.
The evaluation of each model M k is performed on the validation set X train (here, \ denotes the set difference operator). To further improve robustness, we include the option to introduce more variation across the K models, by randomly varying the number of objects drawn from X (k) train between the models within user-specified limits. For each feature n in X train , n = 1, . . . N , we observe the trained weights β k,n throughout models M k , k = 1, . . . , K. For the purpose of feature selection, we acquire relevant information about the importance of feature n across all models from β n = (β 1,n , . . . , β K,n ). All such vectors β n , n = 1, . . . , N , are aggregated in a matrix B of dimension (K ×N ). Since all models comprise L1 regularization terms, the vectors of feature weights β n are typically sparse. However, entries are not constant due to 1) variations in the training subsets and 2) numerical deviations in the parameter optimization. Hence, a straightforward measure of feature relevance is the relative frequency c(β n ), counting how often a feature was selected on average across the K models or, in other words, calculating the relative frequency as an estimate of the probability for the parameter of the n-th feature to be non-zero: Furthermore, we observe two other empirical summary statistics of the feature parameter estimate distributions in the rows of B: the feature-specific mean µ(β n ) and variance σ 2 (β n ) of the feature weights In general, we consider the n-th feature to be a candidate for selection in RENT if 1. c(β n ) is large, i.e. the feature is selected in many of the K elastic net models; 2. the estimates in β n resulting from the K models do not alternate much between positive and negative signs (stability); 3. the mean of distribution resulting from the K parameter estimates in β n is significantly non-zero.
These three simple and transparent requirements may be formulated in corresponding mathematical expressions, to form three quality metrics for assessing a feature n: where t K−1 (.) denotes the cumulative distribution function of Student's t-distribution with K − 1 degrees of freedom.
Considering the second quality metric τ 2 (β n ), the ideal case for feature n would be that all weights have the same sign-either all positive or all negative. In case of constant signs among all weights, τ 2 (β n ) equals τ 1 (β n ). Though, for a considerably large K, we should expect that at least slight sign variations for some features may occur. τ 2 (β n ) simply allows the user to define a required minimum proportion of the parameter estimates to have the same sign. The third quality metric τ 3 (β n )-identifying consistently high model parameter estimates-is chosen such that it corresponds to the well-known statistical Student's t-test with rejection of the null hypothesis H 0 : µ(β n ) = 0. In case that the null hypothesis holds, the test statistic will follow a Student's t-distribution with K − 1 degrees of freedom. The deployed term evaluates the probability of the test statistic under the H 0 -distribution and thus, provides a thresholding at the chosen level of significance.
In order to define feature selection criteria from quality metrics τ 1 (β n ), τ 2 (β n ) and τ 3 (β n ), we introduce corresponding cutoff values t 1 , t 2 , t 3 ∈ [0, 1]. Specifically, a feature n ∈ F is added to the selected feature set F , if it satisfies all three criteria: τ i ≥ t i , ∀i ∈ {1, 2, 3}. Further criteria can be included by the user if necessary. In the provided setup, these quality metrics may be considered as hyper-parameters of the RENT method, allowing the user to regulate the feature selector, by tuning the thresholds t 1 , t 2 and t 3 . The cardinality of the selected features F will increase, if any of these thresholds are decreased and vice versa. All three metrics, τ 1 , τ 2 and τ 3 , are bounded by the interval [0, 1], which facilitates the specification of appropriate thresholds. Since τ 3 can be associated with a Student's t-test, the threshold t 3 for a 5% or 1% significance level, corresponds to the thresholds t 3 = 0.95 and t 3 = 0.99, respectively.

Hyperparameter selection
RENT involves hyperparameters at different stages of the method: before training the elementary models, regularization parameters γ and α control the restrictiveness of the feature selection in the ensemble, followed by the parameters t 1 , t 2 , t 3 determining the final feature set. Thereby, the latter cutoff parameters are (a) dependent on the choice of the regularization parameters, and (b) have mutual dependencies.
Hyperparameter selection is commonly performed using an additional validation dataset or cross-validation-both options are not optimal for RENT, since a validation subset would reduce the number of objects in a high-dimensional dataset even further, and cross-validation would add a substantial computational burden to the procedure. Thus, we deploy an alternative approach from statistical model selection: the Bayesian information criterion (BIC) delivers a trade-off between the information content (quantified as the likelihood) of the model and the model complexity in terms of the number of estimated parameters [30]. BIC is defined as whereL denotes the estimated likelihood of the predictive model, and ρ denotes the number of estimated model parameters. In contrast to similar information measures like the Akaike information criterion (AIC), BIC is known for stronger penalization of model complexity leading to a lower number of selected features, which is favorable in the case of RENT. By minimizing BIC, models with high information content and low complexity are favored. In ordinary linear regression models and other standard GLMs, the number of estimated model parameters equals the number of variables, i.e., features, plus one parameter for the offset β 0 ; thus, we set ρ = |F | + 1. The likelihoodL can be determined from the distribution assumptions of the GLM model, such as the normal distribution of errors in the ordinary least squares regression model, resulting in the sum of squared errors (SSE) as negative log-likelihood function.
RENT uses a two-step hyperparameter estimation procedure: a grid search for regularization parameters α and λ with BIC as target function is performed on the full training dataset first (step 1). Then, the RENT ensemble is trained given the best regularization parameter combination. Finally, in step 2 another grid search for cutoff parameters t 1 , t 2 , t 3 is performed using the same concept as in step 1.

Training runtime complexity of RENT
Since RENT is an ensemble method built on GLMs as elementary models, the runtime complexity of RENT is expressed as a multiple of the runtime complexity of GLMs, denoted by O GLM . In essence, O GLM depends on the applied type of GLM, the parameter optimization algorithm, and the implementation. For instance, a runtime complexity of O GLM = O(N 3 + I train · N 2 ) is reported for Lasso by reducing the computation to solving a least squares regression problem [31]. Variants using iterative algorithms are rather judged by the overall experimental runtime and the runtime complexity per update cycle, while the number of iterations is hard to determine a priori-such information is provided for GLMs with elastic net regularization in [32].
Given the first variant, RENT runs an ensemble comprising K independent GLMs, each trained on a number of N features, which delivers a complexity of train < I train denotes the sample size of each subset during RENT training. In addition, hyper-parameter tuning requires training c GLMs, where c is a constant given by the number of level combinations for regularization and cutoff parameters, resulting in O cN 2 · (N + I train ) .
In total, an upper bound to the full runtime complexity of RENT is given by

Experiments
We demonstrate the potential of RENT as a feature selection method through experiments on multiple datasets. First, we verify the overall concept in a validation study in Section 3.4. Second, we evaluate the performance of RENT in comparison with seven feature selection methods and a baseline elastic net regularized model, in Section 3.5. Based on one dataset, we illustrate how the stability of RENT behaves compared to the stability of established ensemble methods based on the number of unique elementary models K ∈ N.

Experimental Setup and Datasets
Experiments are conducted on multivariate datasets from various domains, including real-world data and synthetic data for both binary classification (class.) and regression tasks (reg.); datasets are listed in Table 1. The size of each dataset is denoted via the number of features (#f eat) and the number of objects (#obj) divided into train and test sets (train/test). Train-test-splits are performed by stratified random sampling. Further, the class balance indicates the percentage of class representation for each classification dataset (train/test).
The broad selection of use cases, including high-dimensional datasets, demonstrates the flexibility and applicability of RENT. Simulated datasets (c0) and (r0) were produced using scikit-learn [33] functions make classification and make regression, respectively. For the MNIST dataset, two binary classification problems are defined by restricting the classes: MNIST cl1,cl2 indicates that only instances from classes cl 1 and cl 2 , where cl 1 , cl 2 ∈ {0, . . . , 9}, were used, ignoring objects from other classes.
A feature selector is trained on X train , then the training data X train is projected into the subspace spanned by the selected features. This column-reduced training dataset is denoted by X train . In our experiments we train an unregularized linear/logistic regression model M on X train . Evaluation is based on the predictive performance obtained from M on the previously unseen test data X test . It is important to note, however, that it may be necessary to use regularization for the model M to avoid overfitting, especially if the reduced X train has more features than objects.

Evaluation Metrics
We use two different measures for quantitative evaluation of the prediction performance in classification settings: F1 score (F1) and Matthews correlation coefficient (MCC) [39]. The F1 score represents the harmonic mean of precision (PR) and recall (RC). Denoting the entries of the confusion matrix by TP (true positive), FP (false positive), FN (false negative), and TN (true negative), the performance measures are defined as follows: where P = T P + F P and N = T N + F N are the sums of the predicted positives and negatives, respectively. Note that F1 scores can be calculated for both class labels, depending on which class is considered as "positive". F1 is more appropriate than accuracy for imbalanced class distributions because the larger class dominates the latter. A disadvantage of F1 is that it does not take into account TN. Therefore, MCC provides more representative results if both classes are equally relevant in the prediction problem and the number of TN objects is high. F1 score, precision, and recall are bounded between [0, 1], where 0 represents a complete disagreement between predicted and actual class, and 1 denotes a perfect match. MCC is bounded between [−1, 1], where −1 denotes that all objects are classified incorrectly, 0 indicates complete randomness, and 1 denotes correct classification of each object, respectively.
For regression problems, we evaluate the root mean squared error of prediction (RMSEP) on the test dataset X test with cardinality I test , defined as and the coefficient of determination (R 2 ) [40] where y i represents the true output of object x i ,ŷ i represents the prediction of y i andȳ represents the mean of the outputs y i , i ∈ {1, . . . , I test }. While RM SEP is always non-negative we seek its minimization. R 2 on the other hand may take negative values but has an upper bound of 1 (associated with perfect predictions) and we therefore seek its maximization.
Besides predictive performance, selection stability is assessed using a measure suggested in [24] evaluating the different outcomes of multiple feature selection runs in a combinatorial way. Specifically, the suggested measure computes a ratio between the sample variance of observed feature frequencies and the theoretical variance, given that the feature selector is stable (null hypothesis). The authors clarify that their measure fulfills five consistency criteria and is asymptotically bounded by the interval [0, 1], where 1 denotes optimal stability. Their concept of measuring feature selection stability by aggregating multiple independently trained models underlines the relevance of our ensemble approach and supports the idea to achieve stability by combining K independent feature selection model runs.

RENT Hyperparameter Selection
In Section 2.2 we introduce hyperparameter selection for both the elastic net modeling and the three cutoff parameters based on the BIC. More precisely, we evaluate the elastic net hyperparameter combinations of γ ∈ {1e −2 , 1e −1 , 1} and α ∈ {0, 0.1, 0.25, 0.5, 0.75, 0.9, 1}. To find the best combination concerning BIC, we train a single logistic/linear regression model with each pairwise combination of hyperparameters γ and α on the training dataset. After determining optimal elastic net parameters γ and α for a given dataset in Table 1, all ensemble models M 1 , . . . , M K in RENT are trained with these parameters. Once all K models are fitted on their respective training subsets X (k) train , we select the cutoff hyperparameters t 1 , t 2 and t 3 with BIC, in the same way as for the elastic net hyperparameter search. For this purpose we perform a grid search on t 1 ∈ [0.2, 1] with stepsize 0.05, t 2 within the same range and t 3 ∈ {0.9, 0.95, 0.975, 0.99}, representing different significance levels in the t-test. A comparison of three different hyperparameter settings based on the lowest, median, and highest BIC values is shown for dataset c0 in Fig. 2 for a varying number of elementary models K ∈ {5, 10, 50, 100, 300, 500}. For each hyperparameter setting (t 1 , t 2 , t 3 ) leading to the lowest, median and highest BIC, respectively, 30 independent runs of RENT are carried out. Each run is conducted on the same training dataset but with a distinct (random) model weight initialization. Performance is measured via the MCC; runtimes are given in seconds and refer to one single run for each method. Across all 30 independent runs, the mean is calculated together with the empirical 2.5% and 97.5% quantiles (corresponding to a two-sided 5% confidence interval) for stability, performance, and runtime, respectively.
We can observe that, as expected, the setup of RENT with optimal hyperparameters (lowest BIC) outperforms those settings with median and highest BIC and achieves the highest stability. Especially RENT based on hyperparameter settings (t 1 , t 2 , t 3 ) with the highest BIC is unstable, even though the performance remains in an acceptable range. Regarding runtime, it takes about 600 seconds for RENT with median BIC to run a single model for K = 500, which is much longer than for the other two settings. A reason for this might be a harder optimization task for specific hyperparameter combinations where it takes more steps for the logistic regression model to converge.
In summary, we observe the following behavior of RENT with lowest BIC at an increasing number of models K: • on average, good MCC and stability are achieved simultaneously, even with low K; • as expected, stability increases significantly from 0.75 for K = 5, saturating at a value close to 1; • average MCC shows little change from K = 100 to K = 500; • runtime increases linearly with the number of trained models.
Hence, our results support the analysis in [4] that repeated use of regularized elastic net models is useful to achieve stable and reproducible results, while keeping the predictive performance at a high level. Our observations further indicate, that no major benefit can be achieved by increasing the number of models to more than approximately 100 with respect to the observed metrics on the given dataset. Therefore, K = 100 seems to be a valid default regarding the trade-off between stability and time for the datasets used in this study. If computation time is critical, the user may set K to a lower number but needs to consider that the distribution of the weights may be insufficiently covered and that this may have an impact on the stability of feature selection.
Alternatively, instead of using BIC, the user may set hyperparameters γ and α manually or use cross-validation to obtain a customized trade-off between predictive performance, stability, and the number of selected features. Note that this approach may be more subjective and that the computational cost can be higher than using BIC, especially if cross-validation is used.  For regression datasets, the analog procedure is applied using R 2 as a quality metric. Both tests are conducted at a significance level of 0.05.

Validation Study of Features Selected with RENT
(VS1) Compare a number of ∈ N randomly selected feature sets, representing inefficient feature selections, to the features selected by RENT. The steps of the procedure are: (a) sample independent, random feature subsets from X train , containing ∆ features each, where ∆ corresponds to the number of features selected by the RENT approach (b) train a new model for each of the feature sets by restricting X train to those features (c) predict the labels of X test with each of the models and compute MCCs (d) perform a Student's t-test, assuming as null hypothesis that the MCC value obtained from RENT is drawn from the same distribution (VS2) Compare the predictive performance of a model based on features selected with RENT on the real X test labels, to the predictive performance of randomly permuted labels of X test . The steps of the procedure are: (a) train a model on X train with the features selected with RENT (b) randomly permute y test −times and compute the average MCC over the permutations (c) perform a Student's t-test, assuming as null hypothesis that the MCC value obtained from RENT is drawn from the same distribution Performance results from (VS1) and (VS2) provide a reliable indicator of whether models based on features selected by RENT perform better than models based on randomness. Table 2 shows the average MCC of (VS1) and (VS2) in the columns M CC/R 2 . All corresponding p-values from the Student's t-tests are significantly lower than 0.05, mostly below 1e −15 , where equals 100. Since the standard deviation of the mean decreases with the sample size, a higher explanatory power of the Student's t-tests can be achieved by setting to a larger value. However, the runtime increases linearly in . The estimated densities of (VS1) and (VS2) for the Breast cancer Wisconsin dataset (c3) are plotted in Fig. 3. In general, these two validation studies are not limited to RENT and may be applied to other feature selection methods and metrics, as well. Overall, the null hypotheses in both VS1 and VS2 were rejected for all datasets, indicating that RENT performs significantly better than models based on randomness as described in the validation approaches. The experimental results in Table 2 show that (VS2) is close to zero across all datasets, as one would expect from the experimental setup. (VS1) performance is similar to the performance of RENT models for datasets c1. This fact may be explained by the individual information content of each feature: especially two-class subsets extracted from MNIST contain many mutually or highly correlated features. Therefore, many different feature combinations lead to good predictions.
RENT results for MNIST (datasets c1 and c2) are visualized in Fig. 4. We observe that 1) different features are relevant for distinguishing the class pairs 0-1 and 4-9 and 2). features relevant for 0-1 are typically located in the center of the image, whereas those relevant for 4-9 are more distributed across the image. Overall, distinguishing between 4 and 9 is more complex. Therefore, the number of selected features is much higher in this case than when classifying the numbers 0 and 1.

Comparison of RENT with Established Feature Selectors
The validation study in Section 3.4 showed that RENT is a valid feature selection approach for all datasets used in this study. Hence, we compare RENT to the methods listed in Table 3 as follows: 1) seven established feature selectors applied to classification datasets; 2) five feature selectors applied to regression datasets; 3) a baseline logistic/linear regression model M • with elastic net regularization [17]. For each feature selector, software implementations are publicly available. To compare RENT to traditional filter methods, we consider the Laplacian score (L-score) [8], Fisher score (F-score) [7], mRMR [10], and a representative of the relief family, reliefF [41]. Specifically, we select the top features according to the scores provided by each filter method. Further, we study the behavior of recursive feature elimination (RFE) [15] representing a wrapper based approach. Finally, our comparison also involves two prototypes of state-of-the-art ensemble feature selectors: stability selection (StabSel) [28] and the random forest [42], which can be used for both classification (RFC) and regression (RFR) problems. In contrast to other methods, RENT and M • share the advantage that the user does not have to specify the size of the selected feature set as input. Instead, the number of selected features is indirectly controlled via the elastic net regularization parameters γ and α. Similarly, for StabSel the exact number may be specified optionally. Otherwise, it is determined indirectly by a cutoff and an upper bound per-family error rate (PFER) [43]. For fair performance comparison of the remaining investigated methods, the size of the selected feature set is set to the number of features returned by RENT, denoted as ∆.
For StabSel, we perform a 5-fold cross-validated grid search to estimate adequate parameter settings. The elementary feature selection method is the logistic regression model with L1 regularization; the number of models equals K. Furthermore, we perform a grid search for stability selection on the interval [0.6, 0.9] for the cutoff value and [0.05, 0.95] for the PFER value, with a 0.05 step size each. In our study, the random forest serves as a filter, delivering a ranking of the features. The features with the highest ranks are selected and used as input for M where the number of the selected features corresponds to the number of features ∆ selected by RENT. To fit the random forest model, we set the number of unique trees to K. Other parameters are set to the defaults. Computations are performed on standardized train datasets. All model parameters used for the established methods, such as the neighborhood graph construction in L-score or the step size in RFE, are set to the default values, except for c4 and c5, where the step size is increased to 100 in order to obtain results in a moderate runtime. The regularization parameters γ and α for M • are set to those used for RENT. The results for all datasets and methods are provided in Table 4 for binary classification datasets and in Table 5 for regression datasets. 1  For classification problems, the results achieved with RENT feature selection are competitive with the best results of the other methods, yielding better or equally high F1 scores for predicting class 0 in five out of six datasets. For c0, the performance is only 0.01% below the top value of 0.75. Also, for class 1, RENT achieves the highest performance for four out of six datasets. For dataset c5 the performance is close to the best F1 score. MCC, which is more robust than the F1 score for unbalanced class settings, is highest for five out of six datasets with RENT feature selection. For c5, M • has a higher MCC, but with a much higher number of features (593 features) than RENT (16 features).
With the regression datasets, RENT achieves a performance superior to StabSel, RFR, L-score, mRMR, and RFE and competitive performance to M • for both measures, RMSEP and R 2 .
In many cases, M • is not able to restrict the number of features as efficiently as RENT. Using the same regularization parameters as for RENT, M • selects the following number of features: 290 (c1) vs. 37 for RENT and 593 (c5) vs. 16 for RENT, respectively (see Table 2).
Overall, StabSel achieves good results for all datasets underlining the merits of ensemble feature selection concepts. Note that no results could be obtained for dataset c3 since no feature reached a sufficient selection frequency across all models. The random forest yields competitive results for most datasets in classification (RFC) setups but performs noticeably worse in regression (RFR) setups. L-score and mRMR appear to provide weak performance scores compared to their competing feature selectors. For L-score, the low scores can be explained by its unsupervised setup, which makes it harder to relate the model to any target variable. With mRMR, especially the performances for c1, c2, and c4 are weak. For c1 and c4 this weakness can partly be attributed to the available implementation, which produced an error for these datasets (denoted with superscript b in Table 4). On the other hand, F-score performs well, especially for predicting class 0. The reliefF method achieves good results for c5 but is among the poorest feature selectors for c0, c1, c2 and c3. Neither F-score nor reliefF are applicable to regression problems using the available implementations. Dataset c4 is of particular interest since opposite behavior can be observed among the feature selectors. RENT, StabSel, RFC, and F-score perform well when predicting class 0, whereas the other methods achieve higher scores for class 1. Hence, we assume that the features selected from c4 introduce a bias towards class 0 or class 1, respectively. In terms of MCC, which accounts for both classes in parallel, the best results are achieved by RENT, StabSel, RFC, and F-score.   Fig. 2. While RFC achieves a similar performance as RENT, it is the most unstable ensemble approach for lower number of trees K. Even for high K, RFC never achieves the same stability as RENT and StabSel. Furthermore, RFC has the highest variance in MCC. On the other hand, StabSel reaches higher stability but lower MCC scores than its competitors between K = 50 and K = 100. The provided stability analysis underlines the strong properties of RENT, compared to random forests, which are known to be unstable in multiple scenarios [47].
Regarding the computational costs 2 of the ensemble feature selectors in our study, Fig. 5 demonstrates that the runtime increases linearly in K for all methods. RENT takes longer to compute which might be caused by the fact that the implementation does not yet exploit the full potential for runtime optimization and different implementations and programming languages were used for elementary operations.

Exploratory Post-Hoc Analysis
As an ensemble model approach, RENT offers additional information that can be integrated into exploratory post-hoc analyses. The two post-hoc analyses presented in this section give the user tools to 1) further investigate objects in the dataset and identify which of those are difficult to predict and which not; 2) exploit this information in a principal component analysis model trained on the selected features to understand why some objects are difficult to predict.

Analysis of Training Objects
Based on the ensemble of elementary models in RENT, it is possible to compute summary statistics on a singleobject level. Such information may contribute to improved interpretability of the model in general and single objects in the data in particular. For this purpose, we analyze the predictions of individual objects across all models M k , k = 1, . . . , K. For binary classification problems, we observe the distribution of correct and incorrect classifications of single objects in X k val , and thereby gain insights into the consistency of assigning an object to its true class. From a statistical perspective, this means that we can identify objects with deviating properties belonging to the same class based on the information whether the label of an object is difficult to predict or not. For regression problems, we similarly use the mean absolute errors. Below, we will exemplify the proposed post-hoc analysis for dataset c3.
Given an object x i ∈ X k val , the logistic regression model M k outputs a class probabilityŷ i of x i being assigned class 1 (ProbC1). Among the K models built within RENT, we obtain a ProbC1 value each time an object x i appears in X k val , k = 1, . . . , K. Aggregating this information by object, we can derive statistics and describe the distribution of the ProbC1s for each object x i ∈ X k val by a histogram, as shown in Fig. 6. These results are generated from dataset  c3 (Breast cancer Wisconsin), where we denote a single object in the dataset as a cancer patient. Incorrect predictions provide evidence for patients that are hard to classify or show different characteristics compared to patients from the same class that are easy to classify. We observe that patient 3 belongs to class 0 and that the predicted probabilities of patient 3 are consistently below 0.5, which is the standard decision boundary for logistic regression models. In other words, patient 3 is predicted correctly every time she is part of X k val . Patient 6 belongs to class 0; however, the predicted probabilities are consistently above 0.5, meaning that her class label is always mispredicted. For patients 78 and 102 we observe probabilities both above and below 0.5, indicating that the class predictions of these two patients are rather uncertain, however, to a different degree. Fig. 6 reflects the detailed information provided in Table 6. With a % incorrect of 54.2%, the class predictions for patient 102 are extremely unstable among the 24 models, where this patient was part of the validation set. This type of information on prediction stability may provide a good starting point for detailed studies on how a hard-to-classify object differs from objects that are consistently assigned to the correct class. Thus, it could be of high relevance, inter alia for medical experts, who may identify patients with deviating data characteristics. These difficulties may arise from dominating phenomena in the measured features or measurement errors.
In this way, ensemble based approaches such as RENT allow in-depth analysis of the distribution of class probabilities rather than restricting to single class predictions.

Principal Component Analysis (PCA) on Selected Features
By a PCA [48] of X train , we can obtain a better understanding of the properties of objects and their relation to the features selected by RENT. Note, that unlike in machine learning, where PCA is often only used for feature extraction, visualization (plotting) of PCA scores, PCA loadings, and PCA correlation loadings [49] may be efficient for the purpose of model interpretation. Fig. 7a and Fig. 7b show the scores 3 of the first two principal components (comp 1 and comp 2) applied to the Breast cancer Wisconsin dataset, but with hues based on different information acquired from the ensemble. Every data point in the scores plot represents one object in the data or specifically to this dataset, one patient. The first two principal components explain 92.1% of the total variance in the data that are contained in the selected features. The remaining 7.9% are explained by the remaining principal components. In particular, Fig.  7a shows how the objects are distributed in the sub-space spanned by components 1 and 2. The two classes are well separated, with the objects of class 1 (circles) on the left side of the plot and class 0 (triangles) on the right side. Using results from RENT, the information in this plot may be further enhanced by coloring each object by its true class (class 0 -green triangles; class 1 -red circles) graded according to % incorrect in Table 6. Higher color saturation refers to a higher percentage of incorrect label predictions-suggesting that the object shows an anomalous behavior, which the model cannot sufficiently cover. In this example, objects with ambiguous classes accumulate in the middle area, most of them are close to the intuitive decision boundary. Fig. 7b shows the same scores as seen in Fig. 7a. However, the objects are colored by their average ProbC1 (the average probability of an object belonging to class 1). Again, objects with either a very high or a low value cluster on the right and the left-hand side, respectively. Objects-which we know are difficult to classify-with scores for comp 1 ranging between −0.5 and 1 are located in the center of the image. PCA can also be performed on each class separately, to investigate within-class variations, as shown in extensive Jupyter notebook examples in the RENT GitHub repository.
In addition to the PCA scores, the correlation loadings plot is shown in Fig. 7c, where every point represents one feature in the plane spanned by comp 1 and comp 2. The correlation loadings plot encodes 1) the level of contribution of the selected features to each of the components 1 and 2, and 2) how much of the variance in each feature is explained by the two components. The further away a correlation loading is located from the origin, the higher the amount of explained variance for the feature it represents. The inner and outer circles represent 50% and 100% of explained variances, respectively. Among the four selected features, feature 8, 21 and 28 contribute most to the first component that separates the two classes. It is also evident that these three features are highly correlated, as they are located so close to each other in Fig. 7c. Moreover, they are close to the outer ring, meaning that comp 1 explains nearly 100% of the variance in those features. Feature 22 contributes to both components 1 and 2, but is the feature that contributes most to component 2. By superimposing the scores onto the correlation loadings plots, we can gather information on how the scores and features are interrelated.  The above examples of post-hoc analysis illustrate how combining ensemble information with exploratory analysis by PCA, can provide deeper insight into the data.

Discussion
In summary, RENT performs well on all experimental datasets presented in this study when compared to the other feature selection methods. In particular, a good trade-off between predictive performance and selection stability is achieved. We observe that 1) RENT is consistently among the best performing methods, 2) if outperformed by others, the difference in performance is mostly negligible, and 3) the often lower number of features selected by RENT is a clear benefit. RENT does not fail for any of the presented datasets, whereas other methods show weaknesses on at least one dataset, with regard to either a very large number of selected features or poor predictive quality. In particular, RENT consistently performs well, whether the data are long-thin-many objects compared to the number of featuresor short-wide-relatively few objects compared to the number of features. The initial intention of RENT was to target short-wide datasets, which are particularly challenging when it comes to feature selection. In the presented evaluations, datasets c4 (Dexter text classification), c5 (OVA Lung), and r1 (Milk proteins) clearly fall into this category.
In addition to competitive performance, the number of features selected by RENT for the studied datasets is comparably low, which is a strength of RENT in terms of interpretability of results. Furthermore, the object-wise visualization demonstrated in Section 4.1 can provide previously unseen insights into the properties of the dataset, which may be particularly relevant for medical applications, but also for many other applications in general.
Robustness with regard to noisy data is another strength of RENT, which can be achieved by the extensive use of drawing subsamples from the training set. Particularly the baseline model M • , which is used as a benchmark in the experiments and achieves high performance on multiple datasets, is susceptible to poor initializations and hence, potentially less reliable for the selection of features. Although computationally more intensive than the comparing methods, RENT is less susceptible to poor initializations or convergence issues of optimization routines compared to other approaches.
In total, RENT has five model parameters to adjust by the user: two account for regularization intensity (γ and α) and three cutoffs control the strictness of feature selection (t 1 , t 2 and t 3 ). Both sets of hyperparameters are related, since a softer regularization allows a larger number of features, requiring higher cutoffs (and vice versa). Based on the presented parameter selection procedure using BIC, feature selectors which deliver a low number of features are favored in both stages.
In the current formulation, RENT is applicable for binary classification and regression problems. As introduced in [51], multiclass feature selection is not trivial and will be part of further research. However, a multiclass classification problem can be split into several binary problems, using schemes such as one-vs-one (OVO), one-vs-all (OVA), or error-correcting output coding (ECOC), as described in [52].

Conclusion
In this work, we presented a feature selection technique for binary classification and regression problems. The algorithm builds on the idea of training multiple elastic net regularized models on unique training data subsets. In particular, we define feature importance criteria based on the empirical distribution of feature-wise model weights.
Features are selected if their associated weights are regularly assigned high non-zero values with stable signs across the individual models of the ensemble.
We provided experiments on datasets from different disciplines, demonstrating that RENT is effective with respect to quantitative performance measures and interpretability and robustness. For the presented setups, the stability is very high even with a moderate number of ensemble models and in five out of six binary classification datasets, RENT achieves the highest MCC scores compared to the established feature selectors used in this study. For the regression datasets, RENT performed better or almost equal to the competing approaches. Further, we showed how to utilize information from the ensemble of models in a post-hoc analysis, advancing single-object interpretability.