Przemysław Biecek 
Susanne Dandl
Giuseppe Casalicchio
Marvin N. Wright
The increasing availability of data and software frameworks to create predictive models has allowed the widespread adoption of machine learning in many applications. However, high predictive performance of such models often comes at the cost of interpretability: Many trained models by default do not provide further insights into the model and are often too complex to be understood by humans such that questions like ‘’What are the most important features and how do they influence a prediction?’’ cannot be directly answered. This lack of explanations hurts trust and creates barriers to adapt predictive models, especially in critical areas with decisions affecting human life, such as credit scoring or medical applications.
Interpretation methods are valuable from multiple perspectives:
In this chapter, we focus on some important methods from the field of interpretable machine learning that can be applied post-hoc, i.e. after the model has been trained, and are model-agnostic, i.e. applicable to any model without any restriction to a specific model class. Specifically, we introduce methods implemented in the following three packages:
iml presented in Section 11.2,counterfactualspresented in Section 11.3 that implements counterfactual explanation methods, andDALEX presented in Section 11.4.The iml and DALEX packages offer similar functionality, but they differ in design choices. Both iml and counterfactuals are based on the R6 class system, thus working with them is more similar in style to the mlr3 package. In contrast, DALEX is based on the S3 class system and focuses on comparing multiple predictive models, usually of different types, on the same plot. This is a model comparison-oriented approach, often referred to as Rashomon perspective. We will briefly introduce the available methods in the respective sections. Most of the methods are described in more detail in the introductory book about interpretable machine learning by Molnar (2022).
Throughout this chapter, we use the rchallenge::german data, which contains corrections proposed by Groemping (2019). The data includes 1000 credit applications from a regional bank in Germany with the aim to predict creditworthiness, labeled as “good” (which is the positive class) and “bad” in the column credit_risk. For the sake of clarity and simplicity, we use only half of the 20 features:
library("mlr3")
task = tsk("german_credit")
task$select(cols = c("duration", "amount", "age", "status", "savings", "purpose",
"credit_history", "property", "employment_duration", "other_debtors"))
task<TaskClassif:german_credit> (1000 x 11): German Credit
* Target: credit_risk
* Properties: twoclass
* Features (10):
- fct (7): credit_history, employment_duration, other_debtors,
property, purpose, savings, status
- int (3): age, amount, durationWe want to fit a random forest model using mlr3. We use 2/3 of the data to train the model and the remaining 1/3 to analyze the trained model using different interpretation methods. A detailed discussion on which data set to use for interpretation is given in Section 11.5.
library("mlr3learners")
set.seed(1L)
split = partition(task)We fit a "classif.ranger" learner, i.e. a classification random forest, on the training data to predict the probability of being creditworthy (by setting predict_type = "prob"):
learner = lrn("classif.ranger", predict_type = "prob")
learner$train(task, row_ids = split$train)iml PackageThe iml package (Molnar, Bischl, and Casalicchio 2018) implements a variety of model-agnostic interpretation methods and is based on the R6 class system just like the mlr3 package. It provides a unified interface to the implemented methods and facilitates the analysis and interpretation of machine learning models. Below, we provide examples of how to use the iml package using the random forest model fitted above with the mlr3 package.
The iml package supports machine learning models (for classification or regression) fitted by any R package. To ensure that this works seamlessly, fitted models need to be wrapped in an iml::Predictor object to unify the input-output behavior of the trained models. An iml::Predictor object contains the prediction model as well as the data used for analyzing the model and producing the desired explanation. We construct the iml::Predictor object using the test data which will be used as the basis for the presented model interpretation methods:
The constructor iml::Predictor$new() has an (optional) input argument predict.function which requires a function that predicts on new data. For models fitted with the packages mlr3, mlr or caret, a default predict.function is already implemented in the iml package and the predict.function argument is not required to be specified. For models fitted with other packages, the model-specific predict function of that package is used by default. Passing a custom predict.function to unify the output of the model-specific predict function might be necessary for some packages. This is especially needed if the model-specific predict function does not produce a vector of predictions (in case of regression tasks or classification tasks that predict discrete classes instead of probabilities) or a data.frame with as many columns as class labels (in case of classification tasks that predict a probability for each class label).
When deploying a model in practice, it is often of interest which features are contributing the most towards the predictive performance of the model. On the one hand, this is useful to better understand the problem at hand and the relationship between the features and the target to be predicted. On the other hand, this can also be useful to identify irrelevant features and potentially remove those using feature importance methods as the basis for feature filtering (see Section 5.1). In the context of this book, we use the term feature importance to describe global methods that calculate a single score per feature reflecting the importance regarding a certain quantity of interest such as the model performance. The calculated feature importance is reported by a single numeric score per feature and allows to rank the features according to their importance.
It should be noted that there are also other notions of feature importance not based on a performance measure such as the SHAP feature importance, which is based on Shapley values (see Section 11.2.5 and Lundberg, Erion, and Lee (2019a)) or the partial dependence-based feature importance introduced in Greenwell, Boehmke, and McCarthy (2018), which is based on the variance of a PD plot (see Section 11.2.3).
One of the most popular methods is the permutation feature importance (PFI), originally introduced by Breiman (2001a) for random forests and adapted by Fisher, Rudin, and Dominici (2019) as a model-agnostic feature importance measure, which they called model reliance. The PFI quantifies the importance of a feature regarding the model performance (measured by a performance measure of interest). Specifically, the PFI measures the change in the model performance before and after permuting a feature (the former refers to the original model performance). Permutation means that the observed feature values in a dataset are randomly shuffled. This destroys the dependency structure of the feature with the target variable and all other features while maintaining the marginal distribution of the feature. The intuition of PFI is simple: if a feature is not important, permuting that feature should not affect the model performance. However, the higher the measured change in model performance is, the more important a feature is deemed. Therefore, the PFI summarizes the change in model performance after destroying the information of each feature due to permutation in a numeric feature importance score (that either reflects the difference or the ratio of the change in performance). In general, repeating the permutation process and aggregating the individual importance values is recommended because the results can be unreliable due to the randomness of the permutation process. The PFI method requires a labeled dataset (for which both the feature and true target are given), a trained model, and a performance measure to quantify the decrease in model performance after permutation.
We now analyze which features were the most important ones according to the PFI for our fitted random forest model to classify creditworthiness. We initialize an iml::FeatureImp object with the model and the classification error ("ce") as the performance measure and visualize the results with the $plot() method.
importance = FeatureImp$new(predictor, loss = "ce")
importance$plot()
By default, iml::FeatureImp repeats the permutation 5 times, and in each repetition, the importance value corresponding to the change in the classification error is calculated. The $plot() method of the iml::FeatureImp object shows the median of the 5 resulting importance values (as a point) and the boundaries of the error bar in the plot refer to the 5 % and 95 % quantiles of the importance values.
The number of repetitions can be set with the input argument n.repetitions of iml::FeatureImp$new(). By default, iml::FeatureImp uses the ratio of the model performance before and after permutation as an importance value. Alternatively, the difference instead of the ratio of the two performance measures can be used by setting compare = "difference" in iml::FeatureImp$new().
We see that status is the most important feature. That is, if we permute the status column in the data, the classification error of our model increases by a factor of around 1.17.
Feature effect methods describe how or to what extent a feature contributes towards the model predictions by analyzing how the predictions change when changing a feature. For example, these methods are useful to identify whether the model estimated a non-linear relationship between a feature of interest and the target. In general, we distinguish between local and global feature effect methods. Global feature effect methods refer to how a prediction changes on average when a feature is changed. In contrast, local feature effect methods address the question of how a single prediction of a considered observation changes when a feature value is changed. To a certain extent, local feature effect methods can also reveal interactions in the model that become visible when the local effects are heterogeneous, i.e., if changes in the local effect are different across the observations.
A popular global method to visualize feature effects is the partial dependence (PD) plot (Friedman 2001). It visualizes how the model predictions change on average if we vary the feature values of a certain feature of interest. Later, Goldstein et al. (2015) discovered that the global feature effect as visualized by a PD plot can be disaggregated into local feature effects associated with a single observation. Specifically, they introduced individual conditional expectation (ICE) curves, a visual tool for local feature effects, and showed that the PD plot is just the average of ICE curves. ICE curves display how the prediction of a single observation changes when varying a feature of interest while all other features stay constant (ceteris paribus). Hence, for each observation, the values of the feature of interest are replaced by multiple other values (e.g., on an equidistant grid of the feature’s value range) to inspect the changes in the model prediction. Specifically, an ICE curve visualizes on the x-axis the set of feature values used to replace the feature of interest and on the y-axis the prediction of the model after the original feature value of the considered observation has been replaced. Hence, each ICE curve is a local explanation that assesses the feature effect of a single observation on the model prediction. A heterogeneous shape of the ICE curves, i.e., the ICE curves are not parallel and might cross each other, indicates that the model may have estimated an interaction involving the considered feature.
Feature effects are very similar to regression coefficients \(\beta\) in linear models which offer interpretations such as ‘’if you increase this feature by one unit, your prediction increases on average by \(\beta\) if all other features stay constant’’. However, feature effects cannot only convey linear effects but also more complex ones (similar to splines in generalized additive models) and can be applied to any type of predictive model.
Now, we inspect how the feature amount influences the predictions of the creditworthiness classification. For this purpose, we compute feature effects using PD plots and ICE curves. We generally recommend accompanying PD plots with ICE curves because showing the PD plot alone might obfuscate heterogeneous effects and interactions. To plot PDP and ICE, we initialize an iml::FeatureEffect object with the feature name and select method = "pdp+ice". Again, we used the $plot() method to visualize the results:
effect = FeatureEffect$new(predictor,
feature = "amount", method = "pdp+ice")
effect$plot()
As we have a binary classification task at hand, we see PD and ICE plots for the predicted probabilities of both predicted classes separately. The plot shows that if the amount is smaller than roughly 10,000 DM (standing for Deutsche Mark, the currency in Germany before the Euro was introduced as cash in 2002), there is on average a high chance that the creditworthiness is good.
Interpretable models such as decision trees or linear models can be used as surrogate models to approximate or mimic a (often very complex) black-box model. Inspecting the surrogate model provides insights into the behavior of a black-box model. For example, the (learned) model structure of a decision tree, or the (learned) parameters (coefficients) of a linear model can be easily interpreted and provide information on how the features formed the prediction.
We differentiate between local surrogate models, which approximate a model locally around a specific data point of interest, and global surrogate models that approximate the model across the entire input space (Ribeiro, Singh, and Guestrin 2016; Molnar 2022). The surrogate model is usually trained on the same data used to train the black-box model or on data having the same distribution to ensure a representative input space. For local surrogate models, the data used to train the surrogate model is usually weighted according to the closeness to the data point of interest. In either case, the predictions obtained from the black-box model are used as the target to train a surrogate model that approximates the black-box model.
With iml::TreeSurrogate, we can fit a tree-based surrogate model, more specifically a partykit::ctree() model from the partykit package, to the predictions of our random forest model we want to analyze.
treesurrogate = TreeSurrogate$new(predictor, maxdepth = 2L)
treesurrogate$plot()
statusis either "0 <= ... < 200 DM", "... >= 200 DM" or "salary for at least 1 year") and either a duration of less or equal than 36 months (first node), or more than 36 months (second node). The third and fourth nodes contain applicants that either have no checking account or a negative balance (status). For the third node the additional requirement was a duration of less or equal to 21 months, for the fourth node it was a duration of more than 21 months.As we have used maxdepth = 2, the surrogate tree performs two binary splits yielding 4 leaf nodes. iml::TreeSurrogate extracts the decision rules found by the tree surrogate (that lead to the leaf nodes), and the $plot() method visualizes the distribution of the predicted outcomes for each leaf node. For example, we can see that the leaf node shown in the top left bar chart only consists of creditworthy applicants and results from the split by duration <= 36 and by status being either 0<= ... < 200 DM or ... >= 200 DM / salary for at least 1 year (i.e., meaning that the balance of the credit applicant’s checking account is positive).
We can access the trained tree surrogate via the $tree field of the fitted iml::TreeSurrogate object and can, therefore, apply any functions provided by the partykit package to the model object, e.g., partykit::print.party() to show the performed splits:
partykit::print.party(treesurrogate$tree)[1] root
| [2] status in no checking account, ... < 0 DM
| | [3] duration <= 21: *
| | [4] duration > 21: *
| [5] status in 0<= ... < 200 DM, ... >= 200 DM / salary for at least 1 year
| | [6] duration <= 36: *
| | [7] duration > 36: *Since the surrogate model only uses the predictions of the black-box model (here, the random forest model) and not the real outcomes of the underlying data, the conclusions drawn from the surrogate model do not apply generally, but only to the black-box model (if the approximation of the surrogate model is accurate).
To evaluate whether the surrogate model approximates the prediction model accurately, we can compare the predictions of the tree surrogate and the predictions of the black-box model. For example, we can produce a cross table comparing the predicted classes of the random forest and predicted classes of the surrogate tree:
pred_surrogate = treesurrogate$predict(credit_x, type = "class")
pred_rf = learner$predict_newdata(credit_x)$response
prop.table(table(pred_rf, pred_surrogate$.class))
pred_rf bad good
good 0.11515152 0.67575758
bad 0.13636364 0.07272727mean(pred_rf == pred_surrogate$.class)[1] 0.8121212Mostly, the black-box predicted class and the surrogate predicted classes overlap (with around 81% matching predictions).
In general, it can be very difficult to accurately approximate the black-box model with an interpretable one in the entire feature space. However, it is much simpler if we only focus on a small area in the feature space surrounding a specific point – the point of interest. For a local surrogate model, we conduct the following steps:
To illustrate this using the iml package, we first select a data point we want to explain. Here, we select an observation in the data set, let us call him Steve.
steve = credit_x[249, ]
steve age amount credit_history duration
249 74 5129 critical account/other credits elsewhere 9
employment_duration other_debtors property purpose
249 >= 7 yrs none real estate car (new)
savings status
249 unknown/no savings account ... < 0 DMFor Steve, the model predicts the class bad with 54.6% probability:
predictor$predict(steve) good bad
1 0.4537167 0.5462833As a next step, we use iml::LocalModel, which fits a locally weighted linear regression model that explains why Steve was classified as having bad creditworthiness by inspecting the estimated coefficients. The underlying surrogate model is L1-penalized such that only a pre-defined number of features per class, say k, will have a non-zero coefficient. The default is k = 3L and shows the three most influential features.
The argument gower.power can be used to specify the size of the neighborhood for the local model (default is gower.power = 1). The smaller the value, the stronger the local model will focus on points closer to the point of interest (here: Steve). If the prediction of the local model and the prediction of the black box model greatly differ, you can change the value of gower.power and, e.g., increase the number of non-zero coefficients k for a more accurate local model.
predictor$class = "good"
localsurrogate = LocalModel$new(predictor, steve, gower.power = 0.01)
localsurrogate$results[, c("beta", "effect", "x.recoded")] beta
amount -1.223788e-05
credit_history=critical account/other credits elsewhere -1.042568e-01
status=... < 0 DM -4.376759e-02
effect x.recoded
amount -0.06276809 5129
credit_history=critical account/other credits elsewhere -0.10425677 1
status=... < 0 DM -0.04376759 1From the table, we see that the three locally most influential features are status, amount, and credit history. A closer look at the coefficients reveals that all three features have a negative effect.
Shapley values were originally developed in the context of cooperative game theory to study how the payout of a game can be fairly distributed among the players that form a team. This concept has been adapted for use in machine learning as a local interpretation method to explain the contributions of each input feature to the final model prediction of a single observation (Štrumbelj and Kononenko 2013). The analogy is as follows: In the context of machine learning, the cooperative game played by a set of players refers to the process of predicting by a set of features for a single observation whose prediction we want to explain. Hence, the features are considered to be the players. The total payout, which should be fairly distributed among the players, refers to the difference between the individual observation’s prediction and the mean prediction.
In essence, Shapley values aim to answer the question of how much each input feature contributed to the final prediction for a single observation (after subtracting the mean prediction). By assigning a value to each feature, we can gain insights into which features were the most important ones for the considered observation. Compared to the penalized linear model as a local surrogate model, Shapley values guarantee that the prediction is fairly distributed among the features.
Shapley values are sometimes misinterpreted: The Shapley value of a feature does not display the difference of the predicted value after removing the feature from the model training. The exact Shapley value of a feature is calculated by considering all possible subsets of features and computing the difference in the model prediction with and without the feature of interest included. Hence, it refers to the average marginal contribution of a feature to the difference between the actual prediction and the mean prediction, given the current set of features.
With the help of iml::Shapley, we now generate Shapley values for Steve’s prediction. Again, the results can be visualized with the $plot() method.
If we focus on the plot, the Shapley values (phi) of the features show us how to fairly distribute the difference of Steve’s probability to be creditworthy to the data set’s average probability to be creditworthy among the given features. Steve’s duration of 9 months has the most positive effect on the probability of being creditworthy, with an increase in the predicted probability of about 5 %. The most negative effect on the probability of being creditworthy has the credit history, with a decrease in the predicted probability of more than 15 %.
counterfactuals PackageCounterfactual explanations try to identify the smallest possible changes to the input features of a given observation that would lead to a different prediction (Wachter, Mittelstadt, and Russell 2017). In other words, a counterfactual explanation provides an answer to the question: “What changes in the current feature values are necessary to achieve a different prediction?”. Such statements are valuable explanations because we can understand which features affect a prediction and what actions can be taken to obtain a different, more desirable prediction.
Counterfactual explanations can have many applications in different areas such as healthcare, finance, and criminal justice, where it may be important to understand how small changes in input features could affect the model’s prediction. For example, a counterfactual explanation could be used to suggest lifestyle changes to a patient to reduce their risk of developing a particular disease, or to suggest actions that would increase the chance of a credit being approved. For the rchallenge::german data, a counterfactual could be: “If the credit applicant would have asked for a lower credit amount and duration, he would have been classified being of good credit risk with a probability of > 50 %” (Figure 11.1). By revealing the feature changes that alter a decision, the counterfactuals reveal which features are the key drivers for a decision.
Many methods were proposed in previous years to generate counterfactual explanations. These methods differ in what targeted properties their generated counterfactuals have (for example, are the feature changes actionable?) and with which method (for example, should a set of counterfactuals be returned in a single run?). The simplest method to generate counterfactuals is the What-If approach (Wexler et al. 2019). For a data point, whose prediction should be explained, the counterfactual is equal to the closest data point of a given data set with the desired prediction.
Another method is the multi-objective counterfactuals method (MOC) introduced by Dandl et al. (2020). Compared to the What-If method, MOC can generate multiple counterfactuals that are not equal to observations in a given dataset but are artificially generated data points. The generation of counterfactuals is based on an optimization problem that aims for counterfactuals that
All these four objectives are optimized simultaneously via a multi-objective optimization method.
Due to the variety of methods, counterfactual explanations were outsourced into a separate package, the counterfactuals package, instead of integrating these methods into the iml package. Currently, three methods are implemented in the R package – including What-If and MOC – but the R6-based interface makes it easy to add other counterfactual explanation methods in the future.
To illustrate the overall workflow of the package, we generate counterfactuals for Steve using the What-If approach (Wexler et al. 2019). The counterfactuals package relies on iml::Predictor() as a model wrapper and can, therefore, explain any prediction model fitted with the mlr3 package, including the random forest model we trained above.
predictor$predict(steve) good
1 0.4537167The random forest classifies Steve with 45.4% as being good creditworthy. The What-If method can be used to answer, e.g., how the features need to be changed such that Steve is classified as good with a probability of more or equal to 60%. Since we have a classification model we initialize a counterfactuals::WhatIfClassif() object. By calling $find_counterfactuals(), we generate a counterfactual for Steve.
library("counterfactuals")
whatif = WhatIfClassif$new(predictor, n_counterfactuals = 1L)
cfe = whatif$find_counterfactuals(steve,
desired_class = "good", desired_prob = c(0.6, 1))cfe is a counterfactuals::Counterfactuals object which offers many visualization and evaluation methods. For example, $evaluate(show_diff = TRUE) shows how the features need to be changed.
data.frame(cfe$evaluate(show_diff = TRUE)) age amount credit_history duration
1 -8 -3603 all credits at this bank paid back duly 3
employment_duration other_debtors property purpose savings
1 <NA> <NA> <NA> <NA> <NA>
status dist_x_interest no_changed dist_train dist_target
1 no checking account 0.242221 5 0 0
minimality
1 4For a probability of more than 60 % for being creditworthy, the age is reduced by 8 years and the amount is reduced by 3603 DM while the duration is enlarged by 3 months. Furthermore, Steve’s credit history must be improved (he must pay back all his credits at his bank duly, instead of having a critical account/other credits elsewhere), and his status must be "no checking account". For the employment duration, other debtors, property, purpose, and savings no changes are required. Additional columns in cfe$evaluate() reveal the quality of the counterfactuals, for example, the number of required feature changes (no_changed) or the distance to the closest training data point (dist_train) which is 0 because the counterfactual is a training point.
Calling the MOC method is very similar to calling the What-If method. Instead of a counterfactuals::WhatIfClassif() object, we initialize a counterfactuals::MOCClassif() object. We set the epsilon parameter to 0 to penalize counterfactuals in the optimization process with predictions outside the desired range. With MOC, we can also prohibit changes in specific features, here age, via the fixed_features argument. For illustrative purposes, we let the multi-objective optimizer only run for 30 generations.
library("counterfactuals")
set.seed(123L)
moc = MOCClassif$new(predictor, epsilon = 0, n_generations = 30L,
fixed_features = c("age"))
cfe_multi = moc$find_counterfactuals(steve,
desired_class = "good", desired_prob = c(0.6, 1))Since the multi-objective approach does not guarantee that all counterfactuals have the desired prediction, we remove all counterfactuals with predictions not equal to the desired prediction via the $subset_to_valid() method.
cfe_multi$subset_to_valid()
cfe_multi5 Counterfactual(s)
Desired class: good
Desired predicted probability range: [0.6, 1]
Head:
age amount credit_history duration
1: 74 5129 all credits at this bank paid back duly 9
2: 74 5129 existing credits paid back duly till now 9
3: 74 5129 no credits taken/all credits paid back duly 9
6 variables not shown: [employment_duration, other_debtors, property, purpose, savings, status]Overall we generated 5 counterfactuals. Compared to the counterfactual of the What-If approach, MOC’s counterfactuals were artificially generated and are not necessarily equal to actual observations of the underlying data set. For a concise overview of the required feature changes, we can use the plot_freq_of_feature_changes() method. It visualizes the frequency of feature changes across all returned counterfactuals.
cfe_multi$plot_freq_of_feature_changes()
We see that credit history and status are the only changed features. To see how the features have been changed, we can visualize the counterfactuals on a 2-dim ICE plot, also called a surface plot.
cfe_multi$plot_surface(feature_names = c("status", "credit_history"))
The colors and contour lines indicate the predicted value of the model when credit history and status differ while all other features are set to the values of Steve. The white point displays Steve, and the black points are the counterfactuals that only propose changes in the two features. We see that the credit histories are better or equal to Steve’s value of having a critical account/other credits elsewhere (changes to all credits at this bank paid back duly, existing credits paid back duly till now, no credits taken/all credits paid back duly) and that the status is one time lowered to no checking account, three times stayed on … < 0 DM and one time increased to … >= 200 DM / salary for at least 1 year.
DALEX PackageThe DALEX (Biecek 2018) package implements a similar set of methods as the iml package presented above, but the architecture of DALEX is oriented towards model comparison. The logic behind working with this package assumes that the process of exploring models is an iterative process, and in successive iterations we want to compare different perspectives, including perspectives presented/learned by different models. This logic is commonly referred to as the Rashomon perspective, first described in “Statistical Modeling: The Two Cultures” paper (Breiman 2001b) and more extensively developed and formalized as interactive explanatory model analysis (Baniecki, Parzych, and Biecek 2023).
You can use the DALEX package with any number of models built with the mlr3 package as well as with other frameworks in R. The counterpart in Python is the library dalex (Baniecki et al. 2021). In this chapter, we present code snippets and a general overview of this package. For illustrative purposes, we reuse the german_credit task built in the Section 11.1 on rchallenge::german data. Since we already presented intuitions behind most of the methods in Section 11.2 we do not duplicate the descriptions below and focus on presenting the syntax and code snippets.
Once you become familiar with the philosophy of working with the DALEX package, you can also use other packages from this family such as fairmodels (Wiśniewski and Biecek 2022) for detection and mitigation of biases, modelStudio (Baniecki and Biecek 2019) for interactive model exploration, modelDown (Romaszko et al. 2019) for the automatic generation of IML model documentation in the form of a report, survex (Krzyziński et al. 2023) for the explanation of survival models, or treeshap for the analysis of tree-based models.
The analysis of a model is usually an interactive process starting with a shallow analysis – usually a single-number summary. Then in a series of subsequent steps, one can systematically deepen understanding of the model by exploring the importance of single variables or pairs of variables to an in-depth analysis of the relationship between selected variables to the model outcome. See Bücker et al. (2022) for a broader discussion of what the model exploration process looks like.
This explanatory model analysis (EMA) process can focus on a single observation, in which case we speak of local model analysis, or for a set of observations, in which case we speak of global model analysis. Below, we will present these two scenarios in separate subsections. See Figure 11.2 for an overview of key functions that will be discussed. An in-depth description of this methodology can be found in the EMA book (Biecek and Burzykowski 2021).
Predictive models in R have different internal structures. To be able to compare models created with different frameworks, an intermediate object – a wrapper – is needed to provide a consistent interface for accessing the model. Working with explanations in the DALEX package always starts with the creation of such a wrapper with the use of the DALEX::explain() function or its specialized variants. For models created with the mlr3 package, it is more convenient to use the dedicated DALEXtra::explain_mlr3() version. Let us use it to wrap the learner object.
library("DALEX")
library("DALEXtra")
ranger_exp = DALEXtra::explain_mlr3(learner,
data = credit_x,
y = as.numeric(credit_y$credit_risk == "bad"),
label = "Ranger Credit",
colorize = FALSE)Preparation of a new explainer is initiated
-> model label : Ranger Credit
-> data : 330 rows 10 cols
-> target variable : 330 values
-> predict function : yhat.LearnerClassif will be used ( default )
-> predicted values : No value for predict function target column. ( default )
-> model_info : package mlr3 , ver. 0.16.0 , task classification ( default )
-> predicted values : numerical, min = 0.005031746 , mean = 0.3008911 , max = 0.8951405
-> residual function : difference between y and yhat ( default )
-> residuals : numerical, min = -0.8951405 , mean = -0.0008911013 , max = 0.9577865
A new explainer has been created! The DALEXtra::explain_mlr3() function performs a series of internal checks so the output is a bit verbose. If any of the internal checks is alarming then the explain function output will show a yellow WARNING message. If any check does not pass, such as the predict function or the residuals cannot be calculated then a red ERROR message will appear. Turn the verbose = FALSE argument to make it less wordy.
ranger_expModel label: Ranger Credit
Model class: LearnerClassifRanger,LearnerClassif,Learner,R6
Data head :
age amount credit_history duration
1 35 6948 no credits taken/all credits paid back duly 36
2 28 1403 no credits taken/all credits paid back duly 15
employment_duration other_debtors property
1 1 <= ... < 4 yrs none building soc. savings agr. / life insurance
2 1 <= ... < 4 yrs none building soc. savings agr. / life insurance
purpose savings status
1 car (new) unknown/no savings account ... < 0 DM
2 others unknown/no savings account no checking accountThe result of the explain function is an object of the class explainer in the S3 system. This object is a list containing at least: the model object, the dataset that will be used for calculation of explanations, the predict function, the function that calculates residuals, name/label of the model name and other additional information about the model. The generic print function for this object outputs three items: the name of the model (a label given by the user), information about the model’s classes (in this case, an object of the LearnerClassifRanger class which inherits from several broader classes), and the first few rows of data used for the explanation count.
The global model analysis aims to understand how a model behaves on average on a set of observations. In the DALEX package, functions for global level analysis have names starting with the prefix model_.
As shown in Figure 11.2, the model exploration process starts by evaluating the performance of a model. This can be done with a variety of tools, but in the DALEX package the simplest is to use the DALEX::model_performance function. Since the explain function checks what type of task is being analyzed, it can select the appropriate performance measures for it. In our illustration, we have a binary classification model, so measures such as precision, recall, F1-score, accuracy, or AUC are calculated in the following snippet.
perf_credit = model_performance(ranger_exp)
perf_creditMeasures for: classification
recall : 0.3939394
precision : 0.5652174
f1 : 0.4642857
accuracy : 0.7272727
auc : 0.7292842
Residuals:
0% 10% 20% 30% 40% 50%
-0.89514048 -0.46432526 -0.34574427 -0.23305126 -0.15721556 -0.10991627
60% 70% 80% 90% 100%
-0.06375159 0.04445254 0.47513492 0.70182035 0.95778651 Each explanation can be drawn with the generic plot() function, for classification the default graphical summary is the ROC curve.
library("ggplot2")
old_theme = set_theme_dalex("ema")
plot(perf_credit, geom = "roc")
Various visual summaries may be selected with the geom argument. In credit risk, the LIFT curve is also a popular graphical summary.
The process of model understanding usually starts with understanding which features are the most important. A popular technique for this is the permutation feature importance introduced in Section 11.2.2. More about this technique can be found in the chapter Variable-importance Measures of the EMA book.
The DALEX::model_parts() function calculates the importance of features and its results can be visualized with the generic plot() function.
ranger_effect = model_parts(ranger_exp)
head(ranger_effect) variable mean_dropout_loss label
1 _full_model_ 0.2707158 Ranger Credit
2 amount 0.2655166 Ranger Credit
3 employment_duration 0.2721151 Ranger Credit
4 purpose 0.2744589 Ranger Credit
5 other_debtors 0.2748437 Ranger Credit
6 property 0.2781801 Ranger Creditplot(ranger_effect, show_boxplots = FALSE)
The bars start in the model loss (here 1 minus AUC) calculated for the selected data and end in a loss calculated for the data after the permutation of the selected feature. The more important the feature, the more the model will lose after its permutation. The plot indicates that status is the most important feature.
The type argument in the model_parts function allows you to specify how the importance of the features is to be calculated, by the difference of the loss functions (type = "difference"), by the quotient (type = "ratio"), or without any transformation (type = "raw").
Once we know which features are most important, we can use feature effect plots to understand its effect and show how the model, on average, changes with changes in selected features. The DALEX::model_profile() function calculates effects and, by default, it calculates the partial dependence profiles (see Section 11.2.3).
ranger_profiles = model_profile(ranger_exp)
ranger_profilesTop profiles :
_vname_ _label_ _x_ _yhat_ _ids_
1 duration Ranger Credit 4 0.2619281 0
2 duration Ranger Credit 5 0.2619281 0
3 duration Ranger Credit 6 0.2619281 0
4 duration Ranger Credit 7 0.2609917 0
5 duration Ranger Credit 8 0.2773760 0
6 duration Ranger Credit 9 0.2788940 0plot(ranger_profiles) +
theme(legend.position = "top") +
ggtitle("Partial Dependence for Ranger Credit model","")
From the figure above, we can read the interesting information that the random forest model has learned a non-monotonic relationship for the features amount and age.
The type argument of the DALEX::model_profile() function also allows marginal profiles (with type = "conditional") and accumulated local profiles (with type = "accumulated") to be calculated.
The local model analysis aims to understand how a model behaves for a single observation. In the DALEX package, functions for local analysis have names starting with the prefix predict_.
We will carry out the following examples using a sample credit application, labeled as Steve.
steve = credit_x[249,]
steve age amount credit_history duration
249 74 5129 critical account/other credits elsewhere 9
employment_duration other_debtors property purpose
249 >= 7 yrs none real estate car (new)
savings status
249 unknown/no savings account ... < 0 DMAs shown in Figure 11.2, the local analysis starts with the calculation of a model prediction. It turns out that for the selected credit application, the chances are quite high that it will not be paid.
predict(ranger_exp, steve) bad
0.5462833 A popular technique for assessing the contributions of features to model prediction is break-down (see the Break-down Plots for Additive Attributions chapter in the EMA book for more information about this method).
The function DALEX::predict_parts() function calculates the attributions of features and its results can be visualized with the generic plot() function.
ranger_attributions = predict_parts(ranger_exp, new_observation = steve)
plot(ranger_attributions) + ggtitle("Break Down for Steve")
Looking at the plots above, we can read that the biggest contributors to the final prediction were for Steve the features credit_history and status.
The order argument allows you to indicate the selected order of the features. This is a useful option when the features have some relative conditional importance (e.g. pregnancy and sex).
By far the most popular technique for local model exploration (Holzinger et al. 2022) is Shapley values introduced in Section 11.2.5. The most popular algorithm for estimating these values is the SHAP algorithm. A detailed description of the method and algorithm can be found in the chapter SHapley Additive exPlanations (SHAP) for Average Attributions of the EMA book.
The function DALEX::predict_parts() calculates SHAP attributions, you just need to set type = "shap". Its results can be visualized with a generic plot() function.
ranger_shap = predict_parts(ranger_exp, new_observation = steve,
type = "shap")
plot(ranger_shap, show_boxplots = FALSE) +
ggtitle("Shapley values for Steve", "")
The results for Break Down and SHAP methods are generally similar. Differences will emerge if there are many complex interactions in the model.
Shapley values can take a long time to compute. This process can be sped up at the expense of accuracy. The parameters B and N can be used to tune this trade-off, where N is the number of observations on which conditional expectation values are estimated (500 by default) and B is the number of random paths used to calculate Shapley values (25 by default).
In the previous section, we have introduced a global explanation – partial dependence plots. Ceteris paribus plots are the local version of that explanation. They show the response of a model when only one feature is changed while others stay unchanged. Ceteris paribus profiles are also known as individual conditional explanations (ICE). More on this technique can be found in Section 11.2.3 and the chapter Ceteris-paribus Profiles of the EMA book.
The function DALEX::predict_profile() calculates ceteris paribus profiles which can be visualized with the generic plot() function. The blue dot stands for the prediction for Steve.
ranger_ceteris = predict_profile(ranger_exp, steve)
plot(ranger_ceteris)
Ceteris paribus profiles are often plotted simultaneously for a set of observations. This helps to see potential deviations from the additivity of the model. To draw profiles for multiple observations, simply point to more rows in the data for the predict_profile function. Deviation from additivity is manifested by profiles that are not parallel.
According to Chapter 3, performance evaluation should not be conducted on the training data but on unseen data to receive unbiased estimates for the performance. Similar considerations play a role in the choice of the underlying data set used for post-hoc interpretations of a model: Should model interpretation be based on training data or test data?
This question is especially relevant when we use a performance-based interpretation method such as the permutation feature importance (PFI, see also Section 11.2.2). For example, if the PFI is computed based on the training data, the reported feature importance values may overestimate the actual importance: the model may have learned to rely on specific features during training that may not be as important for generalization to new data. This comes from the fact that the performance measured on the training data is overly optimistic and permuting a feature would lead to a huge drop in the performance. In contrast, computing the PFI on the test data provides a more accurate estimate of the feature importance as the performance measured on test data is a more reliable estimate. Thus, for the same reasons performance evaluation of a model should be performed on an independent test set (see Chapter 3), it is also recommended to apply a performance-based interpretation method such as the PFI on an independent test set.
For performance evaluation, it is generally recommended to use resampling with more than one iteration, e.g., k-fold cross-validation or (repeated) subsampling, instead of a single train-test split as performed by the hold-out method (see Chapter 3). However, for the sake of simplicity, we will use a single test data set when applying the interpretation methods to produce explanations throughout this chapter. Since produced explanations can be considered as statistical estimates, a higher number of resampling iterations can reduce the variance and result in a more reliable explanation (Molnar et al. 2021).
For prediction-based methods that do not require performance estimation such as ICE/PD (Section 11.2.3) or Shapley values (Section 11.2.5), the differences in interpretation between training and test data are less pronounced (Molnar et al. 2022).
In this chapter, we learned how to gain post-hoc insights into a model trained with mlr3 by using the most popular approaches from the field of interpretable machine learning. The methods are all model-agnostic such that they do not depend on specific model classes. We utilized three different packages: iml, counterfactuals and DALEX. iml and DALEX offer a wide range of (partly) overlapping methods, while the counterfactuals package focuses solely on counterfactual methods. We demonstrated on the rchallenge::german data set how these packages offer an in-depth analysis of a random forest model fitted with mlr3. In the following, we show some limitations of the presented methods.
If features are correlated, the insights from the interpretation methods should be treated with caution. Changing the feature values of an observation without taking the correlation with other features into account leads to unrealistic combinations of the feature values. Since such feature combinations are also unlikely part of the training data, the model will likely extrapolate in these areas (Molnar et al. 2022; Hooker and Mentch 2019). This distorts the interpretation of methods that are based on changing single feature values such as PFI, PD plots, and Shapley values. Alternative methods can help in these cases: conditional feature importance instead of PFI (Strobl et al. 2008; Watson and Wright 2021), accumulated local effect plots instead of PD plots (Apley and Zhu 2020), and the KernelSHAP method instead of Shapley values (Lundberg, Erion, and Lee 2019b).
The explanations derived from an interpretation method can also be ambiguous. A method can deliver multiple equally plausible but potentially contradicting explanations. This phenomenon is also called the Rashomon effect (Breiman 2001b). Differing hyperparameters can be one reason, for example, local surrogate models react very sensitively to changes in the underlying weighting scheme of observations. Even with fixed hyperparameters, the underlying data set or the initial seed can lead to disparate explanations (Molnar et al. 2022).
The rchallenge::german is only a low-dimensional dataset with a limited number of observations. Applying interpretation methods off-the-shelf to higher dimensional datasets is often not feasible due to the enormous computational costs. That is why in previous years more efficient methods were proposed, e.g., Shapley value computations based on kernel-based estimators (SHAP). Another challenge is that the high-dimensional IML output generated for high-dimensional datasets can overwhelm users. If the features can be meaningfully grouped, grouped versions of the methods, e.g. the grouped feature importance proposed by Au et al. (2022), can be applied.
Model explanation allows us to confront our expert knowledge related to the problem with relations learned by the model. The following tasks are based on predictions of the value of football players based on data from the FIFA game. It is a graceful example, as most people have some intuition about how a footballer’s age or skill can affect their value. The latest FIFA statistics can be downloaded from kaggle.com, but also one can use the 2020 data available in the DALEX packages (see the DALEX::fifa data set). The following exercises can be performed in both the iml and DALEX packages, we provide solutions for both.
Prepare an mlr3 regression task for the fifa data. Select only features describing the age and skills of footballers. Train a predictive model of your own choice on this task, e.g. regr.ranger, to predict the value of a footballer.
Use the permutation importance method to calculate feature importance ranking. Which feature is the most important? Do you find the results surprising?
Use the partial dependence plot/profile to draw the global behavior of the model for this feature. Is it aligned with your expectations?
Choose one of the football players, for example, a well-known striker (e.g., Robert Lewandowski) or a well-known goalkeeper (e.g., Manuel Neuer). The following tasks are worth repeating for several different choices.
For the selected footballer, calculate and plot the Shapley values. Which feature is locally the most important and has the strongest influence on the valuation of the footballer?
For the selected footballer, calculate the ceteris paribus profiles / individual conditional expectation curves to draw the local behavior of the model for this feature. Is it different from the global behavior?