2.1 Feature importance

With respect to random forests, three types of feature importance scores are well known in the literature. The first one is an impurity-based feature importance. The so-called impurity is quantified by the splitting criterion of the collection of contained decision trees. Therefore, it is likely to overestimate the importance of large numerical features (if the dataset is not standardized). Furthermore, it is possible that features that may not be predictive on unseen data are found to be important in the case of overfitting. For these reasons, a second type of feature importance, the so-called permutation importance, has found its way into the literature and is to be preferred. It is defined as the decrease in model performance when a single feature’s values are randomly shuffled. A similar possibility to measure a feature’s importance is based on exclusion of a variable. More precisely, to measure the importance of a feature, a second random forest is trained on the same data excluding the feature under consideration. The importance is, again, measured as the decrease in the model’s performance. Of course, the permutation-based and the exclusion-based approach have their shortcomings if there are clusters of (highly) correlated features (Breiman 2001). Indeed, due to the availability of the modified feature’s values via the correlation, the importance score can be highly misleading. One approach to overcome this problem, which is often used in the process of feature extraction, is to keep only one variable per cluster (Chen et al. 2020; Louppe 2014; Guyon et al. 2002). If the ultimate goal is to design a decent prediction model with as few features as possible, this is the state of the art. But in some cases, researchers are actually more interested in estimating the importance of each feature to determine which features influence the dependent variable and how strongly. In this setting, it may be convenient to treat the correlations differently. There are at least two decorrelation techniques that are usually used either for clustering data or for designing well-performing prediction models: the Gram-Schmidt decorrelation technique (Zhang & Chan 2006) or residual-based decorrelation (Dezfouli et al. 2019). The main idea in both cases is to subtract the information from a given feature F_i given by F₁,…, F_i-1, F_i+1,… F_d and use this residual to train the model.

2.2 Trends

We compare three different ways to define trends in the dataset. The simplest way one might think of examining a trend between the values of a feature Y and the predicted variable is to use the correlation coefficient r(X,Y) = Cov(X,Y)/(σ(X)σ(Y)), which reflects linear trends. A more general correlation coefficient that handles any monotone trends are various types of rank correlation coefficients such as the Spearman correlation coefficient ρ(X,Y) = r(R(X), R(Y)) where R(·) denotes the rank function. This method of finding trends is well established and only considers the observable raw data.

Another approach does not look at the raw data but fits a model and looks for trends in that model. Many practitioners tend to identify trends in multivariate tasks by fitting a linear model to the data and interpreting the sign and corresponding p-value of the coefficient of a feature as a trend. We will denote this coefficient by r_LM(F). However, we will see that this can be very misleading, even for very simple datasets.

In recent years, an old concept from mathematical game theory, called Shapley values, has been used to interpret machine learning models (Lundberg & Lee 2017). In particular, they are well understood mathematically for tree-based models and random forests. The Shapley value of a feature with respect to a data point measures how much the feature value contributes to the prediction compared to the average prediction and is defined as the average marginal contribution of the feature value among all possible combinations of features. For a formal definition, see Shapley’s original paper (1953), and for a detailed discussion of how to use the concept in machine learning, see Janzing et al. (2020) and Sundararajan and Najmi (2020). Shapley scores can be used to determine trends. More precisely, if a low value of a feature X induces a decrease of the Shapley value, then the model predicts a negative trend of the predictive variable with respect to feature X.

2.3 Studied datasets

To test the performance of our estimators in practice, we use two very well-known real datasets, called Kaggle fish market dataset (FISH) (Pyae 2019) and California housing data (HOUSING) (Nugent 2017). In addition, we create multiple different synthetic datasets to explore certain aspects of the estimators.

FISH contains the records of seven different common fish species in fish market sales. The features are species, weight, vertical length, diagonal length, transverse length, height, and width for each fish. Of these characteristics, we used weight, height, and width to predict vertical length. The California housing data refers to the houses found in a given California county and summary statistics based on 1990 census data. The features are longitude, latitude, median age of the house, total number of rooms, total number of bedrooms, population, number of households, median income, and ocean proximity for each county with median house value as the prediction target. We transformed the ocean proximity feature into an ordinal scale.

The first synthetic datasets (SYN1(a) and SYN1(b)) are derived from a base dataset B consisting of 1000 samples and 10 features, three of which are informative. The base dataset is standardized by removing the mean and scaling to unit variance. It is then combined with a Gaussian noise dataset N standardized the same way and with the same structure but no informative features, resulting in dataset SYN1(a). For SYN1(b), the same procedure is applied using uniform white noise rather than Gaussian noise. A family of datasets is obtained by:

SYN1(a) and SYN1(b) each consist of the combination of the base dataset with 250 different random noise datasets. These datasets are used to compare the robustness of trend estimators.

The second synthetic dataset (SYN2(a)) consists of 100 samples with independently generated features:

Furthermore, given X₀, we define:

The true label is given by:

Thus, the real labels depend on X₀, X₁, X₂ and A₀, A₁, A₂ can be considered as noisy instances of X₀ with different types of dependencies. Similarly, we define a set SYN2(b), in which the standard Gaussian 𝒩(0, 1) is replaced by a uniform [–1,1] random variable. This is an instance of white noise applied to data. SYN2(a) and SYN2(b) are used to compare different notions of feature importance.

The last synthetic datasets (SYN3(a and b)) consist of 100 samples with only one informative feature X₀, defined as previously. Moreover, A₀, A₁, A₂ are defined as above. The true labels are now given by $Y=4X_0^{1.5}$M11 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[Y = 4X_0^{1.5}\] \end{document} . SYN3 (a and b) are used to compare the notions of feature importance on a cluster of correlated features, in direct comparison to SYN2, in which two additional, uncorrelated, informative features are present.

3.1 Finding trends in a dataset

We compare the commonly used regression coefficients r, ρ, the linear model-based trend estimator, a Shapley-based trend estimator, and propose two novel estimators based on random forest regression to determine the trends of features. For this purpose, we simply define the Shapley-based trend of a feature as the correlation between its values X and its Shapley values s(X), so that we obtain the estimators r(X, s(X)) and ρ(X, s(X)), respectively. Moreover, we report the Shapley values themselves visually to compare the actual marginal contributions to the corresponding trend estimators.

The two proposed trend estimators are the absolute and the relative traversal rate. The random forest regression model uses an ensemble of uncorrelated decision trees. At each node, the current dataset is partitioned into two partition classes based on the values of the node’s feature. We assume without loss of generality that the data in the left partition class belong to small feature values and the data in the right partition class belong to large feature values. To determine the trend of a feature F, we compare the mean of the features in the left and right partition classes per node. If the average value of the predicted variable in the left tree is smaller than in the right tree, this corresponds to a positive correlation with the feature F. The intuition behind this is the following. As we find only data points with a smaller value of F in the left tree than in the right tree, and the model predicts a smaller value of the predicted variable, we expect a positive correlation between the feature F and the dependent variable. More formally, let {F_j}_{j =1…n} denote the set of nodes in the random forest in which the data is partitioned with respect to feature F. The corresponding partition classes are called L(F_n) and R(F_n). If the feature F is clear from the context, we abbreviate these classes to L_n and R_n. Furthermore, for a subset A of the values of the predicted variable, we define $\mathrm{AVG}(A)=|A{|}^{-1}\underset{a\in A}{\Sigma }a$M12 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[AVG(A) = |A{|^{ - 1}}\mathop \Sigma \limits_{a \in A} a\] \end{document} as the average value of the set A. This allows us to define our trend estimators.

Definition 1

Given a random forest$ℛ$M13 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{\cal R}\] \end{document} , let {F_j}_{j= 1 …n} denote the set of nodes in the random forest in which the data is split with respect to feature F. The absolute traversal rate of feature F is defined as:

$\mathrm{ATR}(ℛ,F)=n^{-1}{\displaystyle \sum _{i=1}^n\left(1\{\mathrm{AVG}(L(F_i))\le \mathrm{AVG}(R(F_i))\}-1\{\mathrm{AVG}(L(F_i))\ge \mathrm{AVG}(R(F_i))\}\right)}.$

M5 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ATR({\cal R},F) = {n^{ - 1}}\sum\limits_{i = 1}^n {\left({{\bf{1}}\{ AVG(L({F_i})) \le AVG(R({F_i}))\} - {\bf{1}}\{ AVG(L({F_i})) \geq AVG(R({F_i}))\} } \right)}.\] \end{document}

Moreover, the relative traversal rate of feature F is defined as:

$\mathrm{RTR}(ℛ,F)=2{\displaystyle \sum _{i=1}^n\frac{R(F_i)-L(F_i)}{|L(F_i)+R(F_i)|}}.$

M6 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[RTR({\cal R},F) = 2\sum\limits_{i = 1}^n {\frac{{R({F_i}) - L({F_i})}}{{{\rm{|}}L({F_i}) + R({F_i})|}}}.\] \end{document}

The ATR formalizes the idea that we have a trend when the feature with a higher value causes the model to return a higher value. It is a number between –1 and 1. If it takes a value close to 1, this means that in almost every split, the smaller feature’s values induce a smaller model prediction, while a value close to –1 means that in most splits the smaller feature’s values induce a larger model prediction. A value around zero means that about half of the splits induce a positive, and half of the splits a negative correlation.

The RTR is a weighted variant of this idea. If the relative difference between the average prediction in the left and right sub-tree is large, the split can be said to induce a stronger trend than in the case in which the average predictions are close to each other in both trees.

Both ATR and RTR estimate the impact of a feature by assessing how it splits the tree into partition classes (see Figure 1). In particular, the impact of a feature is estimated by comparing those partition classes, which can be done without knowing any of the other features. This is in contrast to the Shapley-based trend estimator, which assesses the impact of a feature by calculating the average marginal contribution of a feature over different coalitions of features.

Figure 1 

Each occurrence of feature F splits the dataset into two parts. In the example, F₁ creates partition classes L₁ = {2, 4, 7, 3} and R₁ = {8, 12, 4, 6}. The split at F₂ creates classes L₂ = {7} and R₂ = {3}, whereas the split at F₃ defines L₃ = {8, 12} and R₃ = {4, 6}. The model is agnostic to any features other than F.

Trend Estimator. For the empirical analysis of the datasets, we wrote a trend estimator module, which trains both a random forest regression model as well as a linear model on the input data. For each feature, the trend estimator then outputs

1. the coefficient of the linear model
2. Spearman’s rank correlation and Pearson correlation coefficient between the target and
   1. the Shapley value
   2. the feature
3. the absolute traversal rate (ATR)
4. the relative traversal rate (RTR).

3.2 Measures of Importance

We compare five different notions of feature importance. Two definitions, the impurity-based feature importance and the permutation-based feature importance, are well studied concepts (Breiman 2001). We use scikit-learn’s default implementation of these measures. Moreover, the exclusion-based feature importance is a folklore variant to measure a feature’s importance by the model’s decrease in accuracy by excluding a feature. To compare different notions of feature importance with the exclusion-based variant, we denote by MSE(R, F) the mean squared error of the random forest model in which feature F was excluded. We then calculate MSE(R, F)/MSE(R) as a measure of feature importance. We train 50 random forest models R₁, …, R₅₀ independently on 80% of the data each and estimate Ê(F) as the mean over MSE(R_i,F)/MSE(R_i) (for i = 1 … 50). Finally, we define the exclusion-based feature importance of features (F₁, …, F_d) as the normalised value of Ê(F), thus E(F) = Ê(F)/(Ê(F₁) + … + Ê(F_d)).

In addition, we introduce two novel types of feature importance based on residual learning. The idea is that the importance of feature F_i is determined by its residuals given features F₁, …, F_i-1, F_i+1,… F_d. With a slight misuse of notation, we interpret F_i ∈ ℝⁿ as the vector of all values corresponding to feature F_i and denote by Y ∈ ℝⁿ the values of the dependent variable. We denote by ${\mathcal{A}}_j$M14 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{{\cal A}_j}\] \end{document} an arbitrary algorithm that takes F₁, …, F_j-1 as input and outputs a vector in ℝⁿ. Given a fixed permutation π of [d], we denote by i₁, …, i_d the new order under π. To determine the importance of F_i, we determine its importance under all permutations π with the property that i_d = i and weight it by the performance of a model consisting only of the feature F_i. The interpretation is as follows: given all other features, what can be learned from feature $F_{i_d}$M15 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{F_{{i_d}}}\] \end{document} ? The algorithm to compute the importance can now be expressed as follows.

• For all permutations π which map i ⟼ d, do the following
  • – Define $W_1^{\pi }=F_{i_1}$M16 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[W_1^\pi = {F_{{i_1}}}\] \end{document} .
  • – Replace the values of feature $F_{i_j}$M17 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[F_{i_{j}}\] \end{document} with $W_{i_j}^{\pi }=F_{i_j}-{\mathcal{A}}_{{\rangle }_{|}}(W_{i_1}^{\pi },\dots ,W_{i_{j-1}}^{\pi })$M18 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[W_{{i_j}}^\pi = {F_{{i_j}}} - {{\cal A}_{{\rangle _|}}}(W_{{i_1}}^\pi, \ldots, W_{{i_{j - 1}}}^\pi)\] \end{document} (for j = 2 … d).
  • – Train a random forest with features $\left\{W_{i_j}\right\}$M19 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[\{ {W_{{i_j}}}\}\] \end{document} .
  • – Determine the impurity-based feature importance of $W_{i_d}^{\pi }$M20 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[W_{{i_d}}^\pi\] \end{document} .
• Determine the average feature importance of F_i as the mean over all $W_{i_d}^{\pi }$M21 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[W_{{i_d}}^\pi\] \end{document} , call this (FI)_i.
• Train a random forest regressor ${ℛ}_i$M22 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{{\cal R}_i}\] \end{document} with feature F_i and dependent variable Y and measure $r(ℛ(F_i),Y)$M23 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[r({\cal R}({F_i}),Y)\] \end{document} .
• Return ${\tilde{f}}_i=r({ℛ}_i(F_i),Y){(FI)}_i$M24 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{\tilde f_i} = r({{\cal R}_i}({F_i}),Y){(FI)_i}\] \end{document}

After applying this algorithm, we are left with ${\tilde{f}}_1,\dots ,{\tilde{f}}_d$M25 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{\tilde f_1}, \ldots, {\tilde f_d}\] \end{document} . Finally, we define the feature importance based on the residual algorithm $\mathcal{A}$M26 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{\cal A}\] \end{document} as the standardized version of the above estimator, namely:

Formally, the algorithm is given as Algorithm 1. We note the following.

• $f_j(\mathcal{A})$M32 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{f_j}({\cal A})\] \end{document} is a random quantity because it depends on the training of the random forest regressors ${ℛ}_1,\dots ,{ℛ}_d$M33 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{{\cal R}_1}, \ldots, {{\cal R}_d}\] \end{document} and the random forest regressors using the features $\left\{W_{i_j}\right\}$M34 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[\{ {W_{{i_j}}}\}\] \end{document} .
• In applications, it may not be possible to iterate over all permutations π. Instead, the average impurity-based feature importance is estimated by sampling some permutations.
• The algorithm is highly dependent on the residual algorithm $\mathcal{A}$M35 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{\cal A}\] \end{document} .

In this contribution, we empirically analyze the feature importance based on two different residual algorithms: classical residual learning by random forest regression and decorrelation by the Gram-Schmidt method.

4.1 Finding Trends

In the following, we report our empirical results on the performance of the different trend estimators on the HOUSING, FISH, and SYN1 datasets.

4.1.1 SYN1

The datasets SYN1(a) and SYN1(b) were used to test the robustness of the different trend estimators with respect to the mixing of the dataset with noise. To do this, the trend estimator module was applied to D_w for each w ∈ [0.01, …, 1.] for each of the 250 random noise datasets. In our experiment, the aggregated output shows that both ATR and RTR, as well as the Shapley correlation, are more robust than the linear model for the informative features for both Gaussian noise and uniform white noise (Figure 2). The Shapley values are the most robust, followed by RTR and ATR. Interestingly, non-informative features were also assigned large ρ_S, r_S, RTR and ATR values.

Figure 2 

Mean and 95% confidence interval for the different trend estimators on SYN1(a) and SYN(b) for 250 independent trials each. On the x-axis, the proportion of noise is reported. Features 1–3 are informative, whereas features 4–10 are non-informative.

4.1.2 FISH

We performed three experiments on the FISH dataset using the features weight, height and width to predict length. All three selected features are positively correlated with the target (Figure 3).

Figure 3 

Pairplot of the used fish market dataset features (weight, height and width) and the predicted variable (Length).

First, we applied the trend estimation module to the FISH dataset. To control for random effects, we performed 100 bootstrapping iterations, sampling from a subset of 70%. The linear regression model assigned a negative coefficient to the Height feature, while the other trend estimators reported a positive trend (Figure 4).

Figure 4 

Comparison of the trend estimators for FISH. We report the mean and the standard deviation of the different trend estimators over 100 bootstrap iterations, each containing 70% of the data. Relative absolute SHAP values shows the absolute sum of the SHAP values for each run, divided by the highest respective sum.

To evaluate the robustness of the trend estimators to noise, we used a random mixing strategy similar to that used to create SYN1. The FISH data were standardized and mixed with random noise ranging from 0% to 99% noise before being used as input to the trend estimator module. We found that the linear model and the RTR became unstable as the feature-to-target correlation r and ρ decreased, while the ATR and the Shapley measures r_S and ρ_S remained relatively unaffected up to much higher mixing rates (Figure 5).

Figure 5 

Mean and 95% confidence interval w.r.t. 100 independent iterations over noise on FISH. The x-axis reports the proportion of noise mixed to the real data.

4.1.3 Housing

The random forests were trained 100 times on HOUSING, each run using a random initialization and a random subset containing 70% of the training data. The mean and standard deviation of the different trend estimators, as well as the relative absolute SHAP values is shown in Figure 6, i.e. the sum of the absolute SHAP values divided by the maximum absolute sum for each run. The characteristic population is very weakly negatively correlated with the housing price. However, all trend estimators report a significant negative trend for population. The feature total rooms is positively correlated with the target. However, the linear model assigns a negative coefficient to the total number of rooms. All other trend estimators report a positive trend.

Figure 6 

Comparison of the trend estimators on HOUSING. The linear model assigns a negative coefficient to the total number of rooms feature, even though the feature itself is positively correlated to the target.

4.2 Measures of importance

We compare the impurity-based feature importance, the permutation-based feature importance, the exclusion-based feature importance and the feature importance induced by the two described residual algorithms (residual learning and Gram-Schmidt decorrelation). A run consists of fitting a random forest, to measure the feature importance values, multiple independent runs (between 100 and 400) were conducted. To determine the residual-based importance scores, for each feature, 50 permutations that assign this feature to the last position are sampled independently in every run.

First, we compare the different scores on the synthetic datasets SYN2 and SYN3 (see Figure 7). Perhaps the most important observation is that the impurity-based feature importance assigns the same score to all features – in both datasets. This is in strong contrast to all other feature importance scores. It is noteworthy that both residual-based approaches produce very comparable scores on the given datasets. Both residual-based approaches and the permutation-based score assign roughly the same score to X₀ and the slightly noisy variant A₀. However, feature A₁, which is subject to much more noise, receives a significantly higher score under residual-based scoring. Especially on SYN3, the residual-based approaches assign a not too small score to all informative features. The permutation-based score for the informative features X₁ and X₂ is comparatively small. However, all scores assign a higher importance to the noisy instance A₀ of X₀ than to the informative features X₁ and X₂.

Next, we compare the different scores on the real datasets HOUSING and FISH (see Figure 8). For HOUSING, it is most striking that the residual learning, impurity-based, and permutation-based scores assign the largest value to the median income, followed by the proximity to the ocean and the latitude/longitude, while the Gram-Schmidt-based score assigns only a large value to the median income and all other features receive comparable scores. In addition, the population is found to be more important by the impurity-based and permutation-based approaches as opposed to the residual learning-based approach.

For FISH, all measures assign the highest score to weight and all measures assign a non-vanishing score to all three variables. However, the score of width is within one standard deviation in the residual learning-based, impurity-based, and permutation-based approaches. Only the Gram-Schmidt-based score assigns a significantly larger value to weight and considers height to be the second most important feature.