Py: Explainable Models with SHAP#

by Jonathan Tan. Originally published in Actuaries Digital as Explainable ML: A peek into the black box through SHAP.

With data becoming more widely available, there are more and more companies using powerful machine learning models to outperform their competitors. The insurance industry is no exception, and a good model can help the insurer gain a competitive advantage in areas such as(but not limited to):

  • Pricing accurately to reduce the risk of adverse selection

  • Reserving models that allow for optimal use of assets

  • Customer segmentation for campaigns to maximize cost and benefit

  • NLP models for customer sentiment analysis

  • Fraud detection for claims

  • Underwriting automation for life insurance

  • Customer service automation

However, the biggest downside of using these powerful models is losing the ability to easily interpret the results and the driving factors of the predicted outcomes, hence the term “black box models”. While the police department in the Minority Report certainly did not think twice about understanding their predictions before launching a full-on manhunt after their own unit head (played by Tom Cruise), it is highly unlikely that upper management and other stakeholders will react in the same fashion on the first instance of looking at a set of model results.

While an actuary may not necessarily be involved in all the areas listed above, they would most probably be exposed to at least one of them in their day job. GLMs and GAMs offer this transparency and unambiguity of model results, but limits the model choice, which may sometimes cause an issue when the underlying interactions are more complex and require a more complicated black box model to accommodate.

Thankfully, there are some tools and techniques (called explanation models) available that aims to explain the driving factors of a prediction from a powerful black box model, with some example use cases (depending on the model):

  • Drivers of unexpectedly low/high price of a policy to external stakeholders

  • Providing evidence to regulators that underwriting decisions are data driven

  • Providing the driving factors to why a claim was deemed to be fraudulent

  • Reasons a particular customer profile is preferable for a targeted campaign

  • Checking and debugging the model to make sure that the driving factors of a certain prediction are in line with expectations

This article will provide a gentle introduction to SHAP, one of the explanation models mentioned above. It was introduced as a unified approach to model interpretation and is currently the only explanation model that relies on solid and well understood economic theory instead of heuristic methods, which other popular methods (like LIME) are based on.

While this article will not provide a universal solution for the actuary to start using these black box models imprudently, it aims to encourage the exploration of these powerful black box models and give the reader a starting point for its interpretation.

Run this from the browser#

Although this article will make use of the Python implementation, there is also a R wrapper for the SHAP package called shapper, and ml3, which works with the DALEX framework. For more information on how to get started with DALEX, check out Jacky Poon’s top 10 R packages for data analysis.

The python code shown this snippet can be run from the browser through this link.

Note that some loading time is to be expected, particularly when performing grid search and cross-validation for the machine learning models. However, the binder does not require the reader to download the dataset, or any other dependencies to run the models. Furthermore, the reader is encouraged to play around with the code and look at how SHAP explains the different predictions, as well as explore the other graphs that the SHAP package has to offer.

Tools and Packages#

The example uses Python3 and the main packages that we will be using are listed below:

Note that the sklearn pipeline is the main backbone of the process here out of convenience but any other workflow/ pre-processing framework would work just as well.

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import shap
import warnings
from IPython.display import Image
from IPython.core.display import HTML 

from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.metrics import mean_squared_error

from sklearn.ensemble import GradientBoostingRegressor


Importing, Preprocessing, Modelling#

The dataset used in this article includes some policyholder information as follows:

  • age: Policyholder age (integer)

  • sex: Gender of the policyholder (string)

  • bmi: Body mass index (float)

  • children: Number of children/ dependents of policyholder (integer)

  • smoker: Smoking state of policyholder (string)

  • region: Residential area of the policyholder in the US (string)

As well as a claims cost response column. The data can be downloaded from the Github repository here.

Note that although all the usual steps in the modelling process are involved here, this article will mainly focus on using the SHAP package and interpreting the model predictions. With that in mind, the rest of this section will not be covered in detail, but includes:

  • Importing the data

  • Splitting the dataset into training and testing sets

  • Creating a “preprocessor” step to one-hot-encode categorical features

  • Building a sklearn pipeline with the pre-processing step and a GBM Regressor model

  • Performing grid search and cross-validation for the GBM model to get the “best” set of hyperparameters

  • 2 very brief/high-level checks to ensure that:

    • Model predictions approximate the actual values

    • The model has not overfitted on the training set


# Import data from GitHub repository
df = pd.read_csv("")

# Get column names for independent and response variables
rating_factors_col = list(df.columns[:-1])
claims_col = df.columns[-1]
age sex bmi children smoker region charges
0 19 female 27.90 0 yes southwest 16884.9240
1 18 male 33.77 1 no southeast 1725.5523
2 28 male 33.00 3 no southeast 4449.4620
# Splitting up the categorical columns and numerical columns
# This is to differentiate the pre-processing functions applied later in the pipeline

num_features = list(df[rating_factors_col].select_dtypes(include=["int64", "float64"]).columns)
cat_features = [col for col in rating_factors_col if col not in num_features]


# Split our dataset into training and testing sets

X_train, X_test, y_train, y_test = train_test_split(
# Create a preprocessor step
# This ignores numerical columns but applies OHE to categorical features

preprocessor = ColumnTransformer(
    transformers = [
        ("numerical", "passthrough", num_features),
        ("categorical", OneHotEncoder(sparse=False, handle_unknown="ignore"), cat_features)
# Some wrangling done to get OHE feature names
# Note that preprocess.get_feature_names() can be used as well but by default produces x prefix which I did not want

ohe_categories =["categorical"].categories_
ohe_categories_concat = [f"{col}__{val}" for col, vals in zip(cat_features, ohe_categories) for val in vals]
rating_factors_encoded = num_features + ohe_categories_concat
# Create a dictionary of hyperparameters to be optimised
# Create a sklearn Pipeline object which takes in the preprocessor step as well as a GBM object

param_grid = {
    "model__learning_rate": [0.01, 0.1],
    "model__max_depth": [2, 3, 4],
    "model__min_samples_leaf": [1, 3, 5],
    "model__min_samples_split": [2, 4],
    "model__n_estimators": [100, 200]

gbm_model = Pipeline([("preprocessor", preprocessor), ("model", GradientBoostingRegressor())])

Modelling: Gradient Boosted Machine#

# Instantiate a gridsearch-cross validation object
# This takes in our previously created gbm pipeline and grid of hyperparameters to be tested
# Using a standard 5-part CV split with MSE as a score as we are doing a regression exercise

gs = GridSearchCV(
# Run the gridsearch, set the hyperparameters to the set which gave us the lowest MSE
# Fit the optimised model on our training data set and get test set predictions

_gs =, y_train)

gbm_model.set_params(**gs.best_params_), y_train)
gbm_pred = gbm_model.predict(X_test)

Reasonableness checks#

# Function created to create a multiple bar chart comparing actual vs predicted claims by deciles

def plot_avepp(df, expected, actual):
    avepp = df.assign(
        model_bands=pd.qcut(df[expected], 10, labels=np.arange(10) + 1),
    ).groupby("model_bands")[[expected, actual]].agg("sum").reset_index()
    n = avepp["model_bands"].max()
    width = 0.35
    ax = plt.subplot(111)
    pred_fig = + width, avepp[expected], width)
    actual_fig =, avepp[actual], width)
    ax.set_title("AvE Decile Plot")
    ax.set_xticks(np.arange(n) + width / 2)
    ax.set_xticklabels(np.arange(n) + 1)
    ax.legend((pred_fig[0], actual_fig[0]), ("Pred", "Actual"))
test_results = pd.concat(
        pd.DataFrame(gbm_pred, columns=["predicted"]),
).drop("index", axis=1)

plot_avepp(test_results, "predicted", "charges")
# Checking RMSE for training and testing sets
# An indication of overfitting would be a very low training RMSE with a high testing RMSE

print("Overfitting Check")
print("Testing RMSE: {}".format(mean_squared_error(y_test, gbm_pred, squared=False)))
print("Training RMSE: {}".format(mean_squared_error(y_train, gbm_model.predict(X_train), squared=False)))
Overfitting Check
Testing RMSE: 3720.6925277868554
Training RMSE: 4309.477401843349

Model Validation#

What is SHAP?#

SHAP stands for SHapley Additive exPlanations and uses a game theory approach (Shapley Values) applied to machine learning to “fairly allocate contributions” to the model features for a given output. The underlying process of getting SHAP values for a particular feature \(f\) out of the set \(F\) can be summarized as follows:

  • Get the Power Set of \(F\), which contains \(2^{F}\) combinations of features

  • Run the model for all possible combinations, treating features outside of the subset as missing

  • Record the marginal contributions to the model output for \(f\)

  • Calculate a weighted sum of \(f\)’s marginal contributions to get the SHAP value of feature \(f\) for a given output. Or in other words, \(f\)’s expected contribution to the model output

Just from looking at the steps mentioned above, we can infer that this approach would have 2 major issues to solve:

  1. Statistical models typically are not able to handle missing values

  2. \(2^{F}\) is exponential complexity \(O(2^{N})\), which is not feasible in most cases

Thankfully, there is a user-friendly API for us to define what “missing” means by creating background dataset that imputes the missing values while running our models. As for (2), SHAP values are generally approximated with sampling subsets of \(2^{F}\). However, in the case of tree-based models and deep neural networks, there are elegant solutions implemented such that we can get exact solutions that run in polynomial \(O(N^{2})\) time!

While the article will not delve into the details on how the SHAP algorithm handles these issues, there are various resources in the appendix with more detail.

SHAP package in Python#

The SHAP python framework provides a variety of visualisations for model validation that can be found here. However, for the purposes of this article, the focus will be on the Waterfall, Force and Dependency plots to interpret model predictions.

We will also look at the results in the context of a particular observation with index=30.

idx = 30
feature = "bmi"

shap_df = pd.DataFrame(gbm_model.named_steps["preprocessor"].transform(X_train), columns=rating_factors_encoded)

Like the LIME package, SHAP works with explainer objects to calculate the results, and provides us with 3 main explainer categories:

  • shap.TreeExplainer

  • shap.DeepExplainer

  • shap.KernelExplainer

The first 2 are model specific algorithms, which makes use of the model architecture for optimizations to compute exact SHAP values as mentioned above. The KernelExplainer on the other hand, is a model agnostic algorithm uses sampling to approximate SHAP values. Since a GBM regressor is the model of choice, the TreeExplainer object is used for explanations in the following sections.

gbm_explainer = shap.TreeExplainer(gbm_model.named_steps["model"])
gbm_shap_values = gbm_explainer(shap_df)

# This line below is a quick workaround to get pass an assert condition in the SHAP 
# plots.waterfall source code. Can be ignored
gbm_shap_values.base_values = gbm_shap_values.base_values[:, 0]

Waterfall Plots (Local)#

The SHAP waterfall plots aims to explain how individual claim predictions are derived.

  1. The Y-axis encodes features, with the values observed for idx=30

  2. The X-axis encodes the range of our response (claims costs) in dollars

  3. Note that \(E[f(X)]=\$13,189.258\) is the result from the null model, or the global average claims cost from our dataset

  4. Also note that \(f(x)=\$14,959.311\) is the model prediction for values observed in (1)

The plot below shows that this policyholder’s expected claims cost is ~13% higher than the average, and both the gender and region did not contribute materially to this outcome. While the higher-than-average BMI did contribute to ~$1.3k more to the costs, this was mostly offset by the fact that she had no children.

The 2 main drivers mostly offset each other as well, but the fact that she is aged outweighed the effect of being a non-smoker by ~$1.4k, which makes up the bulk of the net increase from the average.

shap.plots.waterfall(gbm_shap_values[idx], max_display=14, show=False)

Force plots (“Global”)#

The force plots in the SHAP package can output both local and “global” interpretation graphs. While it does not provide a global explanation in the form of an equation like in our GLMs, it does give us a model-level view of the results to work with. This is done by stacking and sorting all the SHAP values for all predictions into 1 plot as shown below.

While both axes can be customized from the combo box, by default, the Y-axis shows the output value of the model, while the X-axis plots all the samples in the dataset sorted by similarity (sorting it by output value is easier to read). Hovering over an area of the graph, provides a quick summary of the significant rating factors that are driving the modelled costs up (in red) or down (in blue).

At a high level the diagram shows that the model has put a lot of emphasis on the state of smoking and age, which is in line with high BMIs have an interaction effect with smoking, accounting for many of the higher claim predictions.

shap.force_plot(gbm_explainer.expected_value[0], gbm_shap_values.values, shap_df)
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.