midr: French Motor Data Example

Introduction

In modern actuarial science, there is an inherent tension between predictive accuracy and model transparency. While ensemble tree-based models like Gradient Boosting Machines (GBMs) frequently outperform traditional Generalized Linear Models (GLMs), their “black-box” nature presents significant hurdles for model governance, regulatory compliance, and price filing.

This notebook demonstrates a solution using Maximum Interpretation Decomposition (MID) via the {midr} package in R.

WarningCompatibility Notice

This article relies on features introduced in midr (>= 0.6.0). Please ensure your packages are up to date.

What is MID?

MID is a functional decomposition method that deconstructs a black-box prediction function \(f(\mathbf{X})\) into several interpretable components: an intercept \(g_\emptyset\), main effects \(g_j(X_j)\), and second-order interactions \(g_{jk}(X_j, X_k)\), minimizing the squared residuals \(\mathbf{E}\left[g_D(\mathbf{X})^2\right]\):

\[ f(\mathbf{X}) = g_\emptyset + \sum_{j} g_j(X_{j}) + \sum_{j < k} g_{jk}(X_{j},\;X_{k}) + g_D(\mathbf{X}) \]

To ensure the uniqueness and identifiability of each component, MID imposes centering and probability-weighted minimum-norm constraints on the decomposition.

By approximating a black-box model with this surrogate structure, we can derive a representation that retains the superior predictive power of machine learning models without sacrificing actuarial transparency.

Setting Up

We begin by setting up the environment and loading the necessary libraries. We will use a standard actuarial claim dataset (freMTPL2freq in CASdatasets), split into 338,995 training and 338,996 testing observations.

Code 
# data manipulation
library(arrow)
library(dplyr)

# predictive modeling
library(lightgbm)

# surrogate modeling
library(midr)

# visualization
library(ggplot2)
library(patchwork)

# load training and testing datasets
train <- read_parquet("data/train.parquet")
test  <- read_parquet("data/test.parquet")

A key component of our evaluation is the Weighted Mean Poisson Deviance defined as follows:

\[ \text{Loss}(\mathbf{y}, \hat{\mathbf{y}}, \mathbf{w}) = \frac{\sum_{i=1}^n w_i\;d(y_i,\hat{y}_i)}{\sum_{i=1}^n w_i},\quad d(y_i,\hat{y}_i) = 2\left[y_i (\log y_i - \log \hat{y}_i) - (y_i - \hat{y}_i)\right] \]

Code 
# define loss function
mean_poisson_deviance <- function(
    y_true, y_pred, sample_weight = rep(1, length(y_true))
  ) {
  stopifnot(all(y_pred > 0))
  resid <- ifelse(y_true > 0, y_true * log(y_true / y_pred), 0)
  resid <- resid - y_true + y_pred
  2 * sum(resid * sample_weight) / sum(sample_weight)
}

The Black-Box: LightGBM

GBMs such as LightGBM often yield superior predictive power by capturing high-order interactions. However, this accuracy comes at the cost of being a black box.

Code 
# hold out validation dataset
valid_idx <- seq_len(floor(nrow(train) * 0.2))

# create datasets for training
dtrain <- lgb.Dataset(
  data.matrix(select(train[-valid_idx, ], -Frequency, -Exposure)),
  label = train$Frequency[-valid_idx],
  weight = train$Exposure[-valid_idx],
  categorical_feature = c("VehBrand", "VehGas", "Region")
)

dvalid <- lgb.Dataset.create.valid(
  dtrain,
  data.matrix(select(train[ valid_idx, ], -Frequency, -Exposure)),
  label = train$Frequency[ valid_idx],
  weight = train$Exposure[ valid_idx]
)

# model parameters
params_lgb <- list(
  objective = "poisson",
  learning_rate = 0.03188002,
  num_leaves = 30,
  reg_lambda = 0.004201069,
  reg_alpha = 0.2523909,
  colsample_bynode = 0.5552524,
  subsample = 0.5938199,
  min_child_samples = 9,
  min_split_gain = 0.3920509,
  poisson_max_delta_step = 0.8039541
)

set.seed(42)
fit_lgb <- lgb.train(
  params = params_lgb,
  data = dtrain,
  nrounds = 1000L,
  valids = list(eval = dvalid),
  early_stopping_round = 50L,
  verbose = -1L
)
Code 
# evaluate fitted model
pred_fit_lgb <- predict(
  fit_lgb, data.matrix(select(test, -Frequency, -Exposure))
)

deviance <- mean_poisson_deviance(
  y_true = test$Frequency,
  y_pred = pred_fit_lgb,
  sample_weight = test$Exposure
)
cat("Mean Poisson Deviance:", deviance)
Mean Poisson Deviance: 0.4655023

Surrogate Modeling with MID

We use {midr} to replicate the LightGBM model. By including interaction terms in the model formula, we allow the surrogate to capture the joint relationships that the GBM has learned. The goal is to approximate the LightGBM function \(f_{LGB}(\mathbf{x})\) with our interpretable structure \(g^*(\mathbf{x})\):

\[\log f_{LGB}(\mathbf{x}) \approx g^*(\mathbf{x}) = g_\emptyset + \sum_{j} g_j(x_j) + \sum_{j < k} g_{jk}(x_j, x_k)\]

WarningComputational Considerations

Including all second-order interactions using the (...)^2 syntax results in \(p(p-1)/2\) interaction terms. For high-dimensional data, this can be memory-intensive.

Code 
mid_lgb <- interpret(
  Frequency ~ (VehPower + VehAge + DrivAge + LogDensity +
               VehBrand + VehGas + Region)^2,
  data = train,
  lambda = 0.01,
  weights = Exposure,
  link = "log",
  model = fit_lgb,
  pred.fun = function(model, data) {
    newdata <- data.matrix(select(data, -Frequency, -Exposure))
    predict(model, newdata)
  }
)

summary(mid_lgb)

Model Fidelity: The Uninterpreted Ratio (UR)

To measure how successfully our surrogate replicates the LightGBM model, we use the Uninterpreted Ratio (UR) or the R-squared score (\(R^2 = 1 - UR\)) on the link scale (\(\log\) scale).

The MID model achieves a very low uninterpreted ratio. Only 7.16% of the variance in the model predictions on the log-transformed scale is left unexplained.

Code 
# evaluate fitted surrogate on testing set
pred_mid_lgb = predict(mid_lgb, test, type = "response")

# calculate R-squared on testing dataset
R2_test <- weighted.loss(
  x = log(pred_fit_lgb),
  y = log(pred_mid_lgb),
  w = test$Exposure,
  method = "r2"
)
UR_test <- 1 - R2_test

R2_test_response <- weighted.loss(
  x = pred_fit_lgb,
  y = pred_mid_lgb,
  w = test$Exposure,
  method = "r2"
)
UR_test_response <- 1 - R2_test_response

UR_train <- mid_lgb$ratio[1L]
UR_train_response <- mid_lgb$ratio[2L]

cat(
  "Uninterpreted Ratio on the training set:\n",
  sprintf(
    " %.4f in the log-transformed scale\n %.4f in the original scale\n",
    UR_train, UR_train_response
  ),
  "Uninterpreted Ratio on the testing set:\n",
  sprintf(
    " %.4f in the log-transformed scale\n %.4f in the original scale",
    UR_test, UR_test_response
  ),
  sep = ""
)
Uninterpreted Ratio on the training set:
 0.0716 in the log-transformed scale
 0.1815 in the original scale
Uninterpreted Ratio on the testing set:
 0.0705 in the log-transformed scale
 0.2135 in the original scale

Feature Effects

The LightGBM model can now be inspected through a transparent additive approximation. We can visualize each main effect \(g_j\) using ggplot2-based outputs via ggmid(). The model correctly identifies higher risks for young drivers (DrivAge) and highly populated areas (LogDensity).

Code 
wrap_plots(
  lapply(mid.plots(mid_lgb), function(x) x + labs(y = NULL)),
  design = "AAAAABBBBBCCCCCDDDDD\nEEEEEFFFFGGGGGGGGGGG"
)

A key advantage of {midr} is its ability to isolate interaction effects \(g_{jk}\) from main effects \(g_j\). For example, while the main effect of LogDensity is generally increasing, we can learn from the interaction model that the risk remains high even in low-density areas specifically in Île-de-France.

Code 
wrap_plots(
  nrow = 1, widths = c(3, 2),
  ggmid(mid_lgb, "LogDensity:Region", type = "data",
        data = train[1:1e4, ]) +
    labs(y = NULL, subtitle = "Interaction Effect") +
    theme(legend.position = "bottom"),
  ggmid(mid_lgb, "LogDensity:Region", main.effects = TRUE) +
    labs(y = NULL, subtitle = "Total Effect") +
    scale_y_discrete(labels = NULL) +
    theme(legend.position = "bottom")
)

Effect Importance

To rank the influence of each component discovered in the LightGBM model, we calculate the Effect Importance, defined as the average absolute contribution. This is highly useful to know which effects are important and non-negligible.

Code 
imp_lgb <- mid.importance(mid_lgb, data = train, max.nsamples = 2000)
wrap_plots(
  nrow = 1, widths = c(4, 3),
  ggmid(imp_lgb, theme = "bluescale@qual", max.nterms = 20) +
    labs(title = "Effect Importance",
         subtitle = "Average absolute effect per feature") +
    theme(legend.position = "none"),
  ggmid(imp_lgb, type = "boxplot", theme = "mako_r", max.nterms = 20) +
    labs(title = "",
         subtitle = "Distribution of effect per feature") +
    scale_y_discrete(labels = NULL) +
    theme(legend.position = "none")
)

Conditional Expectation

We further explore the model’s behavior using the ICE plot. In the MID framework, the variation in ICE curves for a feature \(j\) is explicitly governed by the interaction terms \(g_{jk}\) identified from the LightGBM model.

Code 
ice_lgb_link <- mid.conditional(mid_lgb, type = "link", variable = "DrivAge")
ice_lgb <- mid.conditional(mid_lgb, variable = "DrivAge")
wrap_plots(
  nrow = 1,
  ggmid(ice_lgb_link, var.color = LogDensity) +
    theme(legend.position = "bottom") +
    labs(y = "Linear Predictor",
         title = "Conditional Expectation",
         subtitle = "Change in linear predictor"),
  ggmid(ice_lgb, type = "centered", var.color = LogDensity) +
    theme(legend.position = "bottom") +
    labs(y = "Prediction", title = "",
         subtitle = "Centered change in original scale")
)

Additive Attribution

Finally, we perform an Additive Breakdown for individual predictions. This provides an exact allocation of the LightGBM’s prediction into the terms of our surrogate model, serving as a powerful tool for instance-level explanation.

Code 
set.seed(42)
row_ids <- sort(sample(nrow(train), 4))

bd_plots <- lapply(
  row_ids,
  function(idx) {
    res <- mid.breakdown(mid_lgb, train, row = idx)
    label <- paste0("Breakdown for Row ", idx)
    ggmid(res, theme = "shap", max.nterms = 10,
          format.args = list(digits = 4)) +
      labs(x = "Linear Predictor", subtitle = label) +
      theme(legend.position = "none")
  }
)

wrap_plots(grobs = bd_plots)

Extensions to Matrix-Valued Predictions

While the previous examples focused on scalar predictions (e.g., expected claim frequency), the MID algorithm can be naturally extended to handle matrix-valued outputs. This extension unlocks several advanced actuarial applications without fundamentally changing the underlying algorithm:

Survival Analysis (SurvMID(t))

In survival analysis, models predict survival curves \(S(t)\) evaluated at \(m\) distinct time points. By applying the matrix-extended MID, we obtain time-dependent main effects \(g_j(x_j, t)\) and interactions \(g_{jk}(x_j, x_k, t)\). This allows actuaries to visualize how a specific risk factor’s influence evolves over the lifetime of a policy, moving beyond static hazard ratios.

Code 
library(survival)
library(aorsf)

fit_orsf <- orsf(
  Surv(time, status) ~ .,
  data = veteran
)

mid_orsf <- interpret(
  Surv(time, status) ~ .,
  data = veteran,
  model = fit_orsf,
  lambda = .1
)

imp_orsf <- mid.importance(mid_orsf)

p1 <- ggmid(imp_orsf, type = "series") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow = 1)) +
  labs(subtitle = "Time-Dependent Feature Importance")
p2 <- ggmid(mid_orsf, term = "karno", type = "series", intercept = TRUE) +
  theme(legend.position = "bottom") +
  labs(subtitle = "Feature Effect of 'karno'", y = expression(S(t)))
p3 <- ggmid(mid_orsf, term = "celltype", type = "series", intercept = TRUE) +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow = 1)) +
  labs(subtitle = "Feature Effect of 'celltype'", y = expression(S(t)))

p1 / (p2 + p3)

Quantile Regression

Instead of predicting the conditional mean, quantile regression predicts values at \(m\) different quantile levels (e.g., 50th, 75th, 95th, and 99th percentiles). Applying MID to this \(n \times m\) output matrix helps identify which variables drive the median behavior versus which variables govern extreme tail risks (e.g., the severity of unusually large claims).

Code 
library(ranger)

fit_rf <- ranger(
  mpg ~ .,
  data = mtcars,
  quantreg = TRUE
)

mid_rf <- interpret(
  mpg ~ .,
  data = mtcars,
  model = fit_rf,
  lambda = 1,
  pred.args = list(
    type = "quantiles",
    quantiles = c(1:19 / 20),
    target = NULL # required for the current version of midr::get.yhat()
  )
)

labels(mid_rf) <- c(1:19 / 20) * 100
p1 <- ggmid(mid_rf, "hp", theme = "shap", intercept = TRUE) +
  labs(color = "quantile", subtitle = "Feature Effect of 'hp' on Quantile Predictions") +
  theme(legend.position = "bottom")
p2 <- ggmid(mid_rf, "wt", theme = "shap", intercept = TRUE) +
  labs(color = "quantile", subtitle = "Feature Effect of 'wt' on Quantile Predictions") +
  theme(legend.position = "bottom")
p1 + p2

Conclusion

In this notebook, we have demonstrated how Maximum Interpretation Decomposition (MID) bridges the gap between predictive performance and model transparency. By using the {midr} package, we successfully transformed a complex LightGBM model into a structured, additive representation.

Key takeaways include:

  • A pragmatic functional decomposition: MID is defined as the minimizer of the Uninterpreted Ratio (UR).
  • High-fidelity, inherently interpretable global surrogate: Fidelity is explicitly quantified by UR, providing a defensible basis for discussions with management, auditors, and regulators.
  • Highly extensible: It extends naturally to matrix-valued outputs, e.g., survival curves (SurvMID(t)).
  • Ready to use today: Available via {midr}, {midnight}, and {midlearn}.

As machine learning models become increasingly prevalent in insurance pricing and reserving, tools like MID will be essential for ensuring that our “black-boxes” remain accountable, reliable, and fundamentally understood.