[This article was first published on r-bloggers on Machine Learning in R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Machine learning models repeatedly outperform interpretable, parametric models like the linear regression model. The gains in performance have a price: The models operate as black boxes which are not interpretable.

Fortunately, there are many methods that can make machine learning models interpretable. The R package iml provides tools for analysing any black box machine learning model:

• Feature importance: Which were the most important features?
• Feature effects: How does a feature influence the prediction? (Partial dependence plots and individual conditional expectation curves)
• Explanations for single predictions: How did the feature values of a single data point affect its prediction? (LIME and Shapley value)
• Surrogate trees: Can we approximate the underlying black box model with a short decision tree?
• The iml package works for any classification and regression machine learning model: random forests, linear models, neural networks, xgboost, etc.

This blog post shows you how to use the iml package to analyse machine learning models. While the mlr package makes it super easy to train machine learning models, the iml package makes it easy to extract insights about the learned black box machine learning models.

Let’s explore the iml-toolbox for interpreting an mlr machine learning model with concrete examples!

## Data: Boston Housing

We’ll use the MASS::Boston dataset to demonstrate the abilities of the iml package. This dataset contains median house values from Boston neighbourhoods.

data("Boston", package  = "MASS")
#>      crim zn indus chas   nox    rm  age    dis rad tax ptratio  black
#> 1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90
#> 2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90
#> 3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83
#> 4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63
#> 5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90
#> 6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12
#>   lstat medv
#> 1  4.98 24.0
#> 2  9.14 21.6
#> 3  4.03 34.7
#> 4  2.94 33.4
#> 5  5.33 36.2
#> 6  5.21 28.7

## Fitting the machine learning model

First we train a randomForest to predict the Boston median housing value:

library("mlr")
data("Boston", package  = "MASS")

# create an mlr task and model
tsk = makeRegrTask(data = Boston, target = "medv")
lrn = makeLearner("regr.randomForest", ntree = 100)
mod = train(lrn, tsk)

We create a Predictor object, that holds the model and the data. The iml package uses R6 classes: New objects can be created by calling Predictor$new(). Predictor works best with mlr models (WrappedModel-class), but it is also possible to use models from other packages. library("iml") X = Boston[which(names(Boston) != "medv")] predictor = Predictor$new(mod, data = X, y = Boston$medv) ## Feature importance We can measure how important each feature was for the predictions with FeatureImp. The feature importance measure works by shuffling each feature and measuring how much the performance drops. For this regression task we choose to measure the loss in performance with the mean absolute error (‘mae’); another choice would be the mean squared error (‘mse’). Once we created a new object of FeatureImp, the importance is automatically computed. We can call the plot() function of the object or look at the results in a data.frame. imp = FeatureImp$new(predictor, loss = "mae")
plot(imp)

imp$results #> feature original.error permutation.error importance #> 1 lstat 0.9546924 4.4301626 4.640408 #> 2 rm 0.9546924 3.4838984 3.649237 #> 3 nox 0.9546924 1.6936588 1.774036 #> 4 crim 0.9546924 1.6168185 1.693549 #> 5 dis 0.9546924 1.5895507 1.664987 #> 6 ptratio 0.9546924 1.5544831 1.628255 #> 7 indus 0.9546924 1.5182089 1.590260 #> 8 age 0.9546924 1.4034536 1.470058 #> 9 tax 0.9546924 1.3047505 1.366671 #> 10 black 0.9546924 1.2429762 1.301965 #> 11 rad 0.9546924 1.0200721 1.068483 #> 12 chas 0.9546924 0.9989070 1.046313 #> 13 zn 0.9546924 0.9587879 1.004290 ## Partial dependence Besides learning which features were important, we are interested in how the features influence the predicted outcome. The Partial class implements partial dependence plots and individual conditional expectation curves. Each individual line represents the predictions (y-axis) for one data point when we change one of the features (e.g. ‘lstat’ on the x-axis). The highlighted line is the point-wise average of the individual lines and equals the partial dependence plot. The marks on the x-axis indicates the distribution of the ‘lstat’ feature, showing how relevant a region is for interpretation (little or no points mean that we should not over-interpret this region). pdp.obj = Partial$new(predictor, feature = "lstat")
plot(pdp.obj)

If we want to compute the partial dependence curves for another feature, we can simply reset the feature. Also, we can center the curves at a feature value of our choice, which makes it easier to see the trend of the curves:

pdp.obj$set.feature("rm") pdp.obj$center(min(Boston$rm)) plot(pdp.obj) ## Surrogate model Another way to make the models more interpretable is to replace the black box with a simpler model – a decision tree. We take the predictions of the black box model (in our case the random forest) and train a decision tree on the original features and the predicted outcome. The plot shows the terminal nodes of the fitted tree. The maxdepth parameter controls how deep the tree can grow and therefore how interpretable it is. tree = TreeSurrogate$new(predictor, maxdepth = 2)
plot(tree)

We can use the tree to make predictions:

lime.explain$results #> beta x.recoded effect x.original feature feature.value #> rm 4.3059097 6.575 28.311356 6.575 rm rm=6.575 #> ptratio -0.5208891 15.300 -7.969603 15.3 ptratio ptratio=15.3 #> lstat -0.4295847 4.980 -2.139332 4.98 lstat lstat=4.98 plot(lime.explain) ## Explain single predictions with game theory An alternative for explaining individual predictions is a method from coalitional game theory named Shapley value. Assume that for one data point, the feature values play a game together, in which they get the prediction as a payout. The Shapley value tells us how to fairly distribute the payout among the feature values. shapley = Shapley$new(predictor, x.interest = X[1,])
plot(shapley)

We can reuse the object to explain other data points:

shapley$explain(x.interest = X[2,]) plot(shapley) The results in data.frame form can be extracted like this: results = shapley$results
#>   feature         phi      phi.var feature.value
#> 1    crim  0.00678256  1.084280662  crim=0.02731
#> 2      zn -0.03610350  0.007022100          zn=0
#> 3   indus  0.02527449  0.267206656    indus=7.07
#> 4    chas -0.02598383  0.009941115        chas=0
#> 5     nox  0.41512115  0.928767187     nox=0.469
#> 6      rm -0.85564962 12.309151401      rm=6.421

The iml package is available on CRAN and on Github.