import random
import matplotlib.pyplot as plt
import numpy as np
import numpy.random as rnd
import pandas as pd
import geopandas as gpd
import contextily as ctx
from shapely.geometry import Point
import plotly.graph_objects as go
import keras
from keras.metrics import SparseTopKCategoricalAccuracy
from keras.models import Sequential, Model
from keras.layers import Dense, Input
from keras.callbacks import EarlyStopping
from sklearn import tree, set_config
from sklearn.preprocessing import LabelEncoder, OrdinalEncoder, StandardScaler, OneHotEncoder
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.model_selection import train_test_split
from sklearn.compose import make_column_transformer
from sklearn.metrics import mean_poisson_deviance, mean_absolute_error
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.inspection import PartialDependenceDisplay
set_config(transform_output="pandas")
import statsmodels.api as sm
import statsmodels.formula.api as smf
from pygam import PoissonGAM, s, f
from anam import (
ANAM, train_anam, TrainingConfig, prepare_data,
poisson_nll, plot_shape_function,
)Interpretability
ACTL3143 & ACTL5111 Deep Learning for Actuaries
The right to an explanation
On your thirtieth birthday, you awake to find two police officers standing over you. They inform you that you’re under arrest. When you ask what for, they tell you that they don’t know, but their boss, the inspector, will be along soon, and he’ll be able to tell you. A few hours later, you’re ushered into a neighboring tenant’s bedroom. The inspector tells you that he doesn’t know what you’ve been charged with, but you’re free to go about your life, except that you must report to a hearing on Sunday. And so you move through the criminal justice system, never knowing what crime you’ve allegedly committed, nor understanding the rules of the process. On the eve of your thirty-first birthday, you’re taken away and executed.
Fortunately, this story is not about you, but is instead about Josef K., the fictional protagonist of Kafka’s The Trial.
— Vredenburgh (2022, p. 1)
Unfortunately, a broadly similar story is probably true about you, thanks to the use of algorithms to aid decision-making in a variety of institutions. One such example comes from Washington, DC. In 2009, the city set out to improve its school system, where only 8 percent of eighth graders were performing at grade level in math. It introduced IMPACT, a model to evaluate teacher performance developed by Mathematica Policy Research. The aim of the model was to identify poorly performing teachers, as input to firing decisions. The model that Mathematica Policy Research developed was complex and proprietary, and teachers couldn’t find anyone to explain to them why they’d been fired. In response to the firing of her colleagues on the basis of the model’s outputs, Sarah Bax, a math teacher, asked: “How do you justify evaluating people by a measure for which you are unable to provide an explanation?”.
— Vredenburgh (2022, p. 1)
Husky vs. Wolf
A well-known anecdote in the explainability literature (Ribeiro et al., 2016).
Run the permutation test
Decision tree
Permutation importance accuracy plot
Adversarial examples
A demonstration of fast adversarial example generation applied to GoogLeNet on ImageNet. By adding an imperceptibly small vector whose elements are equal to the sign of the elements of the gradient of the cost function with respect to the input, we can change GoogLeNet’s classification of the image (Goodfellow et al., 2014)
What is inherent interpretability?
“Interpretability by design is decided on the level of the machine learning algorithm. If you want a machine learning algorithm that produces interpretable models, the algorithm has to constrain the search of models to those that are interpretable. The simplest example is linear regression: When you use ordinary least squares to fit/train a linear regression model, you are using an algorithm that will produce … models that are linear in the input features. Models that are interpretable by design are also called intrinsically or inherently interpretable models.”
— Molnar (2020)
Some processes don’t fit the tree structure
Neural Additive Models (NAMs)
A NAM is a GAM where each shape function f_j is a small neural network, learned jointly with the rest of the model (Agarwal et al., 2021).
Each covariate (or select interactions) receives its own subnetwork, which contributes additively.
Explaining a black box
A black box machine learning model is a formula that is either too complicated for any human to understand, or proprietary, so that one cannot understand its inner workings. Black box models are difficult to troubleshoot, which is particularly problematic for medical data. Black box models often predict the right answer for the wrong reason (the “Clever Hans” phenomenon), leading to excellent performance in training but poor performance in practice. There are numerous other issues with black box models. In criminal justice, individuals may have been subjected to years of extra prison time due to typographical errors in black box model inputs and poorly-designed proprietary models for air quality have had serious consequences for public safety during wildfires…
— Rudin et al. (2022, pp. 2–3)
Determining whether a black box model is fair with respect to gender or racial groups is much more difficult than determining whether an interpretable model has such a bias… Explaining black boxes, rather than replacing them with interpretable models, can make the problem worse by providing misleading or false characterizations, or adding unnecessary authority to the black box. There is a clear need for innovative machine learning models that are inherently interpretable.
— Rudin et al. (2022, pp. 2–3)
Example: ceteris paribus plot
Fit a random forest to predict diabetes progression, take one patient, and vary their bmi while holding every other feature fixed.
Code
from sklearn.datasets import load_diabetes
from sklearn.ensemble import RandomForestRegressor
diab_X, diab_y = load_diabetes(return_X_y=True, as_frame=True, scaled=False)
diab_model = RandomForestRegressor(random_state=42).fit(diab_X, diab_y)
patient = diab_X.iloc[[0]]
bmi_grid = np.linspace(diab_X["bmi"].min(), diab_X["bmi"].max(), 100)
cp_curve = [diab_model.predict(patient.assign(bmi=bmi))[0] for bmi in bmi_grid]
plt.plot(bmi_grid, cp_curve)
plt.scatter(patient["bmi"], diab_model.predict(patient), color="black", zorder=3)
plt.xlabel("bmi")
plt.ylabel("Predicted progression");
Example: ICE plot
Draw the same ceteris paribus curve for every observation in the dataset.
Example: partial dependence plot
Averaging the ICE curves gives the partial dependence plot.
Warning: correlated features
These methods (and permutation importance) change one feature at a time, leaving the others untouched — so they can manufacture inputs that could never exist.
Explain the model on the raw features
The model to explain is the composition f \circ t; perturb the raw features \boldsymbol{x} and re-apply the preprocessing t, rather than perturbing t(\boldsymbol{x}).
Exposure & frequency
Map of claims
Code
# Compute average claims per location
avg = train_raw.groupby(['lat', 'long'], as_index=False)['nclaims'].mean()
# Build GeoDataFrame in Web Mercator
gdf = gpd.GeoDataFrame(
avg,
geometry=[Point(xy) for xy in zip(avg.long, avg.lat)],
crs="EPSG:4326"
).to_crs(epsg=3857)
# Plot
ax = gdf.plot(
column='nclaims',
cmap='OrRd',
markersize=8,
edgecolor='grey',
linewidth=0.5,
alpha=0.8,
legend=True,
figsize=(10, 5)
)
ctx.add_basemap(ax, source=ctx.providers.CartoDB.Positron(r="@2x"))
ax.set_axis_off()
plt.title("Average Claims per Location")
plt.show()
GLM Ceteris Paribus Plots
Code
plt.plot(ageph_range, y_pred_batch, 'r')
real_age = X_train_raw['ageph'].iloc[0:1].item()
orig_prediction = glm_model.predict(X_train_raw.iloc[0:1]).item()
plt.plot(real_age, orig_prediction, 'o', color='r', markersize=5)
plt.xlabel('Age of Policyholder')
plt.ylabel('Predicted Number of Claims')
plt.title('GLM Predictions for Varying Age of Policyholder #1');
Zoomed out
We can see the trend if we zoom out to ages between 0 & 1,000.
Code
ageph_range = np.linspace(0, 1000, 20)
X_train_first = X_train_raw.iloc[0:1].copy()
y_pred_batch = []
for ageph in ageph_range:
X_train_first['ageph'] = ageph
y_pred_ageph = glm_model.predict(X_train_first)
y_pred_batch.append(y_pred_ageph)
plt.plot(ageph_range, y_pred_batch, 'r')
plt.plot(real_age, orig_prediction, 'o', markersize=5, color='r')
plt.xlabel('Age of Policyholder')
plt.ylabel('Predicted Number of Claims')
plt.title('GLM Predictions for Varying Age of Policyholder #1');
Interrogating an intentional bug
We fed expo in as an ordinary covariate — the CP plot exposes the consequence.
Code
expo_range = np.linspace(train_raw['expo'].min(), train_raw['expo'].max(), 20)
X_train_first = X_train_raw.iloc[0:1].copy()
y_pred_batch = []
for expo in expo_range:
X_train_first['expo'] = expo
y_pred_batch.append(glm_model.predict(X_train_first).item())
plt.plot(expo_range, y_pred_batch, 'r', label='GLM ceteris paribus')
plt.plot([0, expo_range[-1]], [0, y_pred_batch[-1]], '--', color='grey',
label='Proportional to exposure')
real_expo = X_train_raw['expo'].iloc[0:1].item()
plt.plot(real_expo, orig_prediction, 'o', markersize=5, color='r')
plt.xlabel('Exposure')
plt.ylabel('Predicted Number of Claims')
plt.title('GLM Predictions for Varying Exposure of Policyholder #1')
plt.legend();
What if we look at multiple policyholders?
Code
ageph_range = np.linspace(train_raw['ageph'].min(), train_raw['ageph'].max(), 20)
# Just overlay the plots for the first 10 training observations
for i in range(5):
X_train_first = X_train_raw.iloc[i:i+1].copy()
y_pred_batch = []
real_age = X_train_first['ageph'].item()
orig_prediction = glm_model.predict(X_train_first).item()
for ageph in ageph_range:
X_train_first['ageph'] = ageph
y_pred_ageph = glm_model.predict(X_train_first)
y_pred_batch.append(y_pred_ageph)
res = plt.plot(ageph_range, y_pred_batch, label=f'Policyholder {i+1}')
plt.plot(real_age, orig_prediction, 'o', markersize=5, color=res[0].get_color())
plt.xlabel('Age of Policyholder')
plt.ylabel('Predicted Number of Claims')
plt.title('GLM Predictions for Varying Age of Policyholders')
plt.legend(loc='upper left', bbox_to_anchor=(1, 1), ncol=1);
Logarithmic scale
Code
ageph_range = np.linspace(train_raw['ageph'].min(), train_raw['ageph'].max(), 100)
# Just overlay the plots for the first 10 training observations
for i in range(5):
X_train_first = X_train_raw.iloc[i:i+1].copy()
log_y_pred_batch = []
real_age = X_train_first['ageph'].item()
orig_prediction = np.log(glm_model.predict(X_train_first).item())
for ageph in ageph_range:
X_train_first['ageph'] = ageph
y_pred_ageph = glm_model.predict(X_train_first)
log_y_pred_batch.append(np.log(y_pred_ageph))
res = plt.plot(ageph_range, log_y_pred_batch, label=f'Policyholder {i+1}')
plt.plot(real_age, orig_prediction, 'o', markersize=5, color=res[0].get_color())
plt.xlabel('Age of Policyholder')
plt.ylabel('Log-Predicted Num. Claims')
plt.title('GLM Log-Predictions for Varying Age of Policyholders')
plt.legend(loc='upper left', bbox_to_anchor=(1, 1), ncol=1);
Ceteris paribus plot for ageph
Code
# Define range of ageph values
ageph_range = np.linspace(X_train_gam['ageph'].min(), X_train_gam['ageph'].max(), 100)
# Start with a copy of the first encoded training row
X_first = X_train_gam.iloc[0:1].copy()
# Set the dtype of 'ageph' to float
X_first['ageph'] = X_first['ageph'].astype(float)
gam_preds = []
real_age = X_first['ageph'].item()
orig_prediction = gam_model.predict(X_first).item()
for ageph in ageph_range:
X_first['ageph'] = ageph
y_pred = gam_model.predict(X_first)[0]
gam_preds.append(y_pred)
plt.plot(ageph_range, gam_preds, label='GAM', color='red')
# Add original GLM prediction for reference
plt.plot(real_age, orig_prediction, 'o', color='red', markersize=5)
plt.xlabel('Age of Policyholder')
plt.ylabel('Predicted Number of Claims')
plt.title('GAM Predictions for Varying Age')
plt.legend();
Pick the main effects
Trains an ensemble of shallow NAMs and scores each covariate by the variance of its shape function.
Pick the interactions
Scores every pairwise term that passes the strong heredity rule (a pair is only considered if both of its main effects are in the top n\_main), measuring how much each candidate improves the validation loss over the mains-only baseline.
Reading the shape functions
Each curve is that feature’s contribution to \log \mu; back on the original scale, \exp(f_j) is a multiplier on the expected claim count.
Code
fig = model_anam.plot_effects(["agec", "lat", "long"])
# Flag agec's sparse, extrapolated tail: 99% of cars are 18 years or younger.
agec_ax = fig.axes[0]
cutoff = anam_train["agec"].quantile(0.99)
agec_ax.axvspan(cutoff, anam_train["agec"].max(), color="#FF8FA9", alpha=0.25)
agec_ax.text(cutoff, agec_ax.get_ylim()[1], " sparse\n extrapolation",
ha="left", va="top", fontsize=7, color="#B03A5B");
A categorical effect: fuel
Continuous features give curves; a categorical feature gives one value per level.
Code
ax = model_anam.plot_effect("fuel")
# Label each bar with its multiplier on the expected claim frequency.
labels, values = model_anam.effect("fuel")
for x, v in enumerate(values):
ax.annotate(f"×{np.exp(v):.2f}", (x, v), ha="center",
va="bottom" if v >= 0 else "top",
xytext=(0, 3 if v >= 0 else -3), textcoords="offset points")
# Expand the bottom margin so the ×0.74 label on the negative bar isn't clipped.
ymin, ymax = ax.get_ylim()
ax.set_ylim(ymin - 0.1 * (ymax - ymin), ymax)
ax.set_title("Effect of fuel type (categorical)");
The same curves, one honest scale
Each panel above had its own y-axis. Put the six continuous effects on a shared scale and the picture changes: agec and power swing hard, while ageph, lat and long flatten to near-zero lines.
Code
feats = ["ageph", "bm", "power", "agec", "lat", "long"]
fig, axes = plt.subplots(2, 3, figsize=(9, 5.5), sharey=True)
fig.subplots_adjust(hspace=0.45)
for ax, feat in zip(axes.ravel(), feats):
grid, values = model_anam.effect(feat)
ax.plot(grid, values)
ax.axhline(0, color="0.7", lw=0.6, zorder=0)
# Decile rug: black ticks at the 10th–90th percentiles of the data.
deciles = model_anam.deciles_[feat]
ax.plot(deciles, np.zeros_like(deciles), "|",
transform=ax.get_xaxis_transform(), color="black",
markeredgewidth=1.2, markersize=12)
ax.set_title(feat)
fig.supylabel("Partial effect (log-rate)");
Which effects actually matter
Measure the standard deviation of each shape function across the training data. Exponentiating gives a typical multiplier on the expected claim frequency.
Code
from matplotlib.patches import Patch
importance = {}
for feat in model_anam.features_:
grid, values = model_anam.effect(feat)
grid = np.asarray(grid)
if np.issubdtype(grid.dtype, np.number):
effect_on_data = np.interp(anam_train[feat], grid.astype(float), values)
else:
lookup = {str(k): v for k, v in zip(grid, values)}
effect_on_data = anam_train[feat].astype(str).map(lookup).to_numpy()
importance[feat] = float(np.std(effect_on_data))
names = sorted(importance, key=importance.get)
fig, ax = plt.subplots(figsize=(7, 4))
# Colour splits the material effects from the ones that do essentially nothing.
colors = ["#FF8FA9" if importance[f] < 0.002 else "#3F9999" for f in names]
bars = ax.barh(names, [importance[f] for f in names], color=colors)
ax.set_xscale("log")
ax.set_xlabel("Importance: std of shape function over the data (log-rate)")
for f, bar in zip(names, bars):
ax.text(bar.get_width() * 1.15, bar.get_y() + bar.get_height() / 2,
f"×{np.exp(importance[f]):.2f}", va="center")
ax.legend(handles=[
Patch(color="#3F9999", label="material effect"),
Patch(color="#FF8FA9", label="negligible (≈ ×1.00)"),
], loc="lower right", frameon=False)
ax.set_title("Importance");
An interaction effect
Code
ax = model_anam.plot_effect("long", "lat")
ax.scatter(anam_train["long"], anam_train["lat"],
s=2, c="white", alpha=0.15, linewidths=0)
# A degree of longitude is shorter than a degree of latitude this far from
# the equator, so use the true geographic proportions (Belgium is wider
# than it is tall) rather than the default aspect ratio.
ax.set_aspect(1 / np.cos(np.deg2rad(np.mean(ax.get_ylim()))));
The interaction over a map
Code
import contextily as ctx
def padded_effect(model, name, other, pad=0.05, n_points=50):
"""
Like ``model.effect(name, other)``, but evaluated on a grid extended a
fraction ``pad`` beyond the training data's range, so the surface
extrapolates slightly past the outermost municipalities and covers the
whole country.
"""
def grid_for(feat):
lo, hi = model.feature_configs_[feat]["bounds"]
margin = pad * (hi - lo)
g_std = torch.linspace(float(lo) - margin, float(hi) + margin, n_points)
g_orig = model.scalers_[feat].inverse_transform(
g_std.numpy().reshape(-1, 1)).ravel()
return g_std, g_orig
grid_a_std, grid_a = grid_for(name)
grid_b_std, grid_b = grid_for(other)
model.module_.eval()
with torch.no_grad():
values = model.module_.get_shape_function(
name, grid_a_std, other, grid_b_std).cpu().numpy()
if f"{name}_{other}" not in model.module_.interaction_subnets:
values = values.T
return grid_a, grid_b, values
# contextily wants x=longitude, y=latitude, so put long on the x-axis.
long_grid, lat_grid, effect = padded_effect(model_anam, "long", "lat")
fig, ax = plt.subplots(figsize=(6, 5))
# Discrete bands read far more clearly than a translucent continuous heatmap
# does; draw them at full opacity. Match the colour limits to the previous
# slide by using the effect's range over the data region --- the padded
# corners beyond it spill into the end bands (via extend="both").
_, _, effect_data = model_anam.effect("long", "lat")
levels = np.linspace(effect_data.min(), effect_data.max(), 8)
cf = ax.contourf(long_grid, lat_grid, effect.T, levels=levels,
cmap="viridis", extend="both")
fig.colorbar(cf, ax=ax, label="Partial effect on log-rate", extendrect=True)
# Draw the country borders and coastline as vector lines over the colours.
# The file is Natural Earth's (public domain) 10m admin-0 boundary lines
# plus coastline, clipped to this corner of Europe.
borders = gpd.read_file("data/country-borders.geojson")
xlim, ylim = ax.get_xlim(), ax.get_ylim()
borders.plot(ax=ax, color="black", linewidth=1.0)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
# Add the place names from CARTO's labels-only tile layer: transparent
# background with white-haloed names, made for overlaying on data like this.
# The "@2x" retina tiles keep the same zoom level (so the same city labels)
# but at double the resolution, so the text stays sharp at print DPI.
ctx.add_basemap(ax, crs="EPSG:4326",
source=ctx.providers.CartoDB.PositronOnlyLabels(r="@2x"))
# A degree of longitude is shorter than a degree of latitude this far from
# the equator, so use the true geographic proportions.
ax.set_aspect(1 / np.cos(np.deg2rad(np.mean(ax.get_ylim()))))
ax.set(xlabel="long", ylabel="lat",
title="lat × long interaction over Belgium");
Build a neural network
ct = make_column_transformer(
(StandardScaler(), num_vars),
(OneHotEncoder(drop="first", sparse_output=False), cat_vars),
verbose_feature_names_out=False
)
ct.fit(X_train_raw)
X_train_nn = ct.transform(X_train_raw)
X_test_nn = ct.transform(X_test_raw)
y_train_nn = y_train_raw.values
y_test_nn = y_test_raw.values
random.seed(123)
nn_model = Sequential(
[
Input(shape=(X_train_nn.shape[1],)),
Dense(128, activation="relu"),
Dense(128, activation="relu"),
Dense(128, activation="relu"),
Dense(1, activation="exponential")
]
)
Permutation importance example
Partial dependence plot
Using scikit-learn’s PDP
(First need to trick scikit-learn a little..)
Code
# Wrap the column transformer + Keras model as a scikit-learn regressor so
# PartialDependenceDisplay can use it.
class SklearnModel(RegressorMixin, BaseEstimator):
def __init__(self, ct=None, model=None, batch_size=8192):
self.ct = ct
self.model = model
self.batch_size = batch_size
self.n_features_in_ = ct.n_features_in_
def predict(self, X):
X_trans = self.ct.transform(X)
preds = self.model.predict(X_trans, verbose=0, batch_size=self.batch_size)
return preds.ravel()
def fit(self, X=None, y=None):
return self
nn_model_skl = SklearnModel(ct, nn_model)
Explaining the ages and exposure
Explaining the bonus-malus and location
Bivariate age effects
SHAP waterfall plot
SHAP waterfall plot
Grad-CAM
See, e.g., Keras tutorial.
Conclusion
“Figure 20.2: Explainability techniques allow strengthening the feedback extracted from a model. A, data and domain knowledge allow building the model. B, predictions are obtained from the model. C, by analyzing the predictions, we learn more about the model. D, better understanding of the model allows better understanding of the data and, sometimes, broadens domain knowledge.” Biecek & Burzykowski (2021)
