<jupyter_start><jupyter_text>Explainable ML: A peek into the black box through with SHAP The insurance industry has always been more conservative when it comes to modelling, and with good reason. The financial and social impact of a wrongly predicted say, fraud detection model denying a policyholder’s entitlement to a claim is huge, compared to a wrongly predicted song recommendation pushed from your favourite music player’s machine learning algorithm. As actuaries doing technical modelling for insurance premiums, some common “workarounds” would be to stick with our favourite GLMs since it is a lot more transparent, or map/tag your GBM results to that of a GLM in some way, just to get some comfort and clarity regarding what the model is doing globally.In this article, we will introduce (at a high level) the underlying concept of SHAP, as well as work through the python implementation using a sample insurance claims dataset.  Run this from your browser The python code shown this snippet can be run from your browser through this link.Note that some loading time is to be expected, particularly when performing gridsearch 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.  Tools and Packages The example uses [Python3](https://www.python.org/downloads/release/python-385/) and the main packages that we will be using are listed below:- [Shap](https://pypi.org/project/shap/) for model validation- [Matplotlib](https://pypi.org/project/matplotlib/) for visualisation as shap graphs uses it as a backend- [Pandas](https://pypi.org/project/pandas/) and [Numpy](https://pypi.org/project/numpy/) for general data manipulation- [Sci-kit learn's](https://pypi.org/project/scikit-learn/) pipeline frameworkNote that the sklearn pipeline is the main backbone of the process here out of convenience but any other framework would work just as well.<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import shap
import warnings

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, r2_score

from sklearn.ensemble import GradientBoostingRegressor

warnings.filterwarnings('ignore')
plt.ioff()<jupyter_output><empty_output><jupyter_text>Importing, Preprocessing, Modelling The dataset used in the article includes some policyholder information as follows: - **age**: Policyholder age (integer)- **sex**: Gender of the policyholder (female=0, male=1)- **bmi**: Body mass index (float)- **children**: Number of children/ dependents of policyholder (integer)- **smoker**: Smoking state of policyholder (non-smoker=0, smoker=1)- **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](https://github.com/sharmaroshan/Insurance-Claim-Prediction). 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 results. 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 "preprocesser" step to one-hot-encode categorical features- Building a sklearn pipeline with the preprocessing and a GBM Regressor model- Performing gridsearch and cross-validation for the GBM model to get the "best" set of hyperparameters  Importing<jupyter_code>df = pd.read_csv("https://raw.githubusercontent.com/sharmaroshan/Insurance-Claim-Prediction/master/insurance.csv")

rating_factors_col = list(df.columns[:-1])
claims_col = df.columns[-1]
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]<jupyter_output><empty_output><jupyter_text>Preprocessing<jupyter_code>X_train, X_test, y_train, y_test = train_test_split(
            df[rating_factors_col],
            df[claims_col],
            test_size=0.2,
            random_state=123
        )
preprocesser = ColumnTransformer(
    transformers = [
        ("numerical", "passthrough", num_features),
        ("categorical", OneHotEncoder(sparse=False, handle_unknown="ignore"), cat_features)
    ]
)
ohe_categories = preprocesser.fit(X_train).named_transformers_["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
param_grid = {
    "model__learning_rate": [0.01, 0.05, 0.1],
    "model__max_depth": [2, 3, 4, 5],
    "model__max_features": [2, 3, 5, 7],
    "model__min_samples_leaf": [3, 4, 5],
    "model__min_samples_split": [5, 8, 10, 12],
    "model__n_estimators": [100, 200, 300, 500]
}

gbm_model = Pipeline([("preprocesser", preprocesser), ("model", GradientBoostingRegressor())])<jupyter_output><empty_output><jupyter_text>Modelling: Gradient Boosted Machine<jupyter_code>gs = GridSearchCV(
    gbm_model, 
    param_grid=param_grid, 
    n_jobs=-1, 
    cv=5, 
    scoring="neg_root_mean_squared_error"
)
gs.fit(X_train, y_train)
gbm_model.set_params(**gs.best_params_)
gbm_model.fit(X_train, y_train)
gbm_pred = gbm_model.predict(X_test)<jupyter_output><empty_output><jupyter_text>Model Validation  What is SHAP? SHAP stands for Shapley Additive Explanations and uses a game theory approach (Shapley Values) applied on 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](https://en.wikipedia.org/wiki/Power_set) of $F$, which  contains $2^{F}$ combinations of features- Run the model for all combinations in the set- Record the marginal contributions to the model output for $f$- Calculated 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 contribution to the model outputFor more details, there are various resources in the appendix for the reader to delve deeper into this concept, including the original white paper.  SHAP package in Python The SHAP python framework provides a variety of visualisations for model validation and can be found [here](https://shap.readthedocs.io/en/latest/api.html). However, for our purposes, we will only be concentrating on the Partial Dependency Plot (PDP) and the SHAP Waterfall Plot.We will also look at the results in the context of a particular observation with index=30.<jupyter_code>idx = 30
feature = "bmi"

shap_df = pd.DataFrame(lr_model.named_steps["preprocesser"].transform(X_train), columns=rating_factors_encoded)<jupyter_output><empty_output><jupyter_text>Similar to the LIME package, SHAP works with explainer objects to calculate the results, and because we have opted to use a GBM regressor for this example, we will be using the TreeExplainer object.<jupyter_code>gbm_explainer = shap.TreeExplainer(gbm_model.named_steps["model"])
gbm_shap_values = gbm_explainer(shap_lr_df)

# This line below is a quick workaround to get pass an assert condition in the SHAP 
# waterfall function. Can be ignored
gbm_shap_values.base_values = gbm_shap_values.base_values[:, 0]<jupyter_output><empty_output><jupyter_text>Waterfall Plots (Local) The SHAP waterfall plots aims to explain how individual claim predictions come about.1. On the Y-axis, we have our encoded features, with the values observed for idx=302. On the X-axis, we have the range of our response (claims costs) in dollars3. Note that $E[f(X)]=\$13,189.258$ is the result from the null model, or the global average claims cost from our dataset4. Also note that $f(x)=\$14,959.311$ is the model prediction for values observed in (1)So, we can see from the plot below 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.<jupyter_code>shap.plots.waterfall(gbm_shap_values[idx], max_display=14, show=False)
plt.show()<jupyter_output><empty_output><jupyter_text>Force plots ("Global") The force plots in the SHAP package can output both local and "global" interpretation graphs, but in this example, we have inputted the whole array of SHAP values and data points into the function, to get a model-level view of the results. 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 (I personally like sorting it by output value). By hovering over an area of the graph, we can get a quick summary of the significant rating factors that are driving the modelled costs up (in red) or down (in blue). Just at a high level we can observe 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.<jupyter_code>shap.initjs()
shap.force_plot(gbm_explainer.expected_value[0], gbm_shap_values.values, shap_df)<jupyter_output><empty_output><jupyter_text>SHAP Dependency plots ("Global") The SHAP dependency plot is a very simple graph that shows how the SHAP contributions differ for different values of the feature (BMI in this case). This is like a PDP, except that PDPs show average effects while the SHAP dependency plots show the variance on the y-axis. Here, we can clearly see an interaction effect between BMI values and the smoking state of the policyholder.We see that the pink data points (non-smokers) have a more gradual slope for BMI contributions in contrast to the blue points (smokers), where the SHAP contributions jump significantly at the BMI=~30 point, which is in line with what we could guess from the previous force plot output above.<jupyter_code>shap.dependence_plot(feature, gbm_shap_values.values, shap_df)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07147504414908448, 'fraction_numerical': 0.009294544102611766, 'mean_word_length': 5.351829988193624, 'pattern_counts': {'":': 6, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 10, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '98778', 'license': 'MIT', 'n_tokens_mistral': 3078, 'n_tokens_neox': 2916, 'n_words': 1324, 'path': '.ipynb_checkpoints/shap_article_nb_v2-checkpoint.ipynb', 'repo_name': 'jtsw1990/shap_values', 'size': '10759'}