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
Husky vs. Wolf
A well-known anecdote in the explainability literature (Ribeiro et al., 2016).
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)
Adversarial stickers
Adversarial stickers.
Adversarial text
“TextAttack 🐙 is a Python framework for adversarial attacks, data augmentation, and model training in NLP”
Demo
LLMs are even more opaque
Run the permutation test
Decision tree
Permutation importance accuracy plot
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
Scoring rules
Example from Rudin (2019).
When to choose inherent interpretability?
Neural Additive Models (NAMs)
Each covariate (or select interactions) receives its own subnetwork, which contributes additively.
One variable’s effect on the prediction
Various plots showing the effect of one variable on the prediction; from (Molnar, 2020, sec. 12)
Useful for Hyperparameter Tuning
Partial dependence plots for hyperparameter tuning
“Explaining” to Bob
The bold red cross is the instance being explained. LIME samples instances (grey nodes), gets predictions using f (gamma MDN) and weighs them by the proximity to the instance being explained (represented here by size). The dashed line g is the learned local explanation.
“Again the approximation must be imperfect, otherwise one would throw out the black box and instead use the explanation as an inherently interpretable model.” (Rudin et al., 2022)
SHAP waterfall plot
SHAP waterfall plot
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)
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');
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
rank_main_effects trains an ensemble of shallow NAMs and scores each covariate by the variance of its shape function — a proxy for how much signal it carries — without yet committing to a cutoff.
Pick the interactions
rank_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. anam draws it as a bar chart with the real category labels, still in log-rate units.
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
Reading a swing straight off the grid is unfair: it rewards steep extrapolation into empty regions (like old agec). Instead measure the standard deviation of each shape function across the training data — the variance-based importance anam itself uses. Weighting by where the data really is, 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 from data density, not grid extrapolation");
An interaction effect
The interaction over a map
Code
import contextily as ctx
# contextily wants x=longitude, y=latitude, so put long on the x-axis.
long_grid, lat_grid, effect = model_anam.effect("long", "lat")
fig, ax = plt.subplots(figsize=(6, 5))
# Discrete bands (clipped to the 2nd-98th percentile) read far more clearly
# than a translucent continuous heatmap does; draw them at full opacity.
levels = np.linspace(np.percentile(effect, 2), np.percentile(effect, 98), 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 basemap, then recolour it to opaque black linework on a fully
# transparent background: map each tile pixel's darkness to an alpha value.
# This lays Belgium's borders and city labels over the colours without the
# pale basemap fill washing them out.
ctx.add_basemap(ax, crs="EPSG:4326", source=ctx.providers.CartoDB.Positron)
tiles = ax.images[-1]
rgb = np.asarray(tiles.get_array(), dtype=float)[..., :3]
if rgb.max() > 1.0:
rgb = rgb / 255.0
darkness = 1.0 - rgb @ np.array([0.2126, 0.7152, 0.0722])
lines = np.zeros(darkness.shape + (4,))
lines[..., 3] = np.clip((darkness - 0.07) * 3.5, 0.0, 1.0)
tiles.set_data(lines)
tiles.set_zorder(5) # draw the linework above the filled contours
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 plots
# Make partial dependence plot for ageph by hand
ageph_range = np.linspace(X_train_raw['ageph'].min(), X_train_raw['ageph'].max(), 100)
y_preds = []
X_train_pd = X_train_raw.copy()
for age in ageph_range:
X_train_pd['ageph'] = age
X_train_pd_ct = ct.transform(X_train_pd)
y_pred_batch = nn_model.predict(X_train_pd_ct, verbose=0, batch_size=X_train_pd_ct.shape[0])
y_preds.append(y_pred_batch.mean())
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
Grad-CAM
See, e.g., Keras tutorial.
Criticism II
Multiple conflicting explanations (Rudin, 2019)
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)
