Debiasing SHAP scores in random forests

Loecher, Markus Article — Published Version Debiasing SHAP scores in random forests AStA Advances in Statistical Analysis Provided in Cooperation with: Springer Nature Suggested Citation: Loecher, Markus (2023) : Debiasing SHAP scores in random forests, AStA Advances in Statistical Analysis, ISSN 1863-818X, Springer, Berlin, Heidelberg, Vol. 108, Iss. 2, pp. 427-440, https://doi.org/10.1007/s10182-023-00479-7 This Version is available at: https://hdl.handle.net/10419/308118 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ Vol.:(0123456789) AStA Advances in Statistical Analysis (2024) 108:427–440 https://doi.org/10.1007/s10182-023-00479-7 1 3 ORIGINAL PAPER Debiasing SHAP scores inrandom forests MarkusLoecher1 Received: 29 November 2022 / Accepted: 12 July 2023 / Published online: 22 August 2023 © The Author(s) 2023 Abstract Black box machine learning models are currently being used for high-stakes decision making in various parts of society such as healthcare and criminal justice. While tree-based ensemble methods such as random forests typically outperform deep learning models on tabular data sets, their built-in variable importance algorithms are known to be strongly biased toward high-entropy features. It was recently shown that the increasingly popular SHAP (SHapley Additive exPlanations) values suffer from a similar bias. We propose debiased or "shrunk" SHAP scores based on sample splitting which additionally enable the detection of overfitting issues at the feature level. Keywords Interpretable machine learning· Feature importance· Random forests· SHAP values· Explainable artificial intelligence 1 Introduction The unprecedented success of machine learning (ML) algorithms in diverse fields such as medicine, finance, biology, marketing, public health, image classification and natural language processing is truly revolutionary. In these applications it is often important to have models that are both accurate and interpretable, where being interpretable means that we can understand how the model uses the input features to make predictions. Given a data input, how did an algorithm arrive at the output? Which inputs influenced the decision of the algorithm? More specifically, why did the model, e.g., reject the loan application of a person? Which risk factors led to, e.g., the prediction of a high risk of heart disease? Many approaches have been proposed to address the lack of interpretability in ML. For an overview, see Molnar (2022). In bioinformatics, variable importance measures for random forests have gained popularity, particularly for selecting relevant genetic markers to predict diseases. Additionally, they have served as screening * Markus Loecher [email protected] 1 Department ofEconomics, Berlin School ofEconomics andLaw, 10825Berlin, Germany 428 M.Loecher 1 3 tools in significant applications, emphasizing the importance of reliable and wellunderstood feature importance measures. The overall goal of these approaches is to attribute importance to the features of a model according to their contribution to the model outcome. The importance of features in tabular data can be quantified at a local and global level. That is, how important is a feature for a particular decision (local) and how important is a feature for the overall decision-making process (global). Tree-based ensemble methods such as random forests and gradient boosted trees achieve state-of-the-art performance in many domains and have consistently outperformed deep neural network models on tabular-style datasets so common in medicine and public health (Grinsztajn etal. 2022). The default choice to measure global feature importances in most software implementations of random forests is the mean decrease in impurity (MDI). The MDI of a feature is computed as a (weighted) mean of the individual trees’ improvement in the splitting criterion produced by each variable. A well-known shortcoming of this default measure is its evaluation on the inbag samples which can lead to severe distortions of the reported ranking (Kim and Loh 2001; Strobl etal. 2007a). This bias is shown in the left panels of Fig.1 for a simulation design used by Strobl etal. (2007a) with details as follows. 1.1 Data withvarying cardinality A binary response variable Y is predicted from a set of 5 predictor variables that vary in their scale of measurement and number of categories. The first predictor variable X1 is continuous, while the other predictor variables X2,…,X5 are multinomial with 2, 4, 10,20 categories, respectively. The sample size for all simulation studies was set to n = 120. In the first null case, all predictor variables and the response are sampled independently. We would hope that a reasonable Fig. 1 Distributions of global feature importance scores (MDI and SHAP) for random classification forests using five features of different cardinalities (details of which are explained in the text). Upper row: Null simulation where no feature is informative; lower row: Power simulation where only feature X2 affects the outcome ( r=0.15 ). Both simulations illustrate the severe bias exhibited by MDI and—to a lesser extent— SHAP scores which inflate the importances of high cardinality features X1,3,4,5 MDI SHAP null power x1 x2 x3 x4 x5 x1 x2 x3 x4 x5 0 10 20 0 10 20 predictor score 429 1 3 Debiasing SHAP scores inrandom forests variable importance measure would not prefer any one predictor variable over any other. In the second simulation study, the so-called power case, the distribution of the response is a binomial process with probabilities that depend on a parameter r∈[0;0.5] which is related to the signal-to-noise-ratio (SNR) of x2 , namely P(y=1|X2=1)=0.5 −r,P(y=1|X2=2)=0.5 +r . For all figures, we repeatedly (100 times) fitted a random forest using the R1. ranger2. library with 400 classification trees and default parameters. SHAP values were computed using the R treeshap3. library. 1.2 Conditional feature contributions (CFCs) The conventional wisdom of estimating the impact of a feature in tree-based models is to measure the node-wise reduction of a loss function, such as the variance of the output Y, and compute a weighted average of all nodes over all trees for that feature. By its definition, such a mean decrease in impurity (MDI) serves only as a global measure and is typically not used to explain a per-observation, local impact. Saabas etal. (2019) proposed the novel idea of explaining a prediction by following the decision path and attributing changes in the expected output of the model to each feature along the path: Let f be a decision tree model, x the instance we are going to explain, f(x) the output of the model for the current instance, and fx(S)≈E[f(x)∣xS] the estimated expectation of the model output conditioned on the set S of feature values, then— following Lundberg etal. (2019)—we can define the Saabas value4 for the i’th feature as where Di x is the set of nodes on the decision path from x that split on feature i, and Aj is the set of all features split on by ancestors of j. Equation (1) results in a set of feature attribution values that sum up to the difference between the expected output of the model and the output for the current prediction being explained. When explaining an ensemble model made up of a sum of many decision trees, the CFCs for the ensemble model are defined as the sum of the CFCs for each tree. 1.3 SHAP values Lundberg etal. (2020) point out that CFCs are strongly biased to alter the impact of features based on their distance from the root of a tree. This causes CFC scores to be inconsistent, which means one can modify a model to make a feature clearly (1) 𝜙 s i(f,x)= ∑ j∈Di x fx(Aj∪i)−fx(Aj) , 1 https:// www.rproje ct. org/ 2 https:// github. com/ imbshl/ ranger. 3 https:// github. com/ Model Orien ted/ trees hap. 4 Synonymous with conditional feature contributions. 430 M.Loecher 1 3 more important, and yet the CFC attributed to that feature will decrease. The authors extend the original Shapley values from game theory to explain machine learning models based on SHapley Additive exPlanation (SHAP). For trees, these SHAP values can be computed efficiently in polynomial time. SHAP values can be seen as a solution to this problem, though often incurring potentially unnecessary computational cost as pointed out in Loecher (2022). As explained in Lundberg etal. (2019), Shapley values are computed by introducing each feature, one at time, into a conditional expectation function of the model’s output, fx(S)≈E[f(x)∣xS] , and attributing the change produced at each step to the feature that was introduced; then averaging this process over all possible feature orderings. Shapley values represent the only possible method in the broad class of additive feature attribution methods that will simultaneously satisfy three important properties: local accuracy, consistency, and missingness. They are defined as: where R is the set of all feature orderings, PR i is the set of all features that come before feature i in ordering R, and M is the number of input features for the model. Local accuracy states that for a specific input x, the explanation’s attribution values 𝜙i for each feature i need to sum up to the output f(x): In this paper we are mainly concerned with global importance measures, which is naturally defined as the average of the absolute SHAP values per feature across the data: All of our results pertain to marginal SHAP values. 1.4 Related work Marginal Shapley values are used to calculate the average contribution of a feature to the prediction when all possible coalitions of features are considered. This approach assumes that each feature is independent of the others, which may not be the case in reality. When features are highly correlated, the marginal Shapley value for one feature may be misleading as it doesn’t consider the influence of the other correlated features. On the other hand, conditional Shapley values consider the contribution of a feature based on a specific set of conditions. In other words, it takes into account the correlation among the features and calculates the average contribution of the feature based on the specific values of the correlated features. This approach provides a more discerning understanding of the contribution of each feature to the prediction. (2) 𝜙 i(f,x)= ∑ R∈R 1 M! [ fx(PR i∪i)−fx(PR i) ], (3) f(x)=E(f)+ M ∑ i=1 𝜙i(f,x ) (4) I j=1 N N ∑ i=1| 𝜙j(f,xi) | 431 1 3 Debiasing SHAP scores inrandom forests Janzing etal. (2020) offer a causal interpretation of Shapley values that replaces the conventional conditioning by observation with conditioning by intervention. Sundararajan and Najmi (2020) highlighted several drawbacks of the SHAP method, including the production of counterintuitive explanations when some features are deemed unimportant. However, the Baseline Shapley (BShap) method offers an improvement to this "uniqueness" issue in attribution methods. At first glance similar to our findings, Kwon and Zou (2022) also demonstrate that Shapley values often fail to sort features in order of influence on a model prediction. These attribution mistakes are especially pronounced when different marginal contributions have different signal and noise. The authors also propose a reweighing scheme as a solution which they refer to as WeightedSHAP. However, the underlying reasons for the incorrect ordering of features are very different from our rather narrow situation which combines the greedy tree splitting algorithm with varying feature entropies. Covert et al. (2020) introduce Shapley additive global importance (SAGE) for quantifying a model’s dependence on each feature. These global measures are the result of minimizing a structured loss function and hence much less heuristic than the ones defined by Eq. (4). Sutera etal. (2021) demonstrated that SAGE applied to tree-based models is very similar to the Gini importance. Supplement 1.2 contains limited results on SAGE values for our simulated data introduced in Fig.1. We refrained from delving deeper into this alternative method because (i) the provided code on the repository runs extremely slowly, and (ii) our focus are the original SHAP values which have become extremely popular. Additional approaches have been suggested to compute global SHAP values (Frye etal. 2020; Williamson etal. 2020; Casalicchio etal. 2019) where the latter reference proposes local and global feature importance methods based on the ideas of partial dependence and ICE which are lend themselves particularly well for visualization. Nevertheless, it is important to note that the inherent bias resulting from the underlying tree structure persists across all of these methods, provided that the tree-building algorithm remains unchanged. Yasodhara etal. (2021) conducted extensive experiments to assess the accuracy and stability of global feature importance scores, including SHAP and Gini. Their findings revealed notably low correlations with the true feature rankings, even in scenarios without additional noise. Moreover, when inputs or models were perturbed, these correlations further decreased. However, the authors did not specifically examine the impact of uninformative variables with different levels of cardinality. 2 Smoothed SHAP Less known than the distortions w.r.t MDI is the persistence of the bias of feature importance rankings based on (global) SHAP values, as shown in the right panels of Fig.1. Similar to MDI, SHAP shows a strong bias toward high-entropy variables which in the power case yields distorted rankings of variable impact. The grave consequences of providing such flawed explanations for, e.g., financial decisions or predictions of medical outcomes cannot be taken serious enough. To the best of 432 M.Loecher 1 3 our knowledge, no debiasing method has been successfully applied to SHAP scores (yet). 2.1 Sample splitting It was recognized early on that separating the split variable selection from the split point selection can reduce or eliminate the bias (Loh and Shih 1997; Hothorn etal. 2006). A somewhat similar approach is implemented in so-called “honest trees” (Athey and Imbens 2016). Honest trees are based on sample-splitting: Each node value is estimated using a different sub-sample of the training data than the one used to decide the split point. The sample-splitting successfully eliminates the split bias but comes with a cost of lower performance, as the number of samples for each estimation (split point and prediction) is reduced. It is somewhat surprising that the clear differentiation between validation and training data, which is typically observed for most other metrics in statistical modeling, has not been consistently applied to variable importance measures. Recently, several authors (Li etal. 2019; Loecher 2020b; Zhou and Hooker 2021) effectively eliminated the outlined bias inherent to the tree splitting procedure by including out- of-train samples in order to compute a debiased version of the MDI importance. These post-hoc adjustments of scores derived from tree-based models do not alter the tree structure, as opposed to attempts to prevent overfitting directly during the growth of the tree, as described in, e.g., Adler and Painsky (2022). The main idea is based on sample splitting in order to counteract the severe overfitting of high-cardinality features. Instead of evaluating feature importance scores solely on the same in-bag (IB) samples that were used to grow the trees, we penalize discrepancies with the same scores evaluated on out-of-bag (OOB) samples. While random forests—at least when fitted with the bootstrap—conveniently provide these IB/OOB without any "extra costs," we choose to generalize our presentation by splitting the data into train/test sets instead. This leads to a more model agnostic debiasing approach and hopefully simplifies the clarity of presentation. Whenever applicable, we show results for both train/test splits as well as IB/OOB splits. The initial steps of our proposed algorithm are as follows. 1. Split data into train/test sets. (Note the unusual terminology, since we are also "training" a model on "test" data) 2. Fit two separate random forest models on train/test each: RF train , RF test 3. Compute separate in/out SHAP values for, e.g., the test data (a) SHAP in : RF test predicts on test (b) SHAP out : RF train predicts on test We would expect weak/no correlations between SHAP in and SHAP out for noninformative features and in turn strong correlations for relevant features, which is confirmed by Fig.2. 433 1 3 Debiasing SHAP scores inrandom forests Based upon these two sets of in-sample and out-sample SHAP scores, we propose the following debiasing/smoothing/shrinkage scheme. For each feature j: 1 Fit a linear model SHAP j,out =𝛽j ⋅ SHAPj,in +u 2 Use the estimates � SHAPshrunk j,in = 𝛽j⋅SHAPj, in as local explanations instead of SHAPj,in , We can think of these predictions as "smoothed" version of the original SHAP values. We can avoid the need to fit two random forests by computing SHAP scores separately for each tree which was fitted on a specific bootstrapped subset (inbag). If the number of trees is sufficiently high, the union of these tree-wise disjoint inbag/ oob data sets will "generously" cover the entire original data, i.e., each row would have been part of the oob set sufficiently often. By averaging over all trees we hence obtain the equivalent of RF train and RF test SHAP scores with just one random forest. We expect the correlation between these inbag and oob SHAP values to be high for informative and nearly zero for uninformative features, which motivates the following equivalent but simpler smoothing method. For each feature j 1. Fit a linear model SHAPj,oob =𝛽j ⋅ SHAPj,inbag +u x4 x5 x1 x2 x3 −0.2 −0.1 0.0 0.1 0.2 −0.2−0.10.0 0.10.2 −0.2 −0.1 0.00.1 0.2 −0.10 −0.05 0.00 0.05 0.10 −0.10 −0.05 0.00 0.05 0.10 SHAP_in SHAP_out Fig. 2 In-sample versus out-sample SHAP scores for the power simulation where only feature X2 is informative ( r=0.25 ), illustrating the low correlations for features X1,3,4,5 and a strong relationship for X2 , with values of R2 1,2,3,4,5 =[0.038, 0.91, 0.004, 0.05, 0.044 ] respectively. (sample size n=240 ) 434 M.Loecher 1 3 2. Use the estimates � SHAPj,oob =  𝛽 j ⋅SHAP j,inbag as local explanations instead of the original SHAPj which mix inbag and oob values. The two right panels of Fig.3 display the main results: the non-informative features are shrunk toward zero while the relevant variable x2 is hardly modified. The benefits are most striking for low signal-to-noise ratio, e.g., the null simulation (top row) and for a moderate value of r=0.15 . We notice in passing that Kwon and Zou (2023) takes the idea of using oob data much further and propose them as a general data valuation method. We emphasize that our main goal so far has been to debias global SHAP scores rather than local explanations; hence the emphasis on relative rankings and the loss of an absolute interpretable scale. If local accuracy as expressed by Eq. (3) shall be preserved, we can simply rescale the scores by multiplying with (5) 𝛾 = M ∑ j=1 𝜙j(f,x)∕ M ∑ j=1  𝛽j⋅𝜙j(f,x) , raw shrunk_oob shrunk_test 00.15 0.25 x1x2x3x4x5x1x2x3x4x5x1x2x3x4x5 0.000 0.025 0.050 0.075 0.100 0.00 0.05 0.10 0.15 0.20 0.0 0.1 0.2 predictor SHAP Fig. 3 Raw vs. shrunk SHAP scores for the null (top row) as well as 2 power simulations where only feature X2 is informative ( r=0.15, 0.25 in middle and bottom row, respectively). The oob/inbag sample split (middle column) seems to lead to a similar debiasing effect as the computationally more expensive train/test method (right column)