[8]{.chapter-number}  [Explainable AI]{.chapter-title}

8 Explainable AI

The goal of explainable AI (xAI, aka interpretable machine learning) is to explain why a fitted machine learning model makes certain predictions. A typical example is to understand how important different variables are for predictions. The incentives for doing so range from a better technical understanding of the models, through understanding which data is important for improving predictions, to questions of fairness and discrimination (e.g. to understand if an algorithm uses skin color to make a decision).

8.1 A Practical Example

We will fit a random forest and use the iml package for xAI, see https://christophm.github.io/interpretable-ml-book/.

library(iml)
library(ranger) # different random Forest package!
library(EcoData)
library(cito)
set.seed(123)

load("sdm.RData")

X = scale(sdm_env$train[,2:20])
Y = Y_train[,4] # we analyze the fourth species
df = data.frame(Presence = Y, X)

Meaning of the bioclim variables:

Bioclim variable	Meaning
bio1	Annual Mean Temperature
bio2	Mean Diurnal Range (Mean of monthly (max temp - min temp))
bio3	Isothermality (BIO2/BIO7) (×100)
bio4	Temperature Seasonality (standard deviation ×100)
bio5	Max Temperature of Warmest Month
bio6	Min Temperature of Coldest Month
bio7	Temperature Annual Range (BIO5-BIO6)
bio8	Mean Temperature of Wettest Quarter
bio9	Mean Temperature of Driest Quarter
bio10	Mean Temperature of Warmest Quarter
bio11	Mean Temperature of Coldest Quarter
bio12	Annual Precipitation
bio13	Precipitation of Wettest Month
bio14	Precipitation of Driest Month
bio15	Precipitation Seasonality (Coefficient of Variation)
bio16	Precipitation of Wettest Quarter
bio17	Precipitation of Driest Quarter
bio18	Precipitation of Warmest Quarter
bio19	Precipitation of Coldest Quarter

rf = ranger(as.factor(Presence) ~ ., data = df, probability = TRUE)

The cito package has quite extensive xAI functionalities. However, ranger, like most other machine learning packages, has no extensive xAI functionalities. Thus, to do xAI with ranger, we have to use a generic xAI package that can handle almost all machine learning models.

When we want to use such a generic package, we first have to create a predictor object that holds the model and the data. The iml package uses R6 classes, which means new objects can be created by calling Predictor$new(). (Do not worry if you do not know what R6 classes are, just use the command.)

To make the xAI tools available to many different packages/algorithms, the iml package expects the ML-algorithm-specific predict method to be wrapped in a generic predict function in the form of function(model, newdata) predict(model, newdata), and the function wrapper should return a vector of predictions:

predict_wrapper = function(model, newdata) predict(model, data=newdata)$predictions[,2]

predictor = Predictor$new(rf, data = df[,-1], y = df[,1], predict.function = predict_wrapper)
predictor$task = "classif" # set task to classification
# "Predictor" is an object generator.

8.2 Feature/Permutation Importance

Feature importance should not be confused with random forest variable importance, although they are related. It tells us how important each variable is for prediction, can be computed for all machine learning models, and is based on a permutation approach (see the book):

imp = FeatureImp$new(predictor, loss = "ce")
plot(imp)

bio9 (Mean Temperature of Driest Quarter) is the most important variable.

8.3 Partial Dependencies

Partial dependencies are similar to allEffects plots for normal regressions. The idea is to visualize “marginal effects” of predictors (with the “feature” argument we specify the variable we want to visualize):

eff = FeatureEffect$new(predictor, feature = "Bio14", method = "pdp",
                        grid.size = 10)
plot(eff)

One disadvantage of partial dependencies is that they are sensitive to correlated predictors. Accumulated local effects can be used to account for correlation between predictors.

8.4 Accumulated Local Effects

Accumulated local effects (ALE) are basically partial dependence plots but try to correct for correlations between predictors.

ale = FeatureEffect$new(predictor, feature = "Bio14", method = "ale")
ale$plot()

If there is no collinearity, you shouldn’t see much difference between partial dependencies and ALE plots.

8.5 Friedman’s H-statistic

The H-statistic can be used to find interactions between predictors. However, again, keep in mind that the H-statistic is sensitive to correlation between predictors:

interact = Interaction$new(predictor, "Bio14",grid.size = 5L)
plot(interact)

8.6 Global Explainer - Simplifying the Machine Learning Model

Another idea is to simplify the machine learning model with a simpler model such as a decision tree. We create predictions with the machine learning model for a lot of different input values and then we fit a decision tree on these predictions. We can then interpret the simpler model.

library(partykit)

tree = TreeSurrogate$new(predictor, maxdepth = 2)
plot(tree$tree)

8.7 Local Explainer - LIME Explaining Single Instances (observations)

The global approach is to simplify the entire machine-learning black-box model via a simpler model, which is then interpretable.

However, sometimes we are only interested in understanding how single predictions are generated. The LIME (Local Interpretable Model-agnostic Explanations) approach explores the feature space around one observation and, based on this, fits a simpler model locally (e.g. a linear model):

lime.explain = LocalModel$new(predictor, x.interest = df[1,-1])
lime.explain$results
plot(lime.explain)

8.8 Local Explainer - Shapley

The Shapley method computes the so-called Shapley value, feature contributions for single predictions, and is based on an approach from cooperative game theory. The idea is that each feature value of the instance is a “player” in a game, where the prediction is the reward. The Shapley value tells us how to fairly distribute the reward among the features.

shapley = Shapley$new(predictor, x.interest = df[1,-1])
shapley$plot()

8.9 Uncertainties - the bootstrap

Standard xAI methods do not provide reliable uncertainties on the fitted curves. If you want uncertainties or p-values, the most common method is the bootstrap.

In a bootstrap, instead of splitting up the data into test / validation, we sample from the data with replacement and fit the models repeatedly. The idea is to get an estimate of the variability we would expect if we created another dataset of the same size.

k = 10 # bootstrap samples
n = nrow(df)
error = rep(NA, k)

for(i in 1:k){
  bootSample = sample.int(n, n, replace = TRUE)
  rf = ranger(as.factor(Presence) ~ ., data = df[bootSample,], probability = TRUE)
  error[i] = rf$prediction.error
}

hist(error, main = "uncertainty of in-sample error")

Note that the distinction between bootstrap and validation / cross-validation is as follows:

Validation / cross-validation estimates out-of-sample predictive error
Bootstrap estimates uncertainty / confidence interval on all model outputs (could be prediction and inference).

8.10 Exercises

8.10.0.1 xAI in cito

Dataset can be downloaded from the submission server or from dropbox

Data preparation

library(iml)
library(cito)
library(EcoData)

load("sdm.RData")
X = scale(sdm_env$train[,2:20])
Y = Y_train[,4] # we analyze the fourth species
df = data.frame(Presence = as.factor(Y), X)
# we will subsample data (absences) to reduce runtime
data_sub = df[sample.int(nrow(df), 500),]

cito includes several xAI methods directly out of the box

model = dnn(Presence~., data = data_sub, batchsize = 200L,loss = "binomial", verbose = FALSE, lr = 0.15, epochs = 300)

Try the following commands:

summary(model, n_permute = 10)
PDP(model)
ALE(model)

Moreover, try to refit the model with the option bootstrap = 5. This may take a short while. Observe how the xAI outputs change.

Also, try out other species!

model = dnn(Presence~., data = data_sub, batchsize = 200L, bootstrap = 5L, loss = "binomial", verbose = FALSE, lr = 0.15, epochs = 300)

summary(model, n_permute = 10L)

PDP(model)
ALE(model)

8.10.0.2 Interpret your models

Use the iml package to interpret the models you built for the challenges, and calculate:

Feature importance
Accumulated local effect plots
Partial dependency plots

Here is a predict_wrapper for the ranger package

predict_wrapper = function(model, newdata) predict(model, data=newdata)$predictions[,2]