FRI -- Feature Relevance Intervals for Interpretable and Interactive Data Exploration

TL;DR

FRI is an open source Python library designed to identify all relevant features in linear classification and regression tasks, aiding causal feature discovery in high-dimensional data and supporting interactive data exploration.

Contribution

The paper introduces FRI, a novel tool that provides comprehensive feature relevance intervals for linear models, enhancing interpretability and interaction in high-dimensional data analysis.

Findings

Enables identification of all relevant features in linear models.

Supports interactive and batch data exploration.

Facilitates discovery of causal features in biomedical data.

Abstract

Most existing feature selection methods are insufficient for analytic purposes as soon as high dimensional data or redundant sensor signals are dealt with since features can be selected due to spurious effects or correlations rather than causal effects. To support the finding of causal features in biomedical experiments, we hereby present FRI, an open source Python library that can be used to identify all-relevant variables in linear classification and (ordinal) regression problems. Using the recently proposed feature relevance method, FRI is able to provide the base for further general experimentation or in specific can facilitate the search for alternative biomarkers. It can be used in an interactive context, by providing model manipulation and visualization methods, or in a batch process as a filter method.

Tables

Table 1. TABLE I: Characteristics of simulated datasets. Each set consists of 30 features with 500 samples.

	Strongly relevant	Weakly relevant	Irrelevant
Sim1	4	4	22
Sim2	12	8	10
Sim3	4	0	26
Sim4	18	0	12
Sim5	0	20	10

Table 2. TABLE II: Feature selection score on simulated datasets. Values are showing the performance of each method to classify between the relevance of input features.

score	F1					precision					recall
data	Sim1	Sim2	Sim3	Sim4	Sim5	Sim1	Sim2	Sim3	Sim4	Sim5	Sim1	Sim2	Sim3	Sim4	Sim5
Boruta	0.98	0.82	0.91	0.82	0.98	0.99	1.00	0.87	1.00	1.00	1.00	0.72	0.98	0.70	0.95
EFS	0.96	0.76	0.71	0.84	0.94	0.93	1.00	0.57	1.00	1.00	1.00	0.62	0.98	0.73	0.90
ElasticNet	0.62	0.84	0.44	0.82	0.80	0.46	0.74	0.28	0.69	0.67	1.00	0.98	1.00	1.00	1.00
FRI	0.98	0.98	0.99	0.99	0.99	0.98	1.00	0.98	0.99	1.00	0.99	0.97	1.00	0.98	0.99
StabilitySelection	0.77	0.75	1.00	0.91	0.27	1.00	1.00	1.00	1.00	1.00	0.62	0.60	1.00	0.83	0.16

Table 3. TABLE III: Average training set accuracy. In the case of Boruta the internal RandomForest score was reported. For EFS accuracy is not defined.

	Boruta	EFS	ElasticNet	FRI	SS
Sim1	0.99	-	1.00	0.92	1.00
Sim2	0.97	-	1.00	0.96	1.00
Sim3	0.99	-	1.00	0.96	1.00
Sim4	0.97	-	1.00	0.93	1.00
Sim5	1.00	-	1.00	0.91	1.00
colp.	1.00	-	0.99	0.97	0.99
flip	1.00	-	0.90	0.82	0.90
spectf	1.00	-	0.99	0.92	0.98
t21	1.00	-	0.98	0.93	0.98
wbc	1.00	-	1.00	0.98	1.00

Table 4. TABLE IV: ROC-AUC values of logistic regression model using features selected by listed models. The values are averaged over 50 bootstraps.

	Boruta	EFS	EN	FRI	SS
colposcopy	0.568	0.586	0.640	0.661	0.625
flip	0.804	0.652	0.815	0.743	0.705
spectf	0.871	0.874	0.867	0.880	0.888
t21	0.971	0.977	0.971	0.975	0.978
wbc	0.997	0.998	0.998	0.998	0.999

Table 5. TABLE V: Average selected feature set size. Additionally for FRI the size of the strongly( s ) and weakly ( w ) relevant feature set is available.

	Boruta	EFS	EN	SS	FRI	FRIs	FRIw
Sim1	8.1	8.7	17.8	5.0	8.1	5.1	3.0
Sim2	14.3	12.3	26.6	12.1	19.4	12.4	7.0
Sim3	4.6	7.2	14.8	4.0	4.1	4.0	0.1
Sim4	12.6	13.2	26.2	15.0	17.9	17.9	0.0
Sim5	19.1	17.9	29.7	3.2	19.9	0.0	19.9
colp.	35.1	25.4	46.5	41.5	20.3	5.9	14.4
flip	18.8	8.1	16.9	9.1	8.9	8.8	0.1
spectf	44.0	20.3	43.1	5.9	19.9	5.9	14.0
t21	15.5	7.9	14.2	9.6	9.6	6.6	3.0
wbc	29.9	12.5	26.9	4.7	15.6	4.0	11.6

Equations

$$\left(\tilde{\omega },\tilde{b},\tilde{\xi }\right)\in$$

$$y_i({\omega }^{\top }{x}_i-b)\ge 1-{\xi }_i\forall i$$

$${\xi }_i\ge 0\forall i$$

$$(\mathbb{D},j):\underset{\omega ,b,\xi }{\min }\left|{\omega }_j\right|$$

$$y_i({\omega }^{\top }{x}_i-b)\ge 1-{\xi }_i,{\xi }_i\ge 0\forall i$$

$$\sum _{i=1}^n{\xi }_i\le \rho$$

$${\Vert \omega \Vert }_1\le (1+\delta )\cdot \mu .$$

$$(\mathbb{D},j):\underset{\omega ,b,\xi }{\max }\left|{\omega }_j\right|$$

$$\text{ constraints }(\text{ref}eq:optcons)\text{ hold.}$$

$$PI:={\overline{P}}_n\pm T_{n-1}(p)s_n\sqrt{1+(1/n)}.$$

$$K_l:=\left[K_l^{min},K_l^{max}\right].$$

$$(\mathbb{D},j,K):\underset{\omega ,b,\xi }{\min }\left|{\omega }_j\right|$$

$$\text{ constraints }(\text{ref}eq:optcons)\text{ hold and}$$

$$K_l^{min}\ge \left|{\omega }_l\right|\ge K_l^{max}\forall l\ne j.$$

$$F_1=2\cdot \frac{\mathrm{precision}\cdot \mathrm{recall}}{\mathrm{precision}+\mathrm{recall}}$$

Full text

††publicationid: pubid: 978-1-7281-1462-0/19/$31.00 ©2019 IEEE

1 Background

In recent years, due to an increasing availability of highly sensitive biotechnologies as well as an increasing digitalization of biomedical diagnostics, one could observe a trend towards bigger and more complex machine learning models. These models are used successfully in medical diagnoses [1][2] such as cancer prediction [3] using known biomarkers as input. When specific biomarkers are unknown, many learning models can also perform on nearly raw data without a preselection of variables. While such data representation often enables a high prediction accuracy, it is less suited if the purpose is data exploration and understanding of the underlying causal relationships. For the latter, insight into the model behavior and its relevant driving factors is necessary [4] and feature selection constitutes a first step to unravel the underlying relationships. But even for predictive models, the inference of sparse models can have a significant impact on the model performance, since it helps to mediate the curse of dimensionality thus leading to a better generalization ability and improved computational complexity. Because of this, there is a big need for methods to reveal relevant structure and function in the data itself.

Feature Selection (FS) constitutes one particularly prominent paradigm which enables the inference of sparse and interpretable prediction models [5]. The task is mostly undertaken in conjunction with model selection to improve predictive performance. Then the problem is often defined as finding the minimal subset of all features to achieve the best performance given some objective score. Commonly, one distinguishes filter, wrapper, and embedded methods. Filter methods are based on the information content of features as regards the output class label, disregarding the specific classification method. Hence they are particularly suited as a first screening technology for high dimensional data [6]. Wrapper approaches add (remove) features to a growing (shrinking) candidate set [7] while evaluating an inner model on a predefined metric. Embedded approaches use the internal model weights to find relevant features such as the Lasso [8] which enforces sparsity through its choice of regularization. Since they do not rely on iterative feature selection or weighting, embedded approaches have the benefit that they can effectively take into account interdependencies of groups of features. Further, they are specific to the used model.

In this contribution, we will focus on embedded methods for one particularly prominent group of models, namely linear classifiers or regression models, such as also addressed by the famous Lasso. While linear models can not applied to every problem, they are widely used in medicine [9] and the growing size of data makes their efficiency even more relevant. Furthermore, methods such as Lasso can be accompanied by formal guarantees on their ability to identify the true underlying features in the limit [10]. Yet, these guarantees do not hold if conditions as regards unique representability are violated. The latter is the case if features are high dimensional, highly correlated, or if there do exist different feature sets which yield the same prediction accuracy. In such cases, prediction accuracy in the limit still holds [11] yet the fact that Lasso can identify the true underlying features is violated. In particular, greedy and sparse approaches lead to instabilities of the selected feature set when correlated features are present. They remove similar features which could otherwise be grouped into potential functional units.

Finding a minimal subset (the usual objective of Lasso or related methods) is therefore not suited for giving a complete picture of all relevant features but only those which are sufficient to fulfill the prediction. To gain information about the importance of all features one needs to consider a more general measure: feature relevance[7][12]. The membership of the minimal feature set as a binary measure of relevance is too narrow in most cases. Kohavi et al. [7] coined the term all relevant feature selection (ARFS) in the 90s and added further distinction of relevance. Additional to strongly relevant features, which need to be part of the candidate feature set to achieve the best performance, and irrelevant features which are not beneficial for the considered relationship at all, they also introduced the concept of weakly relevant features i.e. features which can contribute information to a model, but are not necessarily included in every good model.

Weakly relevant features are often highly correlated with each other. They are not limited to pairs but also bigger groups of associated features. Out of these, at least one has to be part of the candidate set. Having knowledge about a group of features which could fulfill the same role in a model can be very important in the design of diagnostic tests where the source of data can differ by the cost or invasiveness of acquisition. Explicit redundancies of feature relevance then enable a practitioner to avoid, e.g., expensive features if they can be substituted by others. Additionally, feature groups could induce novel biological relationships. This is especially useful in gene co-expression or metabolomics experiments where groups of functional units are common [13]. Interestingly, it has been shown that extensions of Lasso which take into account feature correlations such as group Lasso can lead to wrong results due to a selection bias caused by feature correlations [14].

Finding an optimal solution to the ARFS problem is computationally intractable [15] but approximations exist. In 2005 the ElasticNet [16] overcame some of the instability of sparse models [17] by using a combination of multiple different regularization terms which lead to better conservation of weakly relevant features. The statistically equivalent signature (SES) [18] approach proposes a technology which groups mutually equivalent features into groups out of which minimal feature subsets can be constructed. Another proposal called stability selection uses resampling for more robust selection especially in high dimensional problems[19][20]. In 2010 Boruta [21] specifically focused on finding all relevant features using statistical testing of contrast features. Alternatives to the Boruta method are discussed and evaluated in [22], whereby Boruta was identified as best performing technology among the tested ones if used for different dimensionalities of the data. In 2017 Neumann et al. presented the Ensemble Feature Selection (EFS) method which combines multiple FS methods to remove individual biases and give aggregated feature relevance ranges for all of them [23, 24].

In this paper, we focus on the question of how to efficiently uncover a detailed view on the relevance of features in the case of possible feature redundancies. Thereby, we investigate linear models as a particularly relevant setting. The main contribution of this article is an accessible implementation of the feature relevance interval method (FRI) [25]111Source code and package available at github.com/lpfann/fri.. The FRI offers an efficient framework which assigns relevance intervals to features, rather than simple coefficient values. These intervals mirror the coefficient range of a feature when considering all possible models, hence offering detailed and complete information also in the case of feature redundancies. In this contribution, we offer an interactive software tool based on the mathematical framework as derived in [25] and we demonstrate its applicability in the context of biomedical data analysis. More precisely, we show that using these bounds, we can classify each feature into one of the three relevance classes. Furthermore, we can use these relevance bounds in visualizing the model which allows interpretability. Especially the information provided by the discrimination between strongly and weakly relevant features can help in biomedical applications. As an example we look at model design and biomarker discovery, where it could highlight elements which are crucial for the problem at hand while also providing alternative markers which have the same information.

Figure 1 displays the structure of our proposed software pipeline and in part this article. To elaborate: in Chapter 2 we shortly recapitulate the theoretical background and give details of the implementation of the method using linear programming. Specifically in 2.1.1 we show how to obtain a baseline solution and give our definition of relevance bounds in 2.1.2. We then describe how we can classify each feature into three relevance groups and how to reduce false positives using a probe based threshold estimation in 2.1.3. Then we show how we can constrain the use of features to certain relevance values (2.2) to facilitate interactive data exploration and model design. In 3.1 we evaluate the method quantitatively using simulated (3.1.1) and biomedical data (3.1.2). Finally we also provide an application focused qualitative evaluation and example in 3.2.

2 Implementation

2.1 Feature Relevance Intervals

For brevity, in the following definitions and experiments we only consider the task of binary classification, but the method and our implementation can be applied to linear regression and ordinal regression problems as well [26].

For a classification problem we observe data $\mathbb{D}$ as $(x_{1},y_{1}),\dots,(x_{n},y_{n})\ \in\mathbf{R}^{d}\times\{-1,1\}$ with $n$ different samples and $d$ real-valued features which have been tied to a target or response $y$. We assume that all $d$ features have been standardized at mean zero and standard deviation 1.

We propose an interactive pipeline, which enables the practitioner a detailed view of the relevance of features on the classification also in the case of multiple solutions or feature redundancies. More precisely, we propose an efficient framework which determines the relevance of features as regards all possible solutions of the classification problem, we demonstrate how this information can be used to identify all strongly as well as weakly relevant features, and we present a procedure to iteratively investigate feature combinations.

2.1.1 Baseline Solution

It is common practice to evaluate the relevance of a feature for a given classification by means of the weights assigned to the feature by a linear classifier such as a support vector machine [27]. Thereby, sparsity can be emphasized by resorting to, e.g., Lasso or sparse SVM models [8][28]. Yet, provided the solution is not unique, the resulting feature relevance is to some extent arbitrary, i.e. possibly rendering the model interpretation invalid. Here we propose an alternative: instead of taking the weights of a single linear model as an approximation of feature relevances we look into using the complete model class of well-performing sparse linear classifiers. A model class is characterized by its similarity in the quality of the solution i.e. similar generalization ability as characterized by the size of the weight vector of the SVM [29] and similar training loss, which measures wrongly classified samples. Through the use of a class of models, we can approximate the global solution of the ARFS problem as shown in [25].

Assume we are interested in linear classifiers of the form $y\mapsto sgn(\omega^{T}x-b)$ where $\omega$ is the normal vector of the separating hyperplane, $b$ denotes the bias and sgn refers to the sign function. First, we acquire a baseline solution to the problem using an $l_{1}$-regularized soft-margin SVM:

[TABLE]

Through optimization, we acquire a model fully defined by the normal vector of a hyperplane $\omega$ and its offset from the origin $b$. Prediction of samples is based on the signed distance from the plane. As usual, $\xi_{i}$ are slack variables to guarantee the feasibility of the constraint optimization problem in case of unavoidable classification errors. $\mathit{C}$ is a regularization parameter which depends on the datasets distribution. We choose the parameter guided by 3-fold stratified cross-validation and the $F_{1}$ score weighted by each class support to account for class imbalances.

From the model with the best $\mathit{C}$, we obtain constraints for controlling the generalization error of equivalent models: the upper limit on the weight vector $\mu:=\left\lVert\tilde{\omega}{}\right\rVert_{1}$ and error term $\rho:=\sum^{n}_{i=1}{\tilde{\xi}_{i}}$. These values determine the class of equivalent classifiers, which consist of all SVM solutions ($\omega^{\prime}$, $b^{\prime}$, $\xi^{\prime}$) such that $\left\lVert w^{\prime}\right\rVert_{1}\leq\mu$ and $\sum\xi^{\prime}_{i}\leq\rho$. All these alternatives are considered equivalent since they show the same performance for the given classification task as the found solution. Hence all weight vectors associated with an equivalent solution are relevant to determine the relevance of a feature to the given classification problem.

2.1.2 Minimum and Maximum Bounds

Using these constraints we now define feature relevance bounds for every feature $j$ independently i.e. we determine the interval of weight vectors which results if we take into account the weights of all possible equivalent classifiers. Mathematically speaking, we want to find extreme weight values for each feature given a similar error to our baseline. This can be done using linear programming.

For the lower bound, the lowest possible value of feature $j$, we define the problem

[TABLE]

And the upper bound for $j$ is defined as

[TABLE]

The optimization problems can be rewritten as linear optimization problems and solved in polynomial time [25] using appropriate solvers. To account for their numerical inaccuracies, which we encountered in our experiments, we also propose a relaxation factor $\delta=0.001$ in $(\ref{eq:optcons})$ to allow minor deviations from $\mu$. Given that all problems can be solved independently, our implementation makes use of parallel computation.

We end up with a real-valued matrix $\mathbb{RI}\in\mathbb{R}^{2\times d}$ which contains all pairs of relevance bounds. These intervals give an indication about the degree up to which a feature can or must be used in the classification and they can be visualized for interpretative purposes as seen in Fig. 2.

2.1.3 Feature Classification

By intuition one could use the following three rules to map feature relevances to their relevance class: {LaTeXdescription}

A feature is strongly relevant when its lower relevance bound is bigger than zero. The model class defined by its prediction accuracy, is dependent on information from it.

When two or more features are correlated, they can replace each other functionally in the model. These features are characterized by a lower bound equal to zero and an upper bound bigger than zero.

By definition, irrelevant features should have no measured relevance at all. Their upper (and lower) bound should therefore be zero.

In practice, the discrimination between relevant and irrelevant features is challenging. The use of slack variables in the overall model and thus our relevance bounds allow variation in the contribution of features. For relevance bounds specifically, even if feature $j$ is independent we observe $\text{maxRel}_{j}>0$. One could introduce a naive data independent threshold to discriminate between noise and relevant features, but this would lead to bad precision or recall of features in most cases. Instead we try to estimate the distribution of relevances of noise features given the model constraints. We expect for a given model class the same amount of allowed variation in the relevances and therefore a normal distribution with an unknown mean and variance. We propose to estimate the distributions parameters and the corresponding prediction interval (PI) to obtain a data dependent threshold[30]. To estimate this noise distribution we use $n$ permutated ($p$) input features from $\mathbb{D}$ i.e. for each $j^{p}$ we compute $\text{maxRel}((\hat{x}_{i},y_{i})^{n}_{i=1},j^{p})$ where $j\notin\hat{x}$.

The prediction interval is then defined as

[TABLE]

Here $\overline{P}_{n}$ denotes the sample mean and $s_{n}$ the standard deviation and $T$ represents the Student’s t-distribution with $n-1$ degrees of freedom. The size of PI depends on parameter $p$, the expected probability that a new value is included in the interval. We propose default values of $p=0.999$ for a low false positive rate and $n>=50$ to achieve a stable distribution.

To classify a feature $j$ as irrelevant we check if $\text{maxRel}_{j}~{}\in~{}PI$. To discriminate between weak and strong relevance we then additionally check the lower bound as described at the beginning of 2.1.3. An example of the feature classification can be seen in Fig. 2 where the relevance bars are colored according to their relevance class.

Our method provides the vector $\mathbb{AR}\in\{\mathit{0},\mathit{1},\mathit{2}\}^{d}$, which encodes strongly $(\mathit{2})$, weakly $(\mathit{1})$ and irrelevant $(\mathit{0})$ features.

2.2 Feature Constraints

The mathematical formalization of relevance intervals as introduced above also opens up the opportunity to integrate prior knowledge about feature relevances and to interactively explore good solutions by means of integrating additional bounds for specific features. More precisely by solving the problem using linear programs, the addition of constraints is very easy. One way we leverage this is the possibility of adding relevance constraints to the optimization.

Given data $\mathbb{D}$ and feature $l\in\mathbb{D}$ we define a set of additional constraint ranges $K$. For a constrained feature $l$ it is defined as

[TABLE]

Note that $K_{l}\geq 0$ because relevances are by definition positive and that $|K|\neq d$, i.e. not all features have to be constrained. In the case of $K_{l}^{min}=K_{l}^{max}$ we consider the models usage of feature $l$ as fixed to a static value.

Although the values in K could be chosen arbitrarily under the given restrictions, in practice one should orientate himself using the relevance bounds in $\mathbb{RI}$.

To compute relevance bounds including individual feature constraints, we have to extend the set of existing constraints in the optimization from Step 2 in Section 2.1. The minimum relevance bound with constraints is defined as

[TABLE]

The maximum bound is defined analogously.

To rewrite the new absolute term $\left\lvert\omega_{l}\right\rvert$ as a convex problem, we utilize the baseline solution $\tilde{\omega}$, which allows us to use the sign of the coefficient $\tilde{\omega}_{l}$ turning the non-convex absolute term into a simple convex one.

By changing the amount of contribution allowed for one feature, we can observe varying intervals for others and infer potential dependencies between them as in Fig. 2 (b). This is especially useful when designing a set of biomarkers. In our tool we provide the means to easily define ranges or values for all features. Presetting values can be based on the relevance intervals which give a first indication of importance. We also plan to integrate a function to automatically group features based on these constrained relevance bounds in the future.

3 Results

3.1 Feature Selection Evaluation

To evaluate the method in context, we run benchmarks against other established methods: Boruta as a representative of a method with statistical testing[31], Ensemble Feature Selection (EFS) [23], ElasticNet using an equal contribution of $L_{1}$ and $L_{2}$ regularization [16] and* stability selection* (SS)[19][20]. We removed the Lasso from the comparison because of very similar performance to the EN. For all methods, the proposed default parameters are used. Hyperparameters are selected according to a cross validation scheme. For the ElasticNet, we choose the feature set depending on the coefficients $c_{i}$ of the model where $c_{i}>10^{-5}$ counts as selected.

We used two types of sets: generated sets where we had knowledge of the underlying ground truth and real world data stemming from medical studies.

3.1.1 Simulation Data

All our simulation sets are sampled from a binary classification problem. To generate a multidimensional classification problem, we use a randomly generated prototype vector which defines a hyperplane. The defining features of this plane are strongly relevant. Now points are sampled in this feature space and the class is determined by the side of the hyperplane the points lie on. Weakly relevant features are constructed by replacing a feature of the original feature space with its linear combination. The elements of this combination are highly correlated and produce a set of redundant features. By removing the original feature and replacing it with those elements we achieve weak relevance by definition. Irrelevant features are sampled from a standard normal distribution. All simulation sets consist of 30 features and 500 samples. They differ in the density of the relevant feature space which is defined by the amount of strongly, weakly and irrelevant variables which are listed in Table I. According to these parameters 50 sets were generated per configuration and the following evaluation refers to the aggregated scores. Sim1 and Sim3 have a sparse relevant feature space while Sim2 and Sim4 are dense. Additionally, in Sim1 and Sim2 weakly relevant features are present, while they are missing completely in Sim3 and Sim4. Sim5 had all strongly relevant features removed.

Performance on these sets can give clues about the effectiveness of the considered feature selection strategy.

Due to a known ground truth, we can explicitly evaluate the validity of the selected features. We focus on the all-relevant feature selection problem and we use the following measures to evaluate the match of the detected feature set and the known ground truth of all relevant features: precision and recall. Recall is defined by TP / (TP+FN) with TP = true positives and FN = false negatives. It denotes how many of the relevant features were selected which is crucial when looking for the all relevant feature set. Precision is defined by TP / (TP+FP) with FP = false positives and describes what rate of false positives are part of the feature set. To get a balanced measure for a feature selection method, a combination of both is necessary. One can use the $F_{1}$ measure which is the harmonic mean of precision and recall:

[TABLE]

The measures can be seen in Table II.

Before we evaluate the selection measures we confirm that all models had a proper fit. Listed in Table III are the training accuracies. One can see in the table that all classification models had accuracy values over $90\%$ which signifies a sufficient fit of the data for all of them.

To evaluate the feature selection performance we mainly observe the $F_{1}$ score in Table II. Here our proposed FRI takes the lead overall with a near perfect score in all simulation sets. Depending on the presence of weakly relevant features, the other methods show loss of recall which leads to a reduced $F_{1}$ score. This is especially evident for Sim2 and Sim4 in the case of SS and EFS. The worst recall is achieved by SS for Sim5 where it did not select many of the weakly relevant variables at all. SS still achieves slightly better scores in Sim3 where no weakly relevant features are present.

3.1.2 Biomedical Data

The biomedical datasets are gathered from multiple studies and differ in size and type:

• •

t21: This set stems from a series of prenatal examinations of pregnant women. The goal of the examinations is the early diagnosis of chromosomal abnormalities, such as trisomy 21. The study covers sociodemographic, ultrasonographic and serum parameters which result in 18 usable features. The original set contains over 50.000 samples but only a low percentage ($\approx 0.8\%$) of abnormal samples. The data have been collected by the Fetal Medicine Centre at King’s College Hospital and University College London Hospital in London [32].

• •

flip: This set is used for the prediction of fibrosis. The diagnosis of fibrosis is represented as a score which is based on sociodemographic and serum parameters. The set consists of 118 patients and 19 features. The data were provided by the Department of Gastroenterology, Hepatology and Infectiology of the University Magdeburg [33].

• •

spectf: The spectf dataset consists of 44 features describing cardiac Single Proton Emission Computed Tomography (SPECT) images. Each of the 267 patients images were diagnosed as either normal or abnormal.

• •

wbc: The wbc dataset contains 32 markers for cell image based breast cancer diagnostics from 569 patients.

• •

colposcopy: Set with 69 Extracted structural features from videos acquired during colposcopies [34]. Classificiation of practitioners clinical judgment using the Schiller modality.

spectf, wbc and colposcopy were acquired through the UCI Machine Learning Repository [35].

The biomedical datasets are preprocessed before analysis. Samples with over $90\%$ missing values are removed. Sets are split into stratified training and testing subsets. If samples still contain missing feature values, we replace them with the features training set mean in both subsets. Similarly, the z-score transformation is based on the training set and applied to both. In case the original set is imbalanced, we use the Synthetic Minority Over-sampling Technique (Smote) [36] in combination with a Nearest Neighbor cleaning rule [37][38] as described in [39]. In one case (t21) with an extremely large majority class, we only perform downsampling. We perform the tests on 50 bootstrap replicates with sample size $0.7n$ with replacement.

To asses the quality of the feature selection on real-world datasets, we have to rely on the problem performance itself since no ground truth as regards feature relevance is available. We expect a FS method to pick features which contain information and a loss of features with information is signified in a decrease of performance. Instead of looking at each models internal accuracy score, we evaluate the selected feature sets by their discriminative power, whereby the latter is uniformly evaluated by a logistic regression model, which is a very popular model for predictive purposes in medical applications [40]. This model is trained using only the predicted feature set and hyperparameter optimization is performed to select the model with the best cross-validated regularization parameter. Finally, for this selected model the Receiver operation characteristics (ROC) on the holdout validation set are recorded. For the comparison, we look at the area under the curve (ROC-AUC), which is listed in Table IV.

Here the AUC on the five datasets produces no clear overall superior method which is in line with the expectation that the minimal optimal set is the objective of most methods and sufficient for prediction. On the spectf, t21 and wbc datasets most methods produce very similar performing feature sets. In the case of colposcopy the feature set selected by FRI achieves the best performance. SS produces sightly better sets in two cases. The ElasticNet performs solid in all cases based on its very conservative selection method where informative features are not removed often.

This leads the part of the evaluation more concerned with the goal stated in the introduction. In the search for an interpretative and complete feature set we need to take the selected set size into account. Table V lists the average feature set sizes over all experiments. Because FRI provides additional information by not only conserving all weakly relevant features but also by denoting the feature class itself we can explicitly list those as well. As mentioned in the last paragraph, we can easily see that EN is very conservative in its selection. It produces by far the biggest feature sets with many false positives in the case of the Sim sets but also most likely in the real datasets. Similarly for Boruta, which achieves better precision in the generated data but still shows seemingly inflated set sizes. SS on the other hand exhibits very good precision overall. Interestingly the size of the sets chosen by stability selection is very similar to FRIs, the set of strongly relevant features chosen by FRI. This indicates that FRI can find strongly relevant features with high precision, but also highlights the additional information provided by the weakly relevant features contained in FRIw.

But why do we need the information in FRIw?

3.2 Evaluation of Interactive Use

By having additional information available in the set of weakly relevant features FRIw we can gain insights into the structure of the data. We can improve the design of models and diagnostic tests in biomedical applications. Our framework given in 2.2 allows introducing constraints into the model. This makes it possible to limit the contribution of certain features to specific intervals or a fixed value. These limits can come from prior knowledge of the practitioner and represent design goals or existing hypotheses. Depending on the chosen values the model and the resulting relevance bounds change and can be visualized again which lends itself to an iterative and interactive process. In the following, we are going to evaluate that use case on simulated data and on the t21 data set.

The simulated set was generated according to 3.1.1. It consists of 8 features, 4 of which are strongly relevant, 3 of which are weakly relevant and one noise feature. Fig. 3 (a) shows the output of FRI without any constraints. The four strongly relevant features (1-4) are visible as four small rectangles with lower relevance bounds (the bottom part of the rectangle) bigger than zero. The model parameters allow some variation in their contribution to the model. The three weakly relevant features (5-7) are visible as three taller rectangles with equal height because they are perfectly correlated in the normalized space. They can replace each other in the model. This is apparent when we preset one of them (e.g feature 5) to the minimum and maximum relevance bound i.e. we calculate $\text{minRel{C}}(\mathbb{D},j,K)$ and $\text{maxRel{C}}(\mathbb{D},j,K)$ for all $j\neq 5$ and $K=\{K_{5}\}$. $K_{5}$ is then either $\mathbb{RI}^{min}_{5}$ or $\mathbb{RI}^{max}_{5}$, i.e fixed to a static value.

In Fig. 3 (b) feature 5 was set to the minimum bound and the relevance bounds of other features are identical. That is because the other two features are still a correlated pair which allows the same degree of variability in contribution. When feature 5 is set to its maximum relevance bound in (c) we see that feature 6 and 7 have no contribution anymore. Additionally, all other relevance bounds are reduced to single values because the model in this state does not allow any more variability.

After this experiment, it is now interesting to apply this procedure on real data with functional associations between features. In Fig. 2 (a) the normal FRI output of the t21 set is presented. As a reminder, this set consists of samples acquired in prenatal examinations of mothers and their unborn children. Features in the study included socioeconomic factors as well as ultrasound imaging metrics. Notably in the output of FRI are two weakly relevant features 1 and 6. Feature 1 represents the gestational age of the fetus in days (‘GA-d’) and feature 6 the crown rump length (‘CRL’) of the fetus, which is the length as indicated on an ultrasound machine. By intuition, we expect an association between the two measures. If we set one of the two features to its minimum relevance bound (Fig. 2 (b)) we see that feature 6 becomes strongly relevant in the model. This highlights the association between the two which is very useful in cases where it is not clear a priori. Furthermore, we can use this as a design tool to easily select ‘better’ features. If we find functional alternatives, we can exclude more expensive features in future experiments or tests.

3.3 Runtime

The computational runtime of a method is not only an indicator for its feasibility on bigger datasets but also especially important in the interactive use case where an analyst is actively involved in model refinement. Because relevance bounds can be solved independently in parallel, we provide the means to speed up computation by utilizing all available CPU cores on the machine. Additionally, we also tested running our program in conjunction with the distributed computation framework Dask which allows scaling up to any amount of separate computing nodes in a high performance cluster such as Grid Engine or even in cloud backends.

In Fig. 4 we display aggregated mean runtime of the methods used in our evaluations on a single CPU thread. Because only FRI and one other method provided a parallel implementation (SS) the computations were limited to one thread, so an advantage of the parallel processing is not taken into account. EN performed best followed by SS. Both show steady runtimes over all types of data. Boruta’s runtime is very dependent on the density of the feature space and shows some variance in the case of t21. The runtime of FRI is similar to Boruta in most cases but takes a hit in smaller datasets because of the constant factor of sampling permutated features for feature classification. EFS shows the slowest performance in most cases, which clearly stems from its use of multiple complex underlying models at the same time.

4 Conclusion

We have presented the software library FRI to produce all relevant feature sets for general feature selection as well as perform interactive data exploration. We described how we implemented the algorithm from [25] and extended the method to allow a practitioner to include new constraints and experiment. We also proposed a threshold estimation method to reduce false positives which are common in all-relevant selection tasks.

In comparison with other methods, we showed that FRI can detect all relevant features in synthetic datasets while minimizing noise through its threshold estimation. On real datasets we showcased good selection performance and additional information provided by the weakly relevant feature set. Our underlying method ensures to conserve all relevant variables while still maintaining interpretability. This is facilitated by the three relevance classes our method produces as well as the relevance bar representation which should enable better understanding for biological and medical experts in the future. In addition to facilitating understanding we also provide a way to incorporate prior knowledge to manipulate the model itself which should help in the design of new experiments and biomarkers for prediction models.

Acknowledgments

The authors would like to thank Professor Kypros Nicolaides and Dr. Argyro Syngelaki from Fetal Medicine Foundation for making available the t21 data used in this study. We are also grateful for the flip datasets provided by Professor Ali Canbay, Department of Gastroenterology, Hepatology, and Infectiology of the University Hospital Magdeburg.