# Explainable Boosting Machine for Predicting Claim Severity and Frequency in Car Insurance

Marketa Krupova, Nabil Rachdi, Quentin Guibert

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marketa Krupova, Nabil Rachdi, Quentin Guibert. Explainable Boosting Machine for Predicting Claim Severity and Frequency in Car Insurance. 2025. hal-05007738

HAL Id: hal-05007738

https://hal.science/hal-05007738v1

Preprint submitted on 26 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Explainable Boosting Machine for Predicting Claim Severity and Frequency in Car Insurance

Markéta Krúpová*¹, Nabil Rachdi†¹, and Quentin Guibert‡²

¹ADDACTIS France, 46 bis Chemin du Vieux Moulin, 69160 Tassin-la-Demi-Lune, France

²CEREMADE, Université Paris-Dauphine, PSL University, CNRS, 75016 Paris, France

March 27, 2025

## Abstract

In a context of constant increase in competition and heightened regulatory pressure, accuracy, actuarial precision, as well as transparency and understanding of the tariff, are key issues in non-life insurance. Traditionally used generalized linear models (GLM) result in a multiplicative tariff that favors interpretability. With the rapid development of machine learning and deep learning techniques, actuaries and the rest of the insurance industry have adopted these techniques widely. However, there is a need to associate them with interpretability techniques. In this paper, our study focuses on introducing an Explainable Boosting Machine (EBM) model that combines intrinsically interpretable characteristics and high prediction performance. This approach is described as a glass-box model and relies on the use of a Generalized Additive Model (GAM) and a cyclic gradient boosting algorithm. It accounts for univariate and pairwise interaction effects between features and provides naturally explanations on them. We implement this approach on car insurance frequency and severity data and extensively compare the performance of this approach with classical competitors: a GLM, a GAM, a CART model and an Extreme Gradient Boosting (XGB) algorithm. Finally, we examine the interpretability of these models to capture the main determinants of claim costs.

**Keywords:** car insurance, explainable boosting machine, generalized additive model, cyclic gradient boosting, interpretable machine learning, glass-box.

## 1 Introduction

Recent years have been marked by a massive development of artificial intelligence (AI) techniques in data processing and statistical modeling. The transformations introduced by these technologies cover a wide range of industrial applications, including the processing and analysis of big data, digital marketing, and the generation of texts, sounds, images, and videos. The insurance industry and actuaries are naturally affected as they increasingly integrate these tools into their daily work (Eling et al., 2022; Richman, 2024). One of the main challenges in insurance remains the implementation of competitive pricing models that can accurately and finely reflect the risks underwritten. In non-life insurance, this involves determining, in accordance with the principles of actuarial science, the average prediction of the claim burden of policyholders. Taking into account the characteristics of the insured risk, actuaries usually examine the frequency and severity of claims independently, based on the historical data collected by the insurer.

Historically, actuarial practices have moved from the use of traditional generalized linear models (GLM) (Nelder and Wedderburn, 1972; McCullagh and Nelder, 1989) to risk pricing models enriched with the latest statistical machine learning (ML) and deep learning (DL) techniques, see

<footnote>*Email: marketa.krupova@addactis.com.</footnote>

<footnote>†Email: nabil.rachdi@addactis.com.</footnote>

<footnote>‡Email: guibert@ceremade.dauphine.fr.</footnote>
Richman (2021b), Richman (2021a), Blier-Wong et al. (2021), and Embrechts and Wüthrich (2022) for the recent adoption of AI by actuaries and Holvoet et al. (2024) for an overview of the recent literature for insurance pricing models with ML and DL techniques. Indeed, these techniques generally outperform traditional pricing methods on tabular data by capturing complex non-linear relationships between variables. As a result, they improve the predictive performance of insurance models, making them a key determinant of competition in the insurance industry. Among others, Henckaerts et al. (2021), Poufinas et al. (2023) and Power et al. (2024) examine the performance of bagging and boosting methods such as random forest or Extreme Gradient Boosting (XGB) for pricing car insurance contracts. The supervised ML techniques used in these studies improve individual prediction performance with high-dimensional data, and notably with telematics data (Pesantez-Narvaez et al., 2019; Maillart, 2021). DL models, i.e. neural network models with an architecture composed of multiple layers, are also used to learn complex data structures in insurance (Denuit et al., 2019). Noll et al. (2020) compare neural network techniques with other ML techniques on motor third-party liability insurance data. Delong and Kozak (2023) develop autoencoders for categorical data. Holvoet et al. (2024) use neural networks and combined actuarial neural networks (CANN) in frequency-severity modeling to improve tariff performance.

Despite their improved performance, some ML and DL models can have the disadvantage of being black-boxes, making it difficult to interpret and provide transparency in their results, see e.g. Pasquale (2015). This is a major challenge in insurance, given the regulatory obligations related to personal data protection (GDPR) (OCDE, 2020) and individualized risk analysis. Like other highly regulated industries, insurance must ensure the explainability, fairness, and transparency of algorithms (Owens et al., 2022; Henckaerts et al., 2021; Henckaerts et al., 2022). In particular, the ability to understand and explain individual insurance rates is essential to ensure consistency with covered risk classes and to maintain actuarial fairness (Meyers and Van Hoyweghen, 2018). In addition, insurers must be vigilant about potential biases and discrimination that their pricing methods might introduce (Charpentier, 2024). In particular, black-box predictive models are inherently difficult to interpret, which requires careful consideration of the effects of certain factors or interactions between variables to avoid implicit discrimination (Lindholm et al., 2022; Frees and Huang, 2023). Thus, insurers face a trade-off in the selection of predictive methods between performance and interpretability. As a result, fully interpretable models such as GLMs remain widely implemented by insurers despite the superior performance of black-box models with high-dimensional feature effects.

Significant efforts have been made in recent years to enhance the transparency, explainability, and interpretability of models with the development of explainable artificial intelligence (XAI) algorithms, see Molnar et al. (2020) and Owens et al. (2022) for a review. When applied in insurance, these tools help to better understand the relationships, both globally and individually, between the target variable and the input features. They also enable the ranking of predictor importance and the evaluation of interactions between them. In general, two main types of algorithms are distinguished (Delcaillau et al., 2022). On the one hand, intrinsically interpretable models such as linear models or decision trees possess ex-ante interpretability, also referred to as model-based interpretability. On the other hand, post-hoc interpretability methods require ex-post analysis and interpretation tools such as the Accumulated Local Effects (ALE) method, the Locally Interpretable Model-agnostic Explanation (LIME) or SHapley Additive exPlanations (SHAP), see Lorentzen and Mayer (2020), Delcaillau et al. (2022), or Charpentier (2024, Chapter 4) for an introduction to these tools. Another approach for globally analyzing the predictions of a black-box model consists of using an interpretable surrogate model to shed light on the predictions by approximating its behavior, see e.g. Henckaerts et al. (2022) and Holvoet et al. (2024) for applications with surrogate GLMs or Maillart (2021) for surrogate regression trees. However, implementing this technique requires a complex procedure to transition from a black-box model to an interpretable surrogate model. Furthermore, its effectiveness depends on how well the surrogate model can approximate the original predictive model.

In this paper, we examine the implementation of the Explainable Boosting Machine (EBM) model for non-life insurance pricing. Unlike the post-hoc approaches described above, the EBM is initially developed by Lou et al. (2012) and later refined by Caruana et al. (2015) and Nori et al. (2019). It is based on generalized additive models (GAM) (Hastie and Tibshirani, 1990) and bagged-boosted decision trees for fitting shape functions, ensuring interpretability through a linear relationship between predictions and original features. More precisely, it employs a cyclic gradient boosting method, a technique based on the classical cyclic coordinate descent optimization method (Luenberger and Ye, 2021) recently revisited by Delong et al. (2023). This approach offers the advantage of being directly interpretable, making it more suitable for decision-making as recommended by Rudin (2019). EBM is classified as a glass-box model, i.e. a ML or a DL model that
emphasizes interpretability and provides transparency into its decision-making processes. The use of shape functions enables visualization of the local contribution of each variable and facilitates the computation of feature importance scores, which is valuable for model interpretation and variable selection. Moreover, the EBM method incorporates the GA²Ms selection algorithm (Lou et al., 2013) for the identification and selection of pairwise interaction terms.

Our approach aligns with a recent stream of literature focused on the development of AI models that are accurate and directly interpretable (Saleem et al., 2022). Given their flexibility, GLM and GAM models often serve as a foundational framework for developing more complex models on tabular data (Wood et al., 2022). The GAMboostLSS model introduced by Mayr et al. (2012) and implemented by Hofner et al. (2016) is a GAM for Location, Scale, and Shape that employs a cyclic boosting method for fitting. However, this approach is less interpretable than the EBM model, which allows for the construction of independent shape functions for each variable. Richman and Wüthrich (2023) recently proposed the LocalGLMnet model, a deep learning approach based on a fully connected feed-forward neural network (FFN), where the response variable is predicted as a local GLM with coefficients learned from the FFN. Building on this idea, Zakrisson and Lindholm (2025) introduced a varying coefficient model where coefficients are trained locally using a cyclic gradient boosting method, a framework closely related to that of the EBM model. Furthermore, Maillart and Robert (2024) developed a GAM model incorporating shape functions derived from an ensemble of decision trees and a GLM applied to pre-binned variables. Agarwal et al. (2021) and Thielmann et al. (2024) proposed, respectively, a Neural Additive Model and a Neural Additive Model for Location, Scale, and Shape, both inspired by the structure of a GAM. These neural network models are closely related to the EBM and GAMboostLSS models and offer an additive structure, which simplifies feature interpretability. These approaches rely on a specialized network architecture, which achieves performance comparable to that of the EBM model. Another neural network-based approach is the TabNet algorithm (Arik and Pfister, 2021), which employs a sequential attention mechanism to select features across multiple decision steps. This algorithm ensures both global and local interpretability of its predictions but suffers from high computational costs. It was recently applied to telematics insurance data by McDonnell et al. (2023).

The EBM model introduced in this paper investigates the contribution of glass-box models in analyzing claim frequency and severity data in car insurance. Our contribution is structured around three main aspects. First, we introduce and adapt the EBM model to frequency and severity data, thereby enriching the set of interpretable models available in the insurance framework. Despite its relevance and ease of use, this model has, to the best of our knowledge, been rarely applied in the insurance literature. Second, we assess the predictive performance of our model and compare it against four benchmark regression models commonly used for insurance pricing: a Generalized Linear Model (GLM), a classical Generalized Additive Model (GAM), a regression tree model (CART) and an Extreme Gradient Boosting (XGB) algorithm. Third, we examine the interpretability benefits of the EBM model for frequency and severity predictions. Since this model is inherently interpretable, it does not require post-hoc methods or surrogate models to explain its results (Henckaerts et al., 2022), whether at a local or global level. A key advantage of the EBM model lies in its ability to capture and express both individual and pairwise feature contributions. In addition, we compare these results with a post-hoc analysis of the XGB model's outcomes.

This paper is organized as follows. Section 2 recalls the fundamental principles of GLM and GAM models. Section 3 describes the EBM method. Section 4 presents the application of this model to frequency and severity data in car insurance, where we assess its performance compared to other regression models. Section 5 examines its usefulness for interpretability purposes. Finally, Section 6 concludes this study.

## 2 Non-life insurance pricing with GLM and GAM

The starting point of our modeling relies on the structure proposed by generalized additive models (GAMs), which themselves build upon generalized linear models (GLMs). We begin with an introduction to the notation used in Section 2.1 and a brief description of GLMs and GAMs in Section 2.2. The framework introduced by GLMs and GAMs is commonly employed by actuaries for the pricing of non-life insurance products due to their ease of use, interpretability of the model, and explainability of the obtained results.
## 2.1 Notations

In all that follows, we consider a target space $\mathcal{Y}$, e.g. $\mathcal{Y} = \{0,1\}$ for classification or $\mathcal{Y} = \mathbb{R}$ for regression, and a $p$-dimensional feature space $\mathcal{X}$, where $p \in \mathbb{N}^*$. Let $Y \in \mathcal{Y}$ denote the target variable and $X_1, ..., X_p$ the $p$ features in $\mathcal{X}$. We denote $\mathbf{x} = (x_1, ..., x_p)$ an observation of the variable $\mathbf{X} = (X_1, ..., X_p)$ and $y$ an observation of the target variable $Y$. Our context is one of supervised learning with a statistical sample of $n \in \mathbb{N}^*$ observations $\mathcal{D} = \{(\mathbf{x}_i, Y_i); i \in I\}$ of an unknown joint law $\mathcal{P}$ on $\mathcal{X} \times \mathcal{Y}$, where $\mathcal{X} \subset \mathbb{R}^p$, $\mathcal{Y} \subset \mathbb{R}_+$ and $I = \{1, ..., n\}$. The chosen framework is therefore that of regression: for the claim severity, positive real values is predicted, while for the claim frequency, natural integers is targeted. For clarity, bold notation is used for vectors and matrices.

A prediction rule is a measurable function $F: \mathcal{X} \to \mathcal{Y}$ which associates the output $F(\mathbf{x})$ with the input $\mathbf{x} \in \mathcal{X}$. To measure the quality of the prediction, a non-negative convex loss function $L: \mathcal{Y} \times \mathcal{Y} \to \mathbb{R}_+$ such that $L(y, y) = 0$ and $L(y, y') > 0$ for $y \neq y'$ is introduced. Given a loss function $L$ and a prediction rule $F$, the generalization risk or error is defined by $R_P(F) = \mathbb{E}_{(\mathbf{x},Y) \sim \mathcal{P}}[L(Y, F(\mathbf{x}))]$ and describes the average behavior of the loss function.

## 2.2 GLM and GAM models

In this section, we present the principles and notations necessary for applying generalized linear models (GLM) (Section 2.2.1) and generalized additive models (GAM) (Section 2.2.2).

### 2.2.1 Generalized linear models

Generalized linear models (Nelder and Wedderburn, 1972; McCullagh and Nelder, 1989) are characterized by three elements: a random component, a deterministic component and a functional relationship between the random and deterministic components, see e.g. Ohlsson and Johansson (2010) and Wüthrich and Merz (2023) for comprehensive presentation with insurance data.

The random component characterizes the probability distribution of the target $Y$. It is assumed that $Y_1, ..., Y_n$ are independent random variables belonging to the exponential dispersion family (EDF), which includes most common distributions, notably Gaussian, Binomial, Poisson, and gamma distributions. Note that the gamma distribution is particularly suitable to model unimodal variables with positive values, such as life expectancy or the body of the distribution of claim severity. The Poisson distribution is commonly used to model count data, such as claim frequency.

The log-likelihood associated with a single observation $y_i$ of the random variable $Y_i$, $i \in I$, is given by

$$ \ln \mathcal{L}(\lambda_i, \phi, y_i) = \frac{\lambda_i y_i - b(\lambda_i)}{a(\phi)} + c(y_i, \phi), \quad y_i \in \mathcal{Y} \subset \mathbb{R}, \quad \lambda_i \in \Lambda \subset \mathbb{R}, \quad (1) $$

and $-\infty$ if $y_i \notin \mathcal{Y}$, where $a: \mathbb{R} \to \mathbb{R}$, $b: \Lambda \to \mathbb{R}$, and $c: \mathbb{Y} \times \mathbb{R} \to \mathbb{R}$ are known real-valued differentiable measurable functions. The term $\lambda_i$ is called the natural parameter and $\phi$ is the dispersion parameter. The function $b$ is the cumulant function, which is twice differentiable.

A GLM introduces a regression relationship between the expectation of $Y_i$, i.e. $\mathbb{E}[Y_i] = b'(\lambda_i)$, and $\mathbf{x}_i$, which assumes for all $\beta \in \mathbb{R}^p$ and for all $i \in I$

$$ g(b'(\lambda_i(\beta_0, \beta))) = \eta_i = \beta_0 + \langle \mathbf{x}_i, \beta \rangle \\ = \underbrace{\beta_0}_{\text{global average}} + \underbrace{\beta_1 x_{i,1}}_{\text{X}_1 \text{ effect}} + \underbrace{\beta_2 x_{i,2}}_{\text{X}_2 \text{ effect}} + \dots + \underbrace{\beta_p x_{i,p}}_{\text{X}_p \text{ effect}}, \quad (2) $$

where $\eta_i$ is the linear predictor, $\langle \cdot, \cdot \rangle$ denotes the scalar product, and $\lambda_i(\beta_0, \beta)$ stresses the dependence on the finite-dimensional parameter $\beta$ and an intercept $\beta_0$. Assuming $g$ is a bijective function, we have $\lambda_i(\beta) = \tilde{b}(g^{-1}(\eta_i))$ where $\tilde{b} = (b')^{-1}$ is the inverse of $b'$. The link function $g$ is called the canonical link when $g = \tilde{b}$. The parameters $\beta_0, ..., \beta_p$ are estimated using the maximum likelihood method.

The linear relationship expressed in Equation (2) allows for an easy interpretation of the effects of different features in an additive form. In particular, the use of the logarithmic link function in the case of the gamma and Poisson distributions leads to a multiplicative structure traditionally used in actuarial modeling, which corresponds to a solution making the predicted tariff easily interpretable.

### 2.2.2 Inherently interpretable generalized additive models

GAMs (Hastie and Tibshirani, 1990) are an extension of GLMs, where the expectation of the response variable $Y$ is explained as a function of a linear combination of predictors. Usually, individual
contribution of each predictor is considered and the deterministic components (2) related to some continuous features of generalized linear models can be replaced by non-parametric components

$$g(\mathbb{E}[Y_i | \mathbf{X}_i = \mathbf{x}_i]) = \beta_0 + \sum_{j=1}^{p_1} \beta_j x_{i,j} + \sum_{j=p_1+1}^{p_2} f_j(x_{i,j}), \quad \forall i \in I, \qquad (3)$$

with $x_{i,1}, \dots, x_{i,p_1}$ are $p_1 \in \mathbb{N}^*$ continuous or categorical features, and $x_{i,p_1+1}, \dots, x_{i,p_2}$ are continuous features such as $p = p_1 + p_2$. The target $Y$ still belongs to the EDF, but the dependency relationship with the features is complexified by the addition of measurable shape functions $f_j$, $j \in \{1, \dots, p_2\}$, such that $\mathbb{E}[f_j] = 0$ and $\mathbb{E}[f_j^2] < \infty$. More precisely, this GAM formulation introduces one-dimensional components $f_1, \dots, f_p$ which remain inherently interpretable due to their additive form

$$\sum_{j=1}^{p_2} f_j(x_{i,j}) = \underbrace{f_1(x_{i,1})}_{\text{X}_1 \text{ effect}} + \underbrace{f_2(x_{i,2})}_{\text{X}_2 \text{ effect}} + \dots + \underbrace{f_{p_2}(x_{i,p_2})}_{\text{X}_{p_2} \text{ effect}}.$$

Maintaining the inherently interpretable additive structure allows for a parsimonious addition of non-linearity into the model: each feature can have a complex non-linear shape but relates a single attribute to the target (Wood et al., 2022). This unidimensional specification for shape functions is traditionally used in actuarial literature where smooth functions are usually considered for incorporating non-linear effects related to continuous features (Denuit et al., 2019). A variety of smoothing spline functions (Wood, 2006) can be retained in non-life pricing problem, see e.g. Denuit and Lang (2004), or for mortality graduation, see e.g. Denuit and Legrand (2018).

Two main approaches are possible for fitting GAMs with spline functions. First, we can rewrite a spline shape functions as

$$f_j(x_{i,j}) = \sum_{q=1}^{Q_j} \beta_{jq} b_{jq}(x_{i,j}), \quad \forall i \in I, \quad \forall j = 1, \dots, p_2,$$

with $b_{jq}$ as a spline function and $Q_j \in \mathbb{N}^*$. Subsequently, the estimation of the coefficients can be based on maximizing the penalized log-likelihood, e.g. using penalized iteratively reweighted least squares. Alternatively, least square algorithms can be used. The implementation of these methods is available in the mgcv R-package (Wood, 2023).

Similarly to GLM, an interaction component can be added to the formula (3) to capture multivariate relationships. For instance, a bivariate smooth function $f_{(j_1,j_2)}(x_{i,j_1}, x_{i,j_2})$ is used to capture spatial effects, see. e.g. Henckaerts et al. (2018). Such bivariate smooth functions also correspond to the initial idea implemented in the GA$^2$Ms algorithm proposed by Lou et al. (2013) for incorporating pairwise interaction terms.

# 3 The Explainable Boosting Machine methodology

Although it is rarely used in actuarial literature to our knowledge, the usual smooth functions of a GAM can be replaced by more general shape functions, such as regression trees or ensembles of decision trees. The Explainable Boosting Machine (EBM) presented in this section is a GAM with univariate and bivariate components as in Equation (3) where each $f_j$ or $f_{(j_1,j_2)}$ shape function is a machine learning technique. Section 3.1 briefly presents the introduction of shape functions in the GAM framework to obtain the EBM. Then, we explain in Section 3.2 how to train this model based on a cyclic boosting algorithm.

## 3.1 A GAM with machine learning shape functions

In order to devise GAM models closer in accuracy to the full-complexity models such as random forests, boosted trees or neural nets while remaining intrinsically interpretable, univariate shape functions using machine learning techniques as building blocks are first introduced. The accuracy of machine learning models is combined with the inherent intelligibility of GAM, yielding a family of glass-box models. Several specifications for shape functions can be considered based on tree structure (Lou et al., 2012):

* Single decision trees, as illustrated in Figure 1 with CART tree (Breiman et al., 1984),

* Independent learning of several decision trees, then averaging the results using a bagging procedure,
* Adaptive learning of several decision trees, then aggregating the results using a boosting procedure,

* Combination of the bagging and boosting procedures, named bagged-boosted hereafter, with decision trees.

These shape functions are more accurate than spline functions and are not limited to continuous features. In particular, Lou et al. (2012) empirically prove the good accuracy of the bagged-boosted procedure. With these tree ensemble shape functions, the GAM is then fitted using backfitting (Hastie and Tibshirani, 1990) or gradient boosting (Friedman, 2001).

$$g(\mathbb{E}[Y | \mathbf{X} = \mathbf{x}]) = \beta_0 + \underset{\text{ CART}(x_1, Y)}{\boxed{\phantom{\beta_0}}} + \underset{\text{ CART}(x_2, Y)}{\boxed{\phantom{\beta_0}}} + \dots + \underset{\text{ CART}(x_p, Y)}{\boxed{\phantom{\beta_0}}}$$

<caption>Figure 1: GAM combining CART tree shape functions</caption>

In addition, Lou et al. (2013) efficiently extend this approach with two-dimensional interaction shape functions called the GA²Ms method, which remains interpretable by plotting heatmaps. As a shape function can be applied to all features, the model (3) can then be expressed in the general form

$$g(\mathbb{E}[Y_i | \mathbf{X}_i = \mathbf{x}_i]) = \beta_0 + \sum_{j=1}^{p} f_j(x_{i,j}) + \sum_{j_1=1}^{p} \sum_{j_2 \neq j_1}^{p} f_{(j_1, j_2)}(x_{i,j_1}, x_{i,j_2}), \quad \forall i \in I, \quad (4)$$

The authors propose an interaction detection algorithm, called FAST, for the computationally efficient selection of relevant pairwise interactions.

To simplify the notation without loss of generality, we omit the bivariate terms in the model Equation (4). Indeed, the binning procedure described in Section 3.2 operates in the same manner regardless of whether the shape functions are univariate or bivariate. Using the notations from Section 2.1, the prediction rule can then be written as

$$F(\mathbf{x}_i) = g^{-1}(\mathfrak{F}(\mathbf{x}_i)) = g^{-1} \left( \beta_0 + \sum_{j=1}^{p} f_j(x_{i,j}) \right), \quad (5)$$

where $\mathfrak{F}(\mathbf{x}_i)$ corresponds to the score of observation $i \in I$.

## 3.2 Fitting shape functions with a cyclic boosting algorithm

The EBM algorithm (Caruana et al., 2015; Nori et al., 2019) is a fast implementation of the model (4), using bagged-boosted shallow CART trees as shape functions. It is a glass box model, i.e. an intrinsically interpretable machine learning model. By its construction, the EBM algorithm ensures accurate local predictions and, in particular, accounts for rare observations. The training procedure for this model is based on a cyclic boosting algorithm, similar to that proposed by Wick et al. (2019). The cyclic boosting algorithm is an iterative procedure that relies on a cyclic coordinate descent algorithm (Wright, 2015; Luenberger and Ye, 2021), in which each shape function associated with a feature $X_j, j = 1, \dots, p$, is fitted at one time in a round-robin fashion using a low learning rate. Contrary to usual boosting algorithms (Breiman, 2001) where all features are considered simultaneously, this procedure ensures that the order in which the features are estimated does not matter and limits the effects of collinearity.

**Binning procedure** The cyclic boosting algorithm involves segmenting each feature into K-bins. Continuous features are discretized using either a quantile-based or equi-distribution approach to ensure bins of equal size. For categorical variables, the bins correspond to their distinct categories. For bivariate variables resulting from the interaction of two features, the bins represent the interactions of the bins of the two combined variables. Let $b_{j,k}, k \in \{1, \dots, K\}$ denote the k-th bin associated with the j-th feature $X_j$. For each bin, we calculate its weight $w_{j,k}$ as the number of observations included in this bin. The number K of bins is identical for each univariate feature. Note that the bins are not updated during the training procedure.
**Algorithm 1 Cyclic boosting algorithm for EBM with regression tree as shape function**

1: **Input:** training data $(\mathbf{x}_i, y_i)_{i=1,...,n}$, $\nu$, $T$, $d$, $L$

2: **Initialization**

3: Set $f_{j,k}^{(0)}(\mathbf{x}_i) = 0$, for $k \in \{1,..., K\}$, $j \in \{1,..., p\}$

4: Set $F^{(0)}(\mathbf{x}_i) = c^*$, where $c^* = \arg \min_c \sum_{i=1}^n L(y_i, c)$

5: Compute $\beta_0 = g(c)$

6: **for** $t = 1,..., T$ **do**

7: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1,..., p$ **do**

8: &nbsp;&nbsp;&nbsp;&nbsp;**for** $k = 1,..., K$ **do**

9: Compute pseudo-residuals for all $i$ such that $x_{i,j} \in b_{j,k}$

$$ \epsilon_i^{(t)} = - \frac{\partial L(y_i, F(\mathbf{x}_i))}{\partial F(\mathbf{x}_i)} \bigg|_{F(\mathbf{x}_i) = F^{(t-1)}(\mathbf{x}_i)} . $$

10: Fit a regression tree of depth $d$ to pseudo-residuals $(x_{i,j}, \epsilon_i^{(t)})$

11: to obtain regions $R_{j,k,r}^{(t)}$ for $r \in \{1,..., M_{j,k}^{(t)}\}$

12: **for** $r = 1,..., M_{j,k}^{(t)}$ **do**

13: Compute prediction

$$ \gamma_{j,k,r}^{(t)} = \arg \min_{\gamma} \sum_{i:x_{i,j} \in R_{j,k,r}^{(t)} \cap b_{k,j}} L(y_i, g^{-1}(\mathfrak{F}^{(t-1)}(\mathbf{x}_i) + \gamma)). $$

14: **end for**

15: Update bin factor

$$ f_{j,k}^{(t)}(\mathbf{x}_i) = f_{j,k}^{(t-1)}(\mathbf{x}_i) + \nu \sum_{r=1}^{M_{j,k}^{(t)}} \gamma_{j,k,r}^{(t)} \mathbf{1}_{\{x_{i,j} \in R_{j,k,r}^{(t)} \cap b_{k,j}\}}. $$

16: **end for**

17: Update shape function $f_j^{(t)}(\mathbf{x}_i) = \sum_{k=1}^K f_{j,k}^{(t)}(\mathbf{x}_i) \mathbf{1}_{\{x_{i,j} \in b_{j,k}\}}$

18: **end for**

19: Update prediction $F^{(t)}(\mathbf{x}_i) = g^{-1}(\beta_0 + \sum_{j=1}^p f_j^{(t)}(\mathbf{x}_i))$

20: **end for**

21: **Return** $F(\mathbf{x}_i) = F^{(T)}(\mathbf{x}_i)$
**Cyclic boosting algorithm for EBM** Hereto, we outline the main steps of this algorithm. Let $t \in \{1, \dots, T\}$ denote the current iteration and initially set each component of the shape functions to zero. Algorithm 1 is initialized by determining the constant that minimizes the loss function. At each boosting iteration $t$, each feature $X_j$, $j \in \{1, \dots, p\}$, is partitioned into bins $b_{j,k}$, $k \in \{1, \dots, K\}$, allowing shape functions to be trained independently in a cyclical manner. Pseudo-residuals are then computed using the negative gradient of the loss function for all observations $i$ such that $x_{i,j}$ belongs to bin $b_{j,k}$

$$e_i^{(t)} = - \frac{\partial L(y_i, F(\mathbf{x}_i))}{\partial F(\mathbf{x}_i)} \bigg|_{F(\mathbf{x}_i) = F^{(t-1)}(\mathbf{x}_i)}$$

We then fit a bagged-boosted shallow decision tree of depth $d$ to these pseudo-residuals. This results in a partitioning of the training subset defined by $b_{j,k}$ into $M_{j,k}^{(t)}$ regions based on different features. Each factor $f_{j,k}^{(t)}(\mathbf{x}_i)$ represents the contribution to the score of feature $X_j$ for observations $i$ in bin $b_{j,k}$. These factors are updated in each region $R_{j,k,r}^{(t)}$ by adding a constant $\gamma_{j,k,r}^{(t)}$, whose value is determined by locally minimizing the loss function. The learning rate of the method is controlled by the hyperparameter $\nu$, which helps mitigate the risk of overfitting. The factors $f_{j,k}^{(t)}(\mathbf{x}_i)$ are computed for feature $X_j$ at iteration $t$ as

$$f_{j,k}^{(t)}(\mathbf{x}_i) = f_{j,k}^{(t-1)}(\mathbf{x}_i) + \nu \sum_{r=1}^{M_{j,k}^{(t)}} \gamma_{j,k,r}^{(t)} \mathbf{1}_{\{x_{i,j} \in R_{j,k,r}^{(t)} \cap b_{k,j}\}}$$

The update of the shape function corresponding to feature $X_j$ is then expressed as the sum of the contributions from each bin $b_{j,k}$

$$f_j^{(t)}(\mathbf{x}_i) = \sum_{k=1}^{K} f_{j,k}^{(t)}(\mathbf{x}_i) \mathbf{1}_{\{x_{i,j} \in b_{j,k}\}}$$

Finally, the prediction function $F$ is updated using Equation (5) once a full cycle over all variables has been independently completed.

<figure><img src="image_3.png" /><figcaption>Figure 2: Learning procedure in the Explainable Boosting Machine algorithm</figcaption></figure>

This procedure follows the same logic of cycling through the features one at a time and is illustrated in Figure 2. The binning procedure, if performed correctly, enables the prediction of rare events in the data. As optimization is performed locally in each bin $b_{j,k}$, atypical observations are treated separately for each factor $X_j$. The handling of rare effects marks a difference between the EBM algorithm and other machine learning methods, which tend to over-regularize them.

The Algorithm 1 converges when using convex loss functions, as is the case for log-likelihood functions belonging to the exponential dispersion family. Similar to GAMboostLSS (Hofner et al., 2016), the number of iterations in the boosting procedure is controlled by an early stopping
hyperparameter, which halts training after a predefined number of consecutive iterations without improvement determined by a tolerance threshold. This hyperparameter is defined globally, in contrast to Delong et al. (2023) who introduced a dimension-wise early stopping criterion.

# 4 Application with claim frequency and severity data

The objective of this section is to fairly evaluate the contributions of the EBM methodology by applying it with claim frequency and severity data in car insurance. We describe the selected dataset in Section 4.1 and the training procedure of EBM in Section 4.2. Section 4.3 presents our EBM specification. To enable an assessment of EBM, we introduce four benchmark machine learning models in Section 4.4. We provide a nuanced view of the competing benchmark model performance by establishing a model evaluation framework (Section 4.5). The numerical evaluations are performed using the R (R Core Team, 2024) and Python software.

## 4.1 Data description

The modeling work of this article relies on a car insurance dataset used for research and development purposes within ADDACTIS France. The warranty modeled is the full accidental damage cover over the period from 2012 to 2021. Risk is assumed to be homogeneous over this period. However, a time variable is included in the study to capture the non-structural trend in both claim frequency and severity. The claim frequency dataset contains roughly 5 million observations and the claim severity dataset for 100 thousand observations. Data analysis work precedes the exploitation of these datasets and can be consulted in detail in Krúpová (2023). Particular attention is paid to the feature selection process, with correlation analysis, factor analysis, and supervised feature selection. Based on this previous study, we select a limited number of features, relevant from both a statistical and a business standpoint for the claim frequency and severity modeling.

Table 1 describes all the variables of this dataset, as well as their use for claim frequency and severity applications. Summary statistics for each feature can be found in Appendix B for claim frequency and claim severity. The selected features can be schematically divided into four categories:

1. Characteristics of the policyholder (Driver Age, License Seniority, Secondary Novice Driver, CRM Coefficient<fnref n="1" />),

2. Characteristics of the vehicle (Vehicle Price Class, Vehicle Age, Vehicle Use),

3. Geographical zone (Zoning Plan),

4. Temporality (Year).

For each policyholder $i \in I$, the target variable is the claim amount for the claim severity dataset and the claim count for the claim frequency dataset. Additionally, the claim frequency dataset contains a variable for exposure-to-risk.

## 4.2 Training and fine tuning procedures

The training procedure relies on the distinction of training and testing datasets: 80% of the dataset is used for training the model and the remaining 20% is used for evaluating its generalizability. To guarantee objectivity and comparability of results, the training and testing datasets are built to have the same characteristics, i.e. similar distributions for each variable. We denote respectively by $\mathcal{D}_g^{\text{train}}$ and $\mathcal{D}_g^{\text{test}}$ the training and testing sets for $g \in \{\text{freq, sev}\}$ for the frequency and severity datasets.

ML models require simultaneous optimization of several parameters in the training procedure. These parameters encapsulate the complexity of the model and therefore need particular attention. Grid search is an optimization method that exhaustively tries every combination of the provided parameters to find the best model. However, grid search optimization is computationally expensive and inefficient on larger datasets. Random search provides an alternative by selecting combinations of parameters to try (Bergstra and Bengio, 2012).

<footnote n="1">¹French system that rewards or penalizes drivers based on their driving history. The initial coefficient is 1 and for each year without accident, the coefficient is reduced by 5%. After 13 years without accident, the coefficient reaches 0.5 and reduces the initial premium by 50%. Conversely, in the event of an accident, the coefficient is increased depending on the severity of the accident and the driver's responsibility and penalizes the initial premium.</footnote>
<table><caption>Table 1: Features selected for the claim frequency and severity modeling of the full accidental damage cover.</caption><thead><tr><th>Feature</th><th>Description</th><th>Frequency</th><th>Severity</th></tr></thead><tbody><tr><td>Driver Age</td><td>Age of the main driver in years (from ages 17 to 100)</td><td>✓</td><td>✓</td></tr><tr><td>License Seniority</td><td>License seniority of the main driver in years (from 0 to 82 years)</td><td>✓</td><td>✓</td></tr><tr><td>Secondary Novice Driver</td><td>Presence of a secondary novice driver</td><td>✓</td><td>✓</td></tr><tr><td>CRM Coefficient</td><td>Reduction-increase (bonus-malus) coefficient</td><td>✓</td><td>✓</td></tr><tr><td>Vehicle Price Class</td><td>Cover-specific vehicle price class</td><td>✓</td><td>✓</td></tr><tr><td>Vehicle Age</td><td>Age of the vehicle in years</td><td>×</td><td>✓</td></tr><tr><td>Vehicle Use</td><td>Main use of the vehicle</td><td>✓</td><td>×</td></tr><tr><td>Zoning Plan</td><td>Cover-specific zoning plan</td><td>✓</td><td>✓</td></tr><tr><td>Year</td><td>Calendar year</td><td>✓</td><td>✓</td></tr></tbody></table>

## 4.3 EBM specification with car insurance data

In this section, we describe the training procedure of EBM. Additional details about the functionality of the InterpretML package (Nori et al., 2019) are provided in Appendix A, along with its implementation in Python.

### 4.3.1 Loss functions

To implement the EBM method described in Section 3.2, it is necessary to define the loss functions minimized during the training phase of a boosting model. For this purpose, we use the deviance of the GLM models employed to model claim frequency and claim severity. Formally, the loss function calculated on a generic dataset $\mathcal{D}$ of size $n$ is expressed as follows

$$L(\mathbf{y}, F(\mathbf{x})) = D(\mathbf{y}, F(\mathbf{x})) = -2 \ln \left( \frac{\mathcal{L}(F(\mathbf{x}))}{\mathcal{L}(\mathbf{y})} \right). \qquad (6)$$

In practice, we refer to Henckaerts et al. (2021) and use the deviance of a Poisson GLM model as the loss function for modeling claim frequency

$$L_{\text{freq}}(\mathbf{y}, F(\mathbf{x})) = \frac{2}{n} \sum_{i=1}^{n} y_i \ln \left( \frac{y_i}{e_i F(\mathbf{x}_i)} \right) - (y_i - e_i F(\mathbf{x}_i)), \quad y_i \in \mathcal{Y}, \quad \mathbf{x}_i \in \mathcal{X}. \qquad (7)$$

where $e_i$ the exposure-to-risk for policyholder $i$. Similarly for modeling claim severity, we consider the deviance of a gamma GLM as the loss function

$$L_{\text{sev}}(\mathbf{y}, F(\mathbf{x})) = \frac{2}{n} \sum_{i=1}^{n} \alpha_i \left( \frac{y_i - F(\mathbf{x}_i)}{F(\mathbf{x}_i)} - \ln \left( \frac{y_i}{F(\mathbf{x}_i)} \right) \right), \quad y_i \in \mathcal{Y}, \quad \mathbf{x}_i \in \mathcal{X}. \qquad (8)$$

where $\alpha_i$ corresponds to a shape parameter that we ignore in practice.

The Poisson deviance (resp. the gamma deviance) can be implemented with EBM by assigning objective = poisson_deviance (resp. objective = gamma_deviance). Offsetting exposure-to-risk can be achieved during the fitting procedure with the init_score parameter, a per-sample initialization score. Exposure-to-risk is then added to the final score at prediction time. For both EBM gamma and Poisson, the automatically assigned link function $g$ is the logarithmic link function, which is often used in insurance pricing applications.

### 4.3.2 Tuning procedure

By construction, the EBM algorithm requires optimization of several hyperparameters. The general parameters, such as the loss function to minimize, the type, and the number of features, are fixed. The parameters controlling the complexity of decision trees used as building blocks as well as those for the bagging, boosting and binning procedures are tuned by grid search. Optimal hyperparameters for claim severity and frequency models are displayed in Table 2. Note that the
<table><caption>Table 2: Grid search and optimal hyperparameters for the claim severity and frequency EBM models.</caption><thead><tr><th>Hyperparameter</th><th>Grid Severity</th><th>Optimal Severity</th><th>Grid Frequency</th><th>Optimal Frequency</th></tr></thead><tbody><tr><td>smoothing_rounds</td><td>[100, 200, 500]</td><td>500</td><td>[200, 500]</td><td>500</td></tr><tr><td>interactions</td><td>[0.5, 0.75, 0.9]</td><td>0.9</td><td>[0.75, 0.9]</td><td>0.9</td></tr><tr><td>max bins</td><td>[32, 1024, 4096]</td><td>32</td><td>[32, 1024]</td><td>1024</td></tr><tr><td>max interaction bins</td><td>[8, 16, 32]</td><td>32</td><td>×</td><td>×</td></tr><tr><td>outer bags</td><td>[14, 50]</td><td>14</td><td>×</td><td>×</td></tr><tr><td>learning_rate</td><td>[0.005, 0.01, 0.02]</td><td>0.02</td><td>[0.005, 0.01]</td><td>0.005</td></tr><tr><td>max leaves</td><td>[2, 3, 4]</td><td>3</td><td>×</td><td>×</td></tr></tbody></table>

grid considered for the claim frequency model is smaller compared to the one used for the claim severity model in order to ensure reasonable computational time.

The tuned EBM performs overall slightly better than the default one and is not prone to overfitting. Quality of fit can also be assessed feature by feature through a univariate analysis of observed and predicted values. An example for the Driver Age feature is given in Figures 3a and 3b for the claim frequency GLM model, and in Figures 3c and 3d for the claim severity GLM model.

<figure><img src="image_4.png" /><figcaption>Figure 3: Observed and predicted values for the Driver Age feature in the claim severity and frequency models. The bar chart shows the exposition (left y-axis), the black line chart the observed risk and the purple line chart the predicted risk by EBM (right y-axis).</figcaption></figure>

## 4.4 Benchmark machine learning models

To assess the predictive performance of EBM, four benchmark models are constructed for both claim frequency and severity modeling. Firstly, we consider the reference model for non-life insurance pricing: the GLM model (Section 2.2.1). We also implement a GAM model (Section 2.2.2), from which EBM inherits the additive structure, and a CART model (Breiman et al., 1984), whose components are used as building blocks in EBM. Finally, a full-complexity model is implemented: the stochastic gradient boosting, called XGBoost (XGB), proposed by Chen and Guestrin (2016).

The reference GLM model is based on the Poisson distribution for the claim frequency model and the gamma distribution for the claim severity model, see Appendix C.1 for regression coefficients
of these GLMs. The use of the logarithmic link function in both claim frequency and severity modeling results in a multiplicative tariff grid that is immediately interpretable and explainable.

The competing benchmark GAM model combines cubic regression splines penalized by the square second derivative cubic spline penalty. Smoothing parameters are estimated using the Generalized Cross Validation (GCV) criterion in the mgcv R-package (Wood, 2023). In the same way as GLM, the Poisson distribution is used in the claim frequency model, and the gamma distribution in the claim severity model, both with the logarithmic link function.

The CART model is based on a two-phase approach using the rpart R-package (Therneau and Atkinson, 2022). In the first phase, deep trees are grown by neutralizing the complexity parameter (cp = 0), as well as fixing the minimum number of observations in each node (minsplit = 2) and the maximum depth of the final tree (maxdepth = 30). In the second phase, these trees are pruned using the one standard error criterion. For the claim frequency model, Poisson trees are implemented. The response information thus contains the response rate (exposition) and the number of events at the node (claim number).

Regarding the XGB, the Poisson deviance in Equation (7) is used as the loss function for the claim frequency model. The loss function for the claim severity model is the gamma deviance (8). This model is implemented using the xgboost Python-package (Chen and Guestrin, 2016). To leverage maximum performance and prevent overfitting, XGB hyperparameters are tuned by grid search as described in Table 3.

<table><caption>Table 3: Optimal hyperparameters for the claim severity and frequency XGB models in the xgboost Python-package.</caption><thead><tr><th>Parameter</th><th>Explanation</th><th>Grid</th><th>Optimal Severity</th><th>Optimal Frequency</th></tr></thead><tbody><tr><td>learning_rate</td><td>Learning rate for boosting</td><td>[0.1, 0.3, 0.5]</td><td>0.1</td><td>0.1</td></tr><tr><td>n_estimators</td><td>Number of iterations</td><td>[100, 200, 300]</td><td>200</td><td>200</td></tr><tr><td>min_child_weight</td><td>Minimum sum of weights required in a child</td><td>[1, 3]</td><td>1</td><td>1</td></tr><tr><td>max_depth</td><td>Maximum depth of a tree</td><td>[4, 6, 8]</td><td>4</td><td>4</td></tr><tr><td>gamma</td><td>Minimum loss reduction required to make a split</td><td>[0, 0.05]</td><td>0</td><td>0.05</td></tr><tr><td>alpha</td><td>L1 regularization term on weights</td><td>[0, 1]</td><td>1</td><td>1</td></tr><tr><td>subsample</td><td>Subsample ratio of the training instances</td><td>[0.75, 1]</td><td>1</td><td>×</td></tr><tr><td>colsample_bytree</td><td>Subsample ratio of columns</td><td>[0.75, 1]</td><td>0.75</td><td>×</td></tr></tbody></table>

## 4.5 Comparison of predictive performance

We now compare the predictive performance of the EBM model with the four benchmark models. Section 4.5.1 describes our evaluation framework. In Section 4.5.2, we examine and compare the performance results using aggregated metrics. Then, we analyze in Section 4.5.3 the distribution of predictions across our different models. Finally, Murphy diagrams are analyzed in Section 4.5.4 to rank the predictive performance of the models using a deviance-based score.

### 4.5.1 Evaluation framework

To assess the quality of a predictive model, performance metrics are introduced. In the regression framework, the following commonly used metrics are considered on a generic dataset $D$ of size $n$.
**Mean Squared Error (MSE)** The MSE is the arithmetic mean of the squared deviations between observations and model predictions

$$ \text{MSE}(\mathbf{y}, F(\mathbf{x})) = \frac{1}{n} \sum_{i=1}^{n} (y_i - F(\mathbf{x}_i))^2. $$

**Root Mean Squared Error (RMSE)** The RMSE is the square root of the MSE

$$ \text{RMSE}(\mathbf{y}, F(\mathbf{x})) = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - F(\mathbf{x}_i))^2}. $$

**Mean Absolute Error (MAE)** The MAE is the arithmetic mean of the absolute values of deviations

$$ \text{MAE}(\mathbf{y}, F(\mathbf{x})) = \frac{1}{n} \sum_{i=1}^{n} |y_i - F(\mathbf{x}_i)|. $$

In the RMSE, deviations are squared before being averaged, penalizing large errors. In the MAE, individual deviations are weighted equally. Both performance metrics aim to minimize the error. However, in the insurance context, these common metrics are not always suitable. Metrics based on the notion of likelihood are better adapted to the underlying distribution of the data, such as counting processes for claim frequency or right-asymmetric distributions for claim severity.

**Explained Deviance Ratio (EDR)** The EDR is a log-likelihood ratio metric that quantifies the improvement over the null model

$$ \text{EDR}(\mathbf{y}, F(\mathbf{x})) = 1 - \frac{D(\mathbf{y}, F(\mathbf{x}))}{D(\mathbf{y}, F_0(\mathbf{x}))}, $$

where $D$ refers to the deviance defined in Equation (6), and $F_0$ is the prediction rule given by the null model. The EDR increases with the proportion of variability explained by the model.

**Normalized Gini index** The Gini index is introduced by Corrado Gini (Gini, 1912) and is derived from Max Otto Lorenz's curve (Lorenz, 1905). In economics, the Lorenz curve is a graphical representation of the cumulative distribution of income or wealth, where the cumulative proportion of total income is plotted against the cumulative proportion of the population. Mathematically, the Lorenz curve can be generalized as a graphical tool for comparing the distribution of a variable with a hypothetical uniform distribution of that variable. Although the Gini index is not a strictly consistent scoring function, as discussed in Wüthrich (2023), it remains a widely used metric for evaluating and comparing model predictions.

Many formulations and variants of the Gini index have been introduced in the literature (Yitzhaki et al., 2003). In this article, we consider the normalized (or adjusted) Gini index proposed by Xin and Huang (2024). This version of the Gini index is designed to improve interpretability and comparability while limiting potential biases. The normalized Gini index is defined as the ratio of the Gini index computed when observations are sorted in ascending order of predicted values to the Gini index computed when observations are sorted in ascending order of their actual values. Formally, let $(\sigma_{(1)}, \dots, \sigma_{(n)})$ be the ordered indices of $(F(\mathbf{x}_1), \dots, F(\mathbf{x}_n))$, such that $F(\mathbf{x}_{\sigma_{(1)}}) \le F(\mathbf{x}_{\sigma_{(2)}}) \le \dots \le F(\mathbf{x}_{\sigma_{(n)}})$. Similarly, we order the observations $(y_1, \dots, y_n)$, such that $y_1 \le \dots \le y_n$. Then, the normalized Gini index is defined as

$$ \text{Gini}_{\text{norm}}(\mathbf{y}, F(\mathbf{x})) = \frac{\frac{\sum_{i=1}^{n} iy_{\sigma(i)}}{\sum_{i=1}^{n} y_i} - \sum_{i=1}^{n} \frac{n-i+1}{n}}{\frac{\sum_{i=1}^{n} iy_{(i)}}{\sum_{i=1}^{n} y_i} - \sum_{i=1}^{n} \frac{n-i+1}{n}}. $$

**Prediction structure** Similarly to the works of Denuit et al. (2021) and Holvoet et al. (2024), we analyze and compare the prediction structures of various competing models. This analysis examines the behavior of predictions from different models by utilizing histograms to identify differences in prediction patterns and potential outliers.
**Murphy diagrams** Following Ehm et al. (2016) and Holvoet et al. (2024), we implement the Murphy diagrams, which allow for a detailed comparison of the relative merits of the competing models by verifying whether the model performance rankings based on the Poisson and gamma deviance hold across various scoring functions. More precisely, for a model *F* and for all (xᵢ, yᵢ) in the testing dataset, the Murphy diagram displays the graph of points (θ, S<sub>θ</sub>(F(x), y)), where S<sub>θ</sub>(F(x), y) corresponds to the elementary scoring function for a parameter θ ∈ ℝ. This estimated quantity is defined as

$$S_{\theta}(F(\mathbf{x}), y) = \frac{1}{n} \sum_{i=1}^{n} |\theta - y_i| \mathbf{1}_{\{\min(F(\mathbf{x}_i, y_i)) \le \theta < \max(F(\mathbf{x}_i, y_i))\}}$$

According to Ehm et al. (2016), this index compares the predictions of two models. Formally, the predictions from a model *F*<sub>A</sub> are said to dominate the predictions of a model *F*<sub>B</sub>, if and only if, ∀θ ∈ ℝ, S<sub>θ</sub>(*F*<sub>A</sub>(x), y) ≤ S<sub>θ</sub>(*F*<sub>B</sub>(x), y). The dominance of one model over another, measured through the elementary scoring function, has the advantage of being preserved regardless of the loss function.

### 4.5.2 Performance results

As defined in Section 4.5.1, several performance metrics are compared (RMSE, MAE, EDR, and normalized Gini) to analyze prediction performance both in-sample (on the training sets) and out-of-sample (on the testing sets) of the competing models: GLM, GAM, CART, EBM without interactions (denoted as EBM), EBM with bivariate interactions (denoted as EBM²), and XGB. We recall that the training and testing datasets differ for the claim frequency and severity models.

Prediction performance results on the training and testing datasets are summarized in Table 4 for the claim frequency model and in Table 5 for the claim severity model. For each model, the percentage column gives the relative gain (in green) or loss (in red) over the reference GLM model on the testing dataset in a metric-wise manner, i.e. minimization for RMSE and MAE and maximization for EDR and Gini. For both claim frequency and severity data, EBM without interactions ranks between the reference GLM and the XGB model in terms of out-of-sample performance. While the performance metrics of GAM splines and EBM without interactions are similar, EBM outperforms in terms of EDR, both for claim severity and frequency, and in terms of Gini for claim severity and approaches the results obtained by XGB. Note that by adding a small number of pairwise interaction terms, EBM substantially narrows the performance gap between the inherently intelligible models (GLM, GAM splines, shallow CART) and the full-complexity model (XGB). On our testing dataset, EBM with interactions even slightly outperforms XGB in terms of Gini, both for claim severity and frequency, and in terms of EDR for claim severity.

<table><caption>Table 4: In-sample and out-of-sample prediction performance in terms of RMSE and MAE for the claim frequency models.</caption><thead><tr><th>Metric</th><th colspan="3">RMSE</th><th colspan="3">MAE</th></tr><tr><th>Model</th><th>Train</th><th>Test</th><th>%</th><th>Train</th><th>Test</th><th>%</th></tr></thead><tbody><tr><td><strong>GLM</strong></td><td>0.1381</td><td>0.1380</td><td>0</td><td>0.0364</td><td>0.0363</td><td>0</td></tr><tr><td><strong>GAM</strong></td><td>0.1380</td><td>0.1380</td><td>0</td><td>0.0364</td><td>0.0363</td><td>-0.1%</td></tr><tr><td><strong>CART</strong></td><td>0.1380</td><td>0.1380</td><td>0</td><td>0.0363</td><td>0.0363</td><td>-0.1%</td></tr><tr><td><strong>EBM</strong></td><td>0.1380</td><td>0.1380</td><td>0</td><td>0.0363</td><td>0.0363</td><td>-0.1%</td></tr><tr><td><strong>EBM<sup>2</sup></strong></td><td>0.1380</td><td>0.1379</td><td>0</td><td>0.0363</td><td><strong>0.0363</strong></td><td>-0.1%</td></tr><tr><td><strong>XGB</strong></td><td>0.1380</td><td><strong>0.1379</strong></td><td>0</td><td>0.0364</td><td>0.0363</td><td>0</td></tr></tbody></table><table><thead><tr><th>Metric</th><th colspan="3">EDR Poisson</th><th colspan="3">Gini<sub>norm</sub></th></tr><tr><th>Model</th><th>Train</th><th>Test</th><th>%</th><th>Train</th><th>Test</th><th>%</th></tr></thead><tbody><tr><td><strong>GLM</strong></td><td>0.0089</td><td>0.0081</td><td>0</td><td>0.2963</td><td>0.3006</td><td>0</td></tr><tr><td><strong>GAM</strong></td><td>0.0119</td><td>0.0110</td><td>+35.4%</td><td>0.3068</td><td>0.3101</td><td>+3.2%</td></tr><tr><td><strong>CART</strong></td><td>0.0122</td><td>0.0100</td><td>+23.0%</td><td>0.3077</td><td>0.3071</td><td>+2.2%</td></tr><tr><td><strong>EBM</strong></td><td>0.0123</td><td>0.0110</td><td>+36.2%</td><td>0.3078</td><td>0.3102</td><td>+3.2%</td></tr><tr><td><strong>EBM<sup>2</sup></strong></td><td>0.0130</td><td>0.0117</td><td>+43.9%</td><td>0.3097</td><td><strong>0.3119</strong></td><td>+3.8%</td></tr><tr><td><strong>XGB</strong></td><td>0.0139</td><td><strong>0.0117</strong></td><td>+44.9%</td><td>0.3110</td><td>0.3111</td><td>+3.5%</td></tr></tbody></table>
<table><caption>Table 5: In-sample and out-of-sample prediction performance in terms of RMSE, MAE, EDR and normalized Gini for the claim severity models.</caption><thead><tr><th>Metric</th><th colspan="3">RMSE</th><th colspan="3">MAE</th></tr><tr><th>Model</th><th>Train</th><th>Test</th><th>%</th><th>Train</th><th>Test</th><th>%</th></tr></thead><tbody><tr><td>GLM</td><td>2112</td><td>2019</td><td>0</td><td>1326</td><td>1317</td><td>0</td></tr><tr><td>GAM</td><td>2099</td><td>2012</td><td>-0.4%</td><td>1321</td><td>1316</td><td>-0.1%</td></tr><tr><td>CART</td><td>2095</td><td>2043</td><td>+1.2%</td><td>1323</td><td>1322</td><td>+0.3%</td></tr><tr><td>EBM</td><td>2100</td><td>2010</td><td>-0.5%</td><td>1320</td><td>1314</td><td>-0.3%</td></tr><tr><td>EBM<sup>2</sup></td><td>2091</td><td><strong>2004</strong></td><td>-0.8%</td><td>1315</td><td>1308</td><td>-0.7%</td></tr><tr><td>XGB</td><td>2072</td><td>2006</td><td>-0.7%</td><td>1300</td><td><strong>1303</strong></td><td>-1.0%</td></tr></tbody></table><table><thead><tr><th>Metric</th><th colspan="3">EDR Gamma</th><th colspan="3">Gini<sub>norm</sub></th></tr><tr><th>Model</th><th>Train</th><th>Test</th><th>%</th><th>Train</th><th>Test</th><th>%</th></tr></thead><tbody><tr><td>GLM</td><td>0.0730</td><td>0.0759</td><td>0</td><td>0.3044</td><td>0.3087</td><td>0</td></tr><tr><td>GAM</td><td>0.0836</td><td>0.0827</td><td>+8.9%</td><td>0.3207</td><td>0.3214</td><td>+4.1%</td></tr><tr><td>CART</td><td>0.0674</td><td>0.0634</td><td>-16.5%</td><td>0.2730</td><td>0.2757</td><td>-10.7%</td></tr><tr><td>EBM</td><td>0.0840</td><td>0.0838</td><td>+10.5%</td><td>0.3218</td><td>0.3233</td><td>+4.7%</td></tr><tr><td>EBM<sup>2</sup></td><td>0.0898</td><td><strong>0.0885</strong></td><td>+16.7%</td><td>0.3307</td><td><strong>0.3308</strong></td><td>+7.1%</td></tr><tr><td>XGB</td><td>0.1035</td><td>0.0867</td><td>+14.3%</td><td>0.3516</td><td>0.3271</td><td>+6.0%</td></tr></tbody></table>

### 4.5.3 Prediction structure

Now, we compare the prediction structures of the competing models. Figures 4 and 5 display histograms of the predictions on the out-of-sample testing sets from claim frequency and severity datasets. We observe that the predictive distributions of GAM splines, EBM and XGB show a similar dispersion for both the claim severity and frequency models. In contrast, the predictions from the GLM are more concentrated. For the claim severity model, the predictions from CART exhibit relatively different behavior compared to the other models.

<figure><img src="image_5.png" /><figcaption>Figure 4: Histograms of the predictions in the claim frequency models on the testing dataset. All the competing models are represented.</figcaption></figure>

In addition, we conduct a univariate analysis of observed and predicted values and provide a visual inspection of each feature. Figure 6a (resp. Figure 6b) gives an example for the Driver Age feature in the claim frequency model (resp. claim severity). For both claim frequency and severity, the prediction curves of GAM splines, EBM, EBM² and XGB appear to be close and reflect the observed risk relatively well. In contrast, the GLM prediction curve shows a less complex and less nuanced relationship between Driver Age and claim severity and frequency responses.
<figure><img src="image_6.png" /><figcaption>Figure 5: Histograms of the predictions in the claim severity models on the testing dataset. All the competing models are represented.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Observed and predicted values for the Driver Age feature on the testing dataset. The bar chart shows the exposition (left y-axis), the black line chart the observed claim frequency and the colored line charts the competing benchmark models (right y-axis).</figcaption></figure>

### 4.5.4 Predictive dominance analysis

In this section, we compare the EBM and all competing benchmark models in terms of predictive dominance based on Murphy diagrams. Figure 7 shows the Murphy diagram for both claim frequency and severity models. These diagrams are built on testing datasets, and we focus on the range of predicted values where the models show significant differences. For the claim severity model, EBM² with interactions appears to have the lowest value of $S_θ$ for all $θ$, indicating its predictive dominance over the other competing benchmark models. The XGB, EBM² without interactions, and GAM models also seem to dominate the GLM and CART models. However, in the claim frequency case, none of the models is consistently lower than the others. The values of XGB and EBM² with interactions are very close. For the claim severity, the values of CART and GLM are significantly higher than those of the other models, but they are not consistently higher relative to each other.

## 5 Feature interpretation with EBM

After assessing the predictive performance of EBM compared to other standard machine learning models, we now examine how this model provides directly interpretable analyses of its predictions. In Section 5.1, we describe the various interpretability tools offered by the EBM model. To evaluate the robustness of the insights provided, we compare them in Section 5.2 with explanations derived from post-hoc interpretation tools applied to competing models. This approach enables us to assess
<figure><img src="image_8.png" /><figcaption>Figure 7: Murphy diagrams on the testing set for all competing claim frequency and severity models according to the elementary scoring function $S_{\theta}$ over a range of values for the parameter $\theta$. A zoom on the areas of interest is presented.</figcaption></figure>

the relevance of the insights that can be drawn from the EBM model. For the sake of conciseness, we present only the results of the EBM model with interactions and the XGB model.

## 5.1 EBM interpretability

Since the structure of EBM is additive, each feature contributes to the prediction in a modular manner. Consequently, the main effect terms and interaction terms provide a clear description of how each feature influences the prediction. In other words, the final effect can be visualized and understood by analyzing the shape function $f_j$ (or $f_{(j_1, j_2)}$ for pairwise interactions) in a trained model, as defined in Equation (4). To this end, we define interpretability tools based on shape functions for a generic dataset $\mathcal{D}$ of size $n$. All our results are then presented on the training set.

### 5.1.1 Interpretation with shape function

For each $\mathbf{x}_i \in \mathcal{X}$, the shape function for variable $X_j$ $f_j(\mathbf{x}_i)$ with $j \in \{1, \dots, p\}$ locally measures its effect for individual $i$. Similarly to the regression coefficients in a GLM model, the values of the average shape functions can be compared to assess the importance of each feature. Note that these shape functions are centered and normalized to unit variance. In addition, the use of the logarithmic link function in the claim frequency and severity models means that the y-axis of shape functions as well as the scale of heatmaps presented below are in logarithmic space. In this section, we present only the results from the claim severity model. A similar analysis for the claim frequency model is provided in Appendix C.2.

Figure 8 provides a comprehensive visualization of the complete set of features in the claim severity model: Driver Age, License Seniority, Secondary Novice Driver, CRM Coefficient, Vehicle Price Class, Vehicle Age, Zoning Plan, and Year. For each of these features $X_j$, the plot is divided into two sections: the upper part displays the estimated shape function, illustrating the feature's contribution to the prediction (*Score*), while the lower part presents a histogram depicting the feature's distribution (*Exposure*) in the training dataset. Uncertainty in the shape functions is represented by error bars, which account for both the uncertainty due to the available data volume in different regions of the feature space and the uncertainty of the learned model. These errors are estimated using a bagging procedure, where the EBM model is trained on different randomly selected subsets of the data. This process is controlled by the *validation_size* and *outer_bags* parameters, as described in Table 7, which outlines the bagging process in the EBM algorithm. For continuous features, a step function is used to enhance result visualization. To achieve this, the values of variable $X_j$ are segmented, and the predicted shape function value is reported for the median of each segment. The segment names have been anonymized for the variables Driver Age and License Seniority. The number of steps is at most equal to the number of *bins* in the binning procedure. For categorical features, a bar chart is displayed for each component.

We clearly observe that the variables Vehicle Price Class, Year (capturing an inflation effect), and CRM Coefficient are the most influential. They exhibit a generally increasing effect on the response variable. Similarly, the least influential variable, Zoning Plan, also shows an increasing trend. In contrast, the variables Driver Age, License Seniority and Vehicle Age generally have a decreasing effect on claim severity for younger ages and lower seniority levels, but their patterns
<figure><img src="image_9.png" /><figcaption>Figure 8: Prediction of the estimated shape functions $f_j(\mathbf{x}_i)$ in the claim severity model for features Driver Age, License Seniority, Secondary Novice Driver, CRM Coefficient, Vehicle Price Class, Vehicle Age, Zoning Plan and Year. The shaded areas displays an error bar calculated using a bagging procedure. We use the default `validation_size` and `outer_bags` parameters, as defined by the InterpretML algorithm described in Appendix A.</figcaption></figure>

become more complex for higher values. For instance, the effect of Driver Age stabilizes for ages between 40 and 50, as well as between 60 and 70, before increasing for older ages, which is generally consistent with observed claim frequency patterns in car insurance. Finally, the presence of a Secondary Novice Driver increases claim severity. Next, we examine the effect of interactions between variables. Heatmaps of the selected pairwise interaction terms chosen by the GA²Ms algorithm can be found in Figure 9 for the claim severity model. Recall that the interactions selected by the FAST algorithm in EBM as the most effective in the claim severity model primarily combine policyholder characteristics with vehicle characteristics. Several notable interactions emerge. First, the interaction between Driver Age and Vehicle Price Class shows an increased severity risk for very young drivers with expensive vehicles, similar to the interaction between License Seniority and
Vehicle Price Class, which displays a more widespread effect across driving experience. In contrast, very experienced drivers with either older or very expensive vehicles exhibit lower claim severity. The interaction between Driver Age and Vehicle Age reveals an increase in severity for very young drivers with very old vehicles, while showing a decrease for older drivers with older vehicles. Next, the interaction between the CRM Coefficient and Vehicle Age indicates that vehicles with a high CRM Coefficient, i.e. those that have experienced a higher claim frequency in the past, and low vehicle age tend to have higher claim severity. However, beyond a certain Vehicle Age, severity seems to decrease with the CRM Coefficient. Finally, we observe that the Zoning Plan variable interacts with the CRM Coefficient, Driver Age, and License Seniority. However, the coding of the Zoning Plan variable makes it difficult to interpret these relationships clearly.

### 5.1.2 Feature importance score

Now, we evaluate the global feature importance score of variable $X_j$ based on the estimated shape function $f_j$ (or $f_{(j_1, j_2)}$) as follows

$$I_j = \frac{\sum_{k=1}^{K} w_{j,k} |\sum_{i=1}^{n} f_{j,k}(\mathbf{x}_i) \mathbf{1}_{\{x_{i,j} \in b_{j,k}\}}|}{\sum_{k=1}^{K} w_{j,k}}, \quad (9)$$

where $w_{j,k} = \sum_{i=1}^{n} \mathbf{1}_{\{x_{i,j} \in b_{j,k}\}}$ corresponds to the weight of bin $b_{j,k}$, i.e. the number of observations in that bin. A high value of the index $I_j$ indicates that the variable $X_j$ contributes significantly to the model's predictions. Then, these indices can be compared to rank the importance of different features.

In Figure 10, we present, both for univariate and pairwise interaction terms, the relative feature importance for claim frequency and severity models for variable $X_j$.

$$RI_j = \frac{I_j}{\sum_{j=1}^{p} I_j}.$$

A feature characterizing the driver, Driver Age, is the most important feature in the claim frequency model, while a feature characterizing the vehicle, Vehicle Price Class, is the dominant feature in the claim severity model. It is also worth noting that the variables CRM Coefficient and Year play a significant role. Finally, the pairwise interaction terms have a smaller contribution than the univariate variables.

## 5.2 Comparison of EBM and XGBoost interpretability

The strength of the EBM model lies in its ability to be directly interpretable, whether for measuring the effect of a single feature or analyzing pairwise interaction terms. In Section 5.1, we highlighted its usefulness in transparently examining the relationships between features and the response variable. In this section, we compare the EBM model with a black-box model, namely the XGB model introduced in Section 4.4, to extract insights from the data. The XGB model, widely used for pricing insurance products, provides predictive performance comparable to that of the EBM model. However, additional analyses are required to interpret its predictions. In our analysis, we apply several post-hoc interpretation techniques to ensure a fair comparison of both methods. This allows us to demonstrate the advantage of the EBM model over the XGB model, as it inherently provides consistent and interpretable analyses without relying on post-hoc interpretation techniques.

### 5.2.1 Comparison of importance measure

First, we introduce a post-hoc feature importance score, as defined, for example, by Delcaillau et al. (2022) or Charpentier (2024, Chapter 4). This indicator is different from the one introduced in Equation 9. It is based on the permutation method and serves as a global agnostic interpretation tool that helps assess the contribution of each feature to the model's predictions. Formally, this indicator for variable $X_j$, with $j \in \{1, \dots, p\}$, is defined as on a generic dataset $\mathcal{D}$ of size $n$

$$FI_j = \frac{1}{n} \sum_{i=1}^{n} L(y_i, F(\mathbf{x}_{i,-j}, x_{i,j})) - L(y_i, F(\mathbf{x}_{i,-j}, \tilde{x}_{i,j})), \quad (10)$$

where $\mathbf{x}_i = (\mathbf{x}_{i,-j}, x_{i,j})$ and $\tilde{x}_{i,j}$ is a random permutation of $x_j$ observed on $\mathcal{D}$. A large value of $FI_j$ indicates that $X_j$ significantly influences the model output and is therefore considered important.
<figure><img src="image_10.png" /><figcaption>Figure 9: Heatmaps of the predicted pairwise interaction terms $f_{(j_1, j_2)}$ in the claim severity model. The pairs of variables include: (CRM Coefficient, Vehicle Age), (CRM Coefficient, Zoning Plan), (Driver Age, Vehicle Age), (Driver Age, Vehicle Price Class), (Driver Age, Zoning Plan), (License Seniority, Vehicle Age), (License Seniority, Vehicle Price Class) and (License Seniority, Zoning Plan).</figcaption></figure>

Figure 11 compares the relative importance of each variable for the EBM and XGB models for claim frequency and severity data using the relative feature importance metric

$$RFI_j = \frac{FI_j}{\sum_{j=1}^{p} FI_j}.$$

Most features exhibit a similar importance ranking in both EBM and XGB, particularly in the claim severity model, where only an inversion of importance is observed between the features Driver Age
<figure><img src="image_11.png" /><figcaption>Figure 10: Features importance for univariate and pairwise interaction terms for the claim frequency and severity models.</figcaption></figure>

and Year (Figure 11b). Some variations appear in the claim frequency model, but Driver Age, CRM Coefficient, and Year remain the dominant features in both EBM and XGB (Figure 11a). These results suggest that EBM and XGB capture comparable relationships between the target variable and the features.

For the main effect terms in EBM, we observe that the variable importance ranking derived from the permutation method aligns with the ranking based on the shape values in Equation (9), see Figure 10. However, we note significant differences between the two feature importance scores regarding the relative impact of the dominant variables in both the claim frequency and the claim severity models. In particular, it is interesting to observe that the results obtained using Equation (9) are more evenly distributed across variables compared to those computed with Equation (10), where the two most important features dominate the others. These discrepancies could partly be explained by the correlation of these variables with other features, as well as the complex effects of Driver Age, Vehicle Price Class and CRM Coefficient on the response variable.

<figure><img src="image_12.png" /><figcaption>Figure 11: Relative feature importance in EBM and XGB using the permutation method. Features are sorted by descending variable importance in EBM.</figcaption></figure>

### 5.2.2 Partial dependence analysis

Partial dependence analysis is a global agnostic interpretation method designed to assess the average marginal effect of one or more features on model predictions (Delcaillau et al., 2022; Charpentier, 2024). Let $\mathcal{X}_j$ be the vector containing all possible values of the variable $X_j$. For a given model $F$,
the partial dependence function for variable $X_j$ is estimated as follows

$$PD(x_j^*) = \frac{1}{n} \sum_{i=1}^{n} F(x_j^*, \mathbf{x}_{i,-j}), \quad \forall x_j^* \in \mathcal{X}_j, \qquad (11)$$

where the element $x_{i,j}$ in $\mathbf{x}_i$ is replaced by a value $x_j^* \in \mathcal{X}_j$. The partial dependence plot (PDP) represents the value of $PD(x_j^*)$ for each value of $x_j^*$. It corresponds to the average effect across all individuals, in contrast to the Individual Conditional Expectation (ICE) plot, which displays the dependence of the prediction for each individual, i.e. $F(x_j^*, \mathbf{x}_{i,-j})$ (Friedman, 2001). The analysis of PDPs can be compared with the analysis of shape functions produced by the EBM model to verify the consistency between both approaches. A drawback of PDPs lies in the strong assumption that $X_j$ is uncorrelated with the remaining features.

Figure 12 illustrates the partial dependence effect plots for Driver Age and Vehicle Price Class in both the EBM and XGB models trained on the claim severity data. For the Driver Age variable (Figure 12b), the partial dependence effect decreases until around age 40, stabilizes between ages 40 and 50, decreases again, and then remains stable between ages 60 and 70. In the EBM model, the effect increases for older ages, whereas in the XGB model, it remains slightly decreasing. Despite these subtle differences, both models exhibit a comparable partial dependence effect. The analysis of ICE curves reveals that the effect of Driver Age is relatively consistent across individuals and aligns with the PDP results. However, there is significant dispersion in this effect among individuals. For the Vehicle Price Class variable (Figure 12b), the PDP indicates that both models capture an increasing relationship between Vehicle Price Class and predicted claim severity. However, in the XGB model, the increase is steeper for Vehicle Price Class values above 17. The ICE curves of the XGB model suggest that beyond a value of 17 the effect of Vehicle Price Class becomes more heterogeneous across individuals, indicating potential interactions with other variables. These interactions are not as apparent in the EBM model.

Now, let us compare these results with those of the shape functions in Figure 8a and Figure 8e for the EBM model. While we observe similar curves for the variable Driver Age, the shape functions highlight the relationship between severity and Driver Age more distinctly at older ages. Regarding the variable Vehicle Price Class, both the shape function and PDP curves exhibit a similar increasing trend, with the shape function capturing a stronger sensitivity for higher price classes. These analyses support the credibility of using shape functions to visualize the impact of a variable on the response variable. Similar results are presented in Appendix C.3 for the other features used in the claim severity model.

<figure><img src="image_13.png" /><figcaption>Figure 12: Partial Dependence Plot (PDP) in EBM and XGB for the claim severity model. We also show the Individual conditional expectation (ICE) curves for some observations.</figcaption></figure>

### 5.2.3 Comparison of profiles using local interpretability analysis

The interpretability of the EBM model can also be assessed locally by directly analyzing the shape functions. To illustrate the relevance of this approach, we identify, similarly to Holvoet et al. (2024), three risk profiles in our data: low, medium, and high, as defined in Table 6.

For each profile, we analyze in Figure 13 the contribution of individual features to the model's prediction in the EBM model based on the estimated shape functions. A positive score computed
<table><caption>Table 6: Example of a low, medium and high risk profile in the claim frequency model and the prediction obtained by EBM and XGB.</caption><thead><tr><th>Feature</th><th>Low Risk</th><th>Medium Risk</th><th>High Risk</th></tr></thead><tbody><tr><td>Driver age</td><td>50</td><td>29</td><td>82</td></tr><tr><td>License seniority</td><td>23</td><td>11</td><td>62</td></tr><tr><td>Secondary novice driver</td><td>0</td><td>0</td><td>0</td></tr><tr><td>CRM coefficient</td><td>0.5</td><td>0.6</td><td>1</td></tr><tr><td>Vehicle price class</td><td>1</td><td>10</td><td>15</td></tr><tr><td>Vehicle use</td><td>personal</td><td>professional</td><td>personal</td></tr><tr><td>Zoning plan</td><td>4</td><td>5</td><td>15</td></tr><tr><td>Year</td><td>2020</td><td>2020</td><td>2013</td></tr><tr><th colspan="4">Prediction</th></tr><tr><td>EBM</td><td>0.0081</td><td>0.0330</td><td>0.1701</td></tr><tr><td>XGB</td><td>0.0091</td><td>0.0354</td><td>0.1282</td></tr></tbody></table>

from a shape function indicates a higher contribution of the variable relative to the average predicted response, whereas a negative score corresponds to a lower contribution of this variable. By construction, the EBM model can also capture the contribution of pairwise interaction terms. The low-risk profile corresponds to an experienced and cautious driver (CRM Coefficient equal to 0.5) who uses his vehicle for personal purposes and does not share it with a novice secondary driver. These characteristics, combined with a low Vehicle Price Class, a low-risk Zoning Plan and the Year 2020 (marked by the COVID-19 pandemic), contribute negatively to the prediction. The License Seniority for this profile is relatively low compared to the driver's age and has a slight positive contribution to the prediction (Figure 13a). In the medium-risk profile, we observe a 29-year-old driver using his vehicle for professional purposes, with a slightly deteriorated CRM Coefficient, a higher Vehicle Price Class and lower License Seniority. These features have a positive contribution to the prediction, while the other features contribute mostly negatively (Figure 13b). Finally, the high-risk profile describes a Senior Driver with a significantly deteriorated CRM Coefficient, a high Vehicle Price Class, and a high-risk Zoning Plan, leading to a particularly high predicted risk (Figure 13c).

<figure><img src="image_14.png" /><figcaption>Figure 13: Local explanations in EBM for the low, medium and high risk profile in the claim frequency model. The contributions of the different variables are estimated based on the shape functions.</figcaption></figure>

We now compare this analysis with the local interpretations obtained from the XGB model. For this purpose, we use the SHAP method, a local agnostic interpretation tool based on the game-theoretically optimal Shapley values (Shapley et al., 1953), see Charpentier (2024, Chapter 4) for an introduction to this method. The interpretation of Shapley values follows the same logic as that of shape functions: a positive (respectively negative) value corresponds to a higher (respectively lower) prediction compared to the average prediction. For the specific profiles presented in Table 6, EBM and XGB exhibit a very similar understanding of risk regarding univariate variables. In the low-risk profile, all features contribute negatively to the prediction, with the Vehicle Price Class being the most influential factor (Figure 14a). In the medium-risk profile, the CRM Coefficient and Vehicle Use make important positive contributions, whereas Year, Driver Age and Zoning Plan have negative contributions (Figure 14b). Finally, in the high-risk profile, the strongly positive contributions of Driver Age, CRM Coefficient, Year and Vehicle Price Class result in a high predicted risk (Figure 14c). For this profile, the License Seniority variable is interpreted differently between the EBM and XGB models. Indeed, XGB considers a positive contribution to the prediction,
whereas EBM assigns a negative contribution, which is partially offset by a positive contribution when interacting with the Vehicle Price Class variable.

<figure><img src="image_15.png" /><figcaption>Figure 14: Shapley values of the XGB model for the low, medium and high risk profile in the claim frequency model.</figcaption></figure>

The comparison of different post-hoc interpretability techniques between EBM and XGB consistently shows that both models capture similar relationships between the target variable and the features. The key difference could lie in the nature of their interpretability: EBM's global explanations provide an exact and complete representation of the model itself, whereas XGB relies on post-hoc interpretability methods, which approximate the model's behavior without being inherently transparent.

# 6 Conclusion

In this paper, we explore the potential of using the Explainable Boosting Machine (EBM) model for analyzing claim frequency and severity in non-life insurance, which are the two key components in determining technical premiums. We develop an analytical framework applied to car insurance data, allowing us to compare the predictive performance of EBM and to highlight its advantages in terms of interpretability. The EBM model is a glass-box model based on a GAM model, where shape functions are independently trained using a cyclic gradient boosting trees procedure. Due to this additive structure, the contribution of each feature can be directly assessed through its impact on the score function. Furthermore, this method incorporates a GA²M selection algorithm to identify the most relevant pairwise interaction terms.

Our empirical results demonstrate the good out-of-sample performance of the EBM model, particularly when pairwise interaction terms are included. We compare its predictive accuracy with state-of-the-art regression models commonly used in insurance, such as the XGB model. Additionally, we assess the relevance of EBM's predictions through various interpretability tools. This approach enables both local and global consistency checks between post-hoc interpretability analyses, typically conducted by practitioners on black-box models like XGB, and direct interpretation of EBM predictions. Building such a framework is particularly valuable in insurance as it ensures high transparency in prediction interpretation and facilitates variable selection, including the identification of meaningful variable interactions. Unlike black-box models such as XGB or deep learning approaches, directly interpretable models like EBM reduce the need for developing complex surrogate white-box models for pricing insurance contracts. Furthermore, EBM can serve as a powerful tool for guiding the development of high-performing GLMs, either by aiding variable selection or optimizing feature engineering.

Our framework opens several avenues for future research. First, EBM training relies on deviance as the loss function. Exploring alternative loss functions specifically tailored to insurance pricing could further enhance its applicability. Another promising extension would be to develop dimension-wise early stopping, similarly to Zakrisson and Lindholm (2025). Finally, it would be interesting to investigate the relevance of this model for more complex datasets, such as those with spatial interactions or temporal components, as well as other risks, such as home insurance or health insurance.

**Acknowledgments** Quentin Guibert gratefully acknowledges funding from the Institut Europlace de Finance (IEF) as part of the research project "Machine Learning and Explainability".
**Supplementary material** The results in this paper are obtained using R and Python. Supplementary material related to this paper can be found at https://github.com/MarketaKrupova/EBM-Car-Insurance.

**Declarations** The authors declare no potential conflict of interests.

**Open Access** Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Agarwal, R., Melnick, L., Frosst, N., Zhang, X., Lengerich, B., Caruana, R., and Hinton, G. (2021). Neural Additive Models: Interpretable Machine Learning with Neural Nets. *Advances in neural information processing systems* 34, pp. 4699–4711.

Arik, S. O. and Pfister, T. (2021). TabNet: Attentive Interpretable Tabular Learning. *Proceedings of the AAAI Conference on Artificial Intelligence* 35.8, pp. 6679–6687. DOI: 10.1609/aaai.v35i8.16826.

Bergstra, J. and Bengio, Y. (2012). Random search for hyper-parameter optimization. *Journal of machine learning research* 13.2.

Blier-Wong, C., Cossette, H., Lamontagne, L., and Marceau, E. (2021). Machine Learning in P&;C Insurance: A Review for Pricing and Reserving. *Risks* 9.1, p. 4. DOI: 10.3390/risks9010004.

Breiman, L. (2001). Random Forests. *Machine Learning* 45.1, pp. 5–32. DOI: 10.1023/A:1010933404324.

Breiman, L., Friedman, J., Stone, C. J., and Olshen, R. A. (1984). Classification and Regression Trees. New Ed. Boca Raton: Chapman and Hall/CRC.

Caruana, R., Lou, Y., Gehrke, J., Koch, P., Sturm, M., and Elhadad, N. (2015). Intelligible Models for HealthCare: Predicting Pneumonia Risk and Hospital 30-day Readmission. *Proceedings of the 21st ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*. KDD ’15. New York, NY, USA: Association for Computing Machinery, pp. 1721–1730. DOI: 10.1145/2783258.2788613.

Charpentier, A. (2024). Insurance, Biases, Discrimination and Fairness. Springer Actuarial. Cham: Springer Nature Switzerland. DOI: 10.1007/978-3-031-49783-4.

Chen, T. and Guestrin, C. (2016). XGBoost: A Scalable Tree Boosting System. *Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*. Vol. 11. KDD ’16. ACM, 785–794. DOI: 10.1145/2939672.2939785.

Delcaillaud, D., Ly, A., Papp, A., and Vermet, F. (2022). Model transparency and interpretability: survey and application to the insurance industry. *European Actuarial Journal* 12.2, pp. 443–484. DOI: 10.1007/s13385-022-00328-y.

Delong, L. and Kozak, A. (2023). The use of autoencoders for training neural networks with mixed categorical and numerical features. *ASTIN Bulletin: The Journal of the IAA* 53.2, pp. 213–232. DOI: 10.1017/asb.2023.15.

Delong, L., Lindholm, M., and Zakrisson, H. (2023). On Cyclic Gradient Boosting Machines. SSRN Scholarly Paper. Rochester, NY. DOI: 10.2139/ssrn.4352505.

Denuit, M., Charpentier, A., and Trufin, J. (2021). Autocalibration and Tweedie-dominance for insurance pricing with machine learning. *Insurance: Mathematics and Economics* 101, pp. 485–497. DOI: 10.1016/j.insmatheco.2021.09.001.
Denuit, M., Hainaut, D., and Trufin, J. (2019). Effective Statistical Learning Methods for Actuaries I: GLMs and Extensions. Springer Actuarial. Cham: Springer International Publishing. DOI: 10.1007/978-3-030-25820-7.

Denuit, M. and Lang, S. (2004). Non-life rate-making with Bayesian GAMs. *Insurance: Mathematics and Economics* 35.3, pp. 627–647. DOI: 10.1016/j.insmatheco.2004.08.001.

Denuit, M. and Legrand, C. (2018). Risk classification in life and health insurance: extension to continuous covariates. *European Actuarial Journal* 8.1, pp. 245–255. DOI: 10.1007/s13385-018-0171-9.

Ehm, W., Gneiting, T., Jordan, A., and Krüger, F. (2016). Of quantiles and expectiles: consistent scoring functions, Choquet representations and forecast rankings. *Journal of the Royal Statistical Society Series B: Statistical Methodology* 78.3, pp. 505–562.

Eling, M., Nuessle, D., and Staubli, J. (2022). The impact of artificial intelligence along the insurance value chain and on the insurability of risks. *The Geneva Papers on Risk and Insurance - Issues and Practice* 47.2, pp. 205–241. DOI: 10.1057/s41288-020-00201-7.

Embrechts, P. and Wüthrich, M. V. (2022). Recent Challenges in Actuarial Science. *Annual Review of Statistics and Its Application* 9. Volume 9, 2022, pp. 119–140. DOI: 10.1146/annurev-statistics-040120-030244.

Frees, E. W. J. and Huang, F. (2023). The Discriminating (Pricing) Actuary. *North American Actuarial Journal* 27.1, pp. 2–24. DOI: 10.1080/10920277.2021.1951296.

Friedman, J. H. (2001). Greedy function approximation: A gradient boosting machine. *The Annals of Statistics* 29.5, pp. 1189–1232. DOI: 10.1214/aos/1013203451.

Gini, C. (1912). Variabilità e mutabilità: contributo allo studio delle distribuzioni e delle relazioni statistiche. [Fasc. I.] Tipogr. di P. Cuppini.

Hastie, T. J. and Tibshirani, R. J. (1990). Generalized Additive Models. Boca Raton, Fla: Chapman and Hall/CRC.

Henckaerts, R., Antonio, K., Clijsters, M., and Verbelen, R. (2018). A data driven binning strategy for the construction of insurance tariff classes. *Scandinavian Actuarial Journal* 2018.8, pp. 681–705. DOI: 10.1080/03461238.2018.1429300.

Henckaerts, R., Antonio, K., and Côté, M.-P. (2022). When stakes are high: Balancing accuracy and transparency with Model-Agnostic Interpretable Data-driven suRRogates. *Expert Systems with Applications* 202, p. 117230.

Henckaerts, R., Côté, M.-P., Antonio, K., and Verbelen, R. (2021). Boosting Insights in Insurance Tariff Plans with Tree-Based Machine Learning Methods. *North American Actuarial Journal* 25.2, pp. 255–285. DOI: 10.1080/10920277.2020.1745656.

Hofner, B., Mayr, A., and Schmid, M. (2016). gamboostLSS: An R Package for Model Building and Variable Selection in the GAMLSS Framework. *Journal of Statistical Software* 74, pp. 1–31. DOI: 10.18637/jss.v074.i01.

Holvoet, F., Antonio, K., and Henckaerts, R. (2024). Neural networks for insurance pricing with frequency and severity data: a benchmark study from data preprocessing to technical tariff. DOI: 10.48550/arXiv.2310.12671.

Krúpová, M. (2023). Construction d'un modèle de Machine Learning interprétable pour la tarification en assurance non-vie. MA thesis. Paris: Université Paris Dauphine - PSL.

Lindholm, M., Richman, R., Tsanakas, A., and Wüthrich, M. V. (2022). Discrimination-Free Insurance Pricing. *ASTIN Bulletin: The Journal of the IAA* 52.1, pp. 55–89. DOI: 10.1017/asb.2021.23.

Lorentzen, C. and Mayer, M. (2020). Peeking into the Black Box: An Actuarial Case Study for Interpretable Machine Learning. SSRN Scholarly Paper ID 3595944. Rochester, NY: Social Science Research Network. DOI: 10.2139/ssrn.3595944.

Lorenz, M. O. (1905). Methods of measuring the concentration of wealth. *Publications of the American statistical association* 9.70, pp. 209–219.
Lou, Y., Caruana, R., and Gehrke, J. (2012). Intelligible models for classification and regression. *Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining*. Beijing China: ACM, pp. 150–158. DOI: 10.1145/2339530.2339556.

Lou, Y., Caruana, R., Gehrke, J., and Hooker, G. (2013). Accurate intelligible models with pairwise interactions. *Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining*. KDD ’13. New York, NY, USA: Association for Computing Machinery, pp. 623–631. DOI: 10.1145/2487575.2487579.

Luenberger, D. G. and Ye, Y. (2021). Linear and Nonlinear Programming. Vol. 228. International Series in Operations Research & Management Science. Cham: Springer International Publishing. DOI: 10.1007/978-3-030-85450-8.

Maillart, A. (2021). Toward an explainable machine learning model for claim frequency: a use case in car insurance pricing with telematics data. *European Actuarial Journal* 11.2, pp. 579–617. DOI: 10.1007/s13385-021-00270-5.

Maillart, A. and Robert, C. (2024). Distill knowledge of additive tree models into generalized linear models: a new learning approach for non-smooth generalized additive models. *Annals of Actuarial Science* 18.3, pp. 692–711. DOI: 10.1017/S1748499524000241.

Mayr, A., Fenske, N., Hofner, B., Kneib, T., and Schmid, M. (2012). Generalized Additive Models for Location, Scale and Shape for High Dimensional Data—A Flexible Approach Based on Boosting. *Journal of the Royal Statistical Society Series C: Applied Statistics* 61.3, pp. 403–427. DOI: 10.1111/j.1467-9876.2011.01033.x.

McCullagh, P. and Nelder, J. A. (1989). Generalized Linear Models. 2nd. Chapman and Hall.

McDonnell, K., Murphy, F., Sheehan, B., Masello, L., and Castignani, G. (2023). Deep learning in insurance: Accuracy and model interpretability using TabNet. *Expert Systems with Applications* 217, p. 119543. DOI: 10.1016/j.eswa.2023.119543.

Meyers, G. and Van Hoyweghen, I. (2018). Enacting Actuarial Fairness in Insurance: From Fair Discrimination to Behaviour-based Fairness. *Science as Culture* 27.4, pp. 413–438. DOI: 10.1080/09505431.2017.1398223.

Molnar, C., Casalicchio, G., and Bischl, B. (2020). Interpretable Machine Learning – A Brief History, State-of-the-Art and Challenges. *ECML PKDD 2020 Workshops*. Ed. by Koprinska, I. et al. Communications in Computer and Information Science. Cham: Springer International Publishing, pp. 417–431. DOI: 10.1007/978-3-030-65965-3_28.

Nelder, J. A. and Wedderburn, R. W. M. (1972). Generalized Linear Models. *Journal of the Royal Statistical Society. Series A (General)* 135.3, pp. 370–384. DOI: 10.2307/2344614.

Noll, A., Salzmann, R., and Wuthrich, M. V. (2020). Case Study: French Motor Third-Party Liability Claims. SSRN Scholarly Paper. Rochester, NY. DOI: 10.2139/ssrn.3164764.

Nori, H., Jenkins, S., Koch, P., and Caruana, R. (2019). InterpretML: A Unified Framework for Machine Learning Interpretability. DOI: 10.48550/arXiv.1909.09223.

OCDE (2020). The Impact of Big Data and Artificial Intelligence (AI) in the Insurance Sector - OECD. https://www.oecd.org/finance/Impact-Big-Data-AI-in-the-Insurance-Sector.htm.

Ohlsson, E. and Johansson, B. (2010). Non-Life Insurance Pricing with Generalized Linear Models. 1st ed. 2010, Corr. 3rd printing 2014 edition. Heidelberg; New York: Springer.

Owens, E., Sheehan, B., Mullins, M., Cunneen, M., Ressel, J., and Castignani, G. (2022). Explainable Artificial Intelligence (XAI) in Insurance. *Risks* 10.12, p. 230. DOI: 10.3390/risks10120230.

Pasquale, F. (2015). The Black Box Society: The Secret Algorithms That Control Money and Information. Harvard University Press.
Pesantez-Narvaez, J., Guillen, M., and Alcañiz, M. (2019). Predicting Motor Insurance Claims Using Telematics Data—XGBoost versus Logistic Regression. *Risks* 7.2, p. 70. DOI: 10.3390/risks7020070.

Poufinas, T., Gogas, P., Papadimitriou, T., and Zaganidis, E. (2023). Machine Learning in Forecasting Motor Insurance Claims. *Risks* 11.9, p. 164. DOI: 10.3390/risks11090164.

Power, J., Côté, M.-P., and Duchesne, T. (2024). A Flexible Hierarchical Insurance Claims Model with Gradient Boosting and Copulas. *North American Actuarial Journal* 0.0, pp. 1–29. DOI: 10.1080/10920277.2023.2279782.

R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria. URL: http://www.R-project.org.

Richman, R. (2021a). AI in actuarial science – a review of recent advances – part 1. *Annals of Actuarial Science* 15.2, pp. 207–229. DOI: 10.1017/S1748499520000238.

Richman, R. (2021b). AI in actuarial science – a review of recent advances – part 2. *Annals of Actuarial Science* 15.2, pp. 230–258. DOI: 10.1017/S174849952000024X.

Richman, R. (2024). An AI vision for the actuarial profession. *CAS E-Forum*.

Richman, R. and Wüthrich, M. V. (2023). LocalGLMnet: interpretable deep learning for tabular data. *Scandinavian Actuarial Journal* 2023.1, pp. 71–95. DOI: 10.1080/03461238.2022.2081816.

Rudin, C. (2019). Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. *Nature Machine Intelligence* 1.5, pp. 206–215. DOI: 10.1038/s42256-019-0048-x.

Saleem, R., Yuan, B., Kurugollu, F., Anjum, A., and Liu, L. (2022). Explaining deep neural networks: A survey on the global interpretation methods. *Neurocomputing* 513, pp. 165–180. DOI: 10.1016/j.neucom.2022.09.129.

Shapley, L. S. et al. (1953). A value for n-person games.

Therneau, T. and Atkinson, B. (2022). rpart: Recursive Partitioning and Regression Trees. R package version 4.1.19. URL: https://CRAN.R-project.org/package=rpart.

Thielmann, A. F., Kruse, R.-M., Kneib, T., and Säfken, B. (2024). Neural Additive Models for Location Scale and Shape: A Framework for Interpretable Neural Regression Beyond the Mean. *Proceedings of The 27th International Conference on Artificial Intelligence and Statistics*. PMLR, pp. 1783–1791.

Wick, F., Kerzel, U., and Feindt, M. (2019). Cyclic Boosting – an explainable supervised machine learning algorithm. *2019 18th IEEE International Conference On Machine Learning And Applications (ICMLA)*, pp. 358–363. DOI: 10.1109/ICMLA.2019.00067.

Wood, S. N., Goude, Y., and Fasiolo, M. (2022). Interpretability in Generalized Additive Models. *Interpretability for Industry 4.0: Statistical and Machine Learning Approaches*. Springer, Cham, pp. 85–123. DOI: 10.1007/978-3-031-12402-0_4.

Wood, S. (2023). mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation.

Wood, S. N. (2006). Generalized Additive Models: An Introduction with R. 1 edition. Boca Raton, FL: Chapman and Hall/CRC.

Wright, S. J. (2015). Coordinate descent algorithms. *Mathematical Programming* 151.1, pp. 3–34. DOI: 10.1007/s10107-015-0892-3.

Wüthrich, M. V. and Merz, M. (2023). Statistical Foundations of Actuarial Learning and its Applications. Springer Actuarial. Cham: Springer International Publishing. DOI: 10.1007/978-3-031-12409-9.

Wüthrich, M. V. (2023). Model selection with Gini indices under auto-calibration. *European Actuarial Journal* 13.1, pp. 469–477. DOI: 10.1007/s13385-022-00339-9.

Xin, X. and Huang, F. (2024). Antidiscrimination Insurance Pricing: Regulations, Fairness Criteria, and Models. *North American Actuarial Journal* 28.2, pp. 285–319. DOI: 10.1080/10920277.2023.2190528.

Yitzhaki, S. et al. (2003). Gini’s mean difference: A superior measure of variability for non-normal distributions. *Metron* 61.2, pp. 285–316.
Zakrisson, H. and Lindholm, M. (2025). A tree-based varying coefficient model. *Computational Statistics*, pp. 1–30. DOI: 10.1007/s00180-025-01603-8.
# Appendix A Implementation EBM with InterpretML

InterpretML is the first open-source Python package to include an implementation of EBM (Nori et al., 2019). In this appendix, we specify the parameters of this algorithm and explain how to account for interactions between features. The user can control the learning procedure of the algorithm by placing the cursor on the following elements of Table 7:

* GAM structure by selecting the main effect features and interaction features as well as the loss function to optimize (in yellow),

* Tree structure by controlling the complexity of the decision trees used as building blocks (in green),

* Binning procedure by choosing the number of bins (in red),

* Cyclic boosting procedure (in blue)

<table><caption>Table 7: List of parameters for the Explainable Boosting Regressor algorithm in InterpretML.</caption><thead><tr><th>Parameter</th><th>Default</th><th>Explanation</th></tr></thead><tbody><tr><td>feature_names</td><td>None</td><td>List of feature names</td></tr><tr><td>feature_types</td><td>None</td><td>List of feature types</td></tr><tr><td>exclude</td><td>None</td><td>Features to be excluded</td></tr><tr><td>interactions</td><td>0.9</td><td>Interaction terms to be included</td></tr><tr><td>objective</td><td>RMSE</td><td>Objective to optimize</td></tr><tr><td>minSamples_Leaf</td><td>2</td><td>Minimum number of samples in the leaves</td></tr><tr><td>minHessian</td><td>0.0001</td><td>Minimum hessian required to consider a potential split valid</td></tr><tr><td>maxLeaves</td><td>3</td><td>Maximum number of leaves in each tree</td></tr><tr><td>maxBins</td><td>1024</td><td>Maximum number of bins per feature</td></tr><tr><td>maxInteraction bins</td><td>32</td><td>Maximum number of bins per feature for interaction terms</td></tr><tr><td>maxRounds</td><td>25000</td><td>Number of iterations</td></tr><tr><td>smoothingRounds</td><td>200</td><td>Number of initial highly regularized iterations for main effect feature graphs</td></tr><tr><td>interactionSmoothingRounds</td><td>50</td><td>Number of initial highly regularized iterations for interaction feature graphs</td></tr><tr><td>earlyStoppingRounds</td><td>50</td><td>Number of iterations of no improvement to turn on early stopping</td></tr><tr><td>earlyStoppingTolerance</td><td>0.00001</td><td>Tolerance to prevent early stopping</td></tr><tr><td>validationSize</td><td>0.15</td><td>Size of the validation set for early stopping and creation of outer bags</td></tr><tr><td>outerBags</td><td>14</td><td>Number of outer bags for generating error bounds and adding smoothness to graphs</td></tr><tr><td>innerBags</td><td>0</td><td>Subsamples drawn with replacement at the time of growing individual trees during boosting</td></tr><tr><td>learningRate</td><td>0.01</td><td>Learning rate for boosting</td></tr><tr><td>greedyRatio</td><td>1.5</td><td>Proportion of greedy boosting steps relative to cyclic boosting steps</td></tr><tr><td>cyclicProgress</td><td>1</td><td>Proportion of boosting cycles actively improving model's performance</td></tr></tbody></table>

Interaction terms to be included in the model can be passed on to the `interactions` parameter in three different ways. An integer specifies the count of interactions to be automatically selected, a percentage determines the integer count of interactions by multiplying the number of features by this percentage, or a list of tuples explicitly gives the features to be considered as interaction terms. Monotonic constraints for each feature's relationship with the target may be specified during the fitting process using `monotone_constraints` parameter or ex post using `monotonize` function based on isotonic regression. The InterpretML package also offers a visualization platform for global and local explanations obtained with EBM.

# Appendix B Data description

Figure 15 presents a distribution histogram of each feature included in the claim frequency dataset, along with the empirical claim frequency. Figure 16 displays a distribution histogram of each feature included in the claim severity dataset, along with the evolution of the empirical average claim cost.
<figure><img src="image_16.png" /><figcaption>(a) Driver Age</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Secondary Novice Driver</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(e) Vehicle Price Class</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(g) Zoning Plan</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) License Seniority</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(d) CRM Coefficient</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(f) Vehicle Use</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(h) Year</figcaption></figure>

Figure 15: Univariate analysis of features for claim frequency. For each feature, the bar chart shows the exposure-to-risk in thousands (left y-axis) and the line chart displays the empirical claim frequency (right y-axis).
<figure><img src="image_24.png" /><figcaption>(a) Driver Age</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c) Secondary Novice Driver</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(e) Vehicle Price Class</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(g) Zoning Plan</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(b) License Seniority</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(d) CRM Coefficient</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(f) Vehicle Age</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(h) Year</figcaption></figure>

Figure 16: Univariate analysis of features for claim severity. For each feature, the bar chart shows the number of claims in thousands (left y-axis) and the line chart displays the empirical average claim cost (right y-axis).
# Appendix C Further analysis with claim frequency and severity data

This appendix presents various complementary analyses conducted on the claim frequency and severity data. Section C.1 provides the coefficients of the GLM models fitted on claim frequency and severity data. Section C.2 extends the analysis presented in Section 5.1.1. Finally, Section C.3 supplements the analysis conducted in Section 5.2.2.

## C.1 GLM coefficients

GLM coefficients and multipliers are presented in Table 8 and Table 9 for the claim frequency and severity models respectively. For example, the Vehicle Price Class multiplier of 1.03 in the claim severity model indicates that moving up by one vehicle price class increases the predicted severity by 3%. A similar interpretation can be applied to the other features in the model.

<table><caption>Table 8: GLM model coefficients, multipliers, Z-values and p-values for claim frequency on the training set.</caption><thead><tr><td>Feature</td><td>Coefficient</td><td>Multiplier</td><td>Standard Error</td><td>Z-value</td><td>p-value</td></tr></thead><tbody><tr><td>Driver Age</td><td>0.0176</td><td>1.0178</td><td>0.0007</td><td>25.09</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>License Seniority</td><td>-0.0080</td><td>0.9920</td><td>0.0008</td><td>-10.50</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>Secondary Novice Driver 1</td><td>0.4674</td><td>1.5958</td><td>0.0167</td><td>27.96</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>CRM Coefficient</td><td>1.9514</td><td>7.0385</td><td>0.0348</td><td>56.08</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>Vehicle Price Class</td><td>0.0224</td><td>1.0227</td><td>0.0009</td><td>23.98</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>Vehicle Use Professional</td><td>0.2912</td><td>1.3380</td><td>0.0167</td><td>17.44</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>Zoning Plan</td><td>0.0076</td><td>1.0076</td><td>0.0005</td><td>16.00</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>Year</td><td>-0.0370</td><td>0.9637</td><td>0.0013</td><td>-28.25</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr></tbody></table>

<table><caption>Table 9: GLM model coefficients, multipliers, Z-values and p-values for claim severity on the training set.</caption><thead><tr><td>Feature</td><td>Coefficient</td><td>Multiplier</td><td>Standard Error</td><td>Z-value</td><td>p-value</td></tr></thead><tbody><tr><td>Driver Age</td><td>-0.0027</td><td>0.9973</td><td>0.0007</td><td>-3.61</td><td>0.0003</td></tr><tr><td>License Seniority</td><td>-0.0026</td><td>0.9974</td><td>0.0008</td><td>-3.24</td><td>0.0012</td></tr><tr><td>Secondary Novice Driver 1</td><td>0.2165</td><td>1.2417</td><td>0.0170</td><td>12.77</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>CRM Coefficient</td><td>0.7911</td><td>2.2058</td><td>0.0363</td><td>21.79</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>Vehicle Price Class</td><td>0.0305</td><td>1.0310</td><td>0.0010</td><td>31.06</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>Vehicle Age</td><td>-0.0125</td><td>0.9876</td><td>0.0008</td><td>-15.87</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>Zoning Plan</td><td>0.0088</td><td>1.0088</td><td>0.0005</td><td>17.07</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr><tr><td>Year</td><td>0.0286</td><td>1.0290</td><td>0.0013</td><td>21.51</td><td>&lt; 2 &middot; 10<sup>-16</sup></td></tr></tbody></table>

## C.2 Shape function for the claim frequency model

Figures 17 and 18 respectively present all the main effect univariate and pairwise interaction terms identified in the claim frequency model. Interestingly, the interaction heatmap between the Secondary Novice Driver and Driver Age in Figure 18b shows a relative neutrality of modality 0 with respect to Driver Age but clearly highlights a high-risk area for modality 1 when Driver Age ranges from 40 to 45. This pattern can be interpreted as the age when parents lend their vehicle to their children, a commonly observed effect in the car insurance industry.
<figure><img src="image_32.png" /><figcaption>Figure 17: Prediction of the estimated shape functions $f_j(\mathbf{x}_i)$ in the claim frequency model for features Driver Age, License Seniority, Secondary Novice Driver, CRM Coefficient, Vehicle Price Class, Vehicle Use, Zoning Plan and Year. The shaded areas displays an error bar calculated using a bagging procedure. We use the default `validation_size` and `outer_bags` parameters, as defined by the `InterpretML` algorithm described in Appendix A.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>Figure 18: Heatmaps of the predicted pairwise interaction terms $f_{(j_1, j_2)}$ in the claim frequency model. The pairs of variables include: (CRM Coefficient, Second Novice Driver), (Second Novice Driver, Driver age), (Second Novice Driver, License Seniority), (Second Novice Driver, Vehicle Price Class), (Second Novice Driver, Year) and (License Seniority, Zoning Plan).</figcaption></figure>
### C.3 PDP for the claim severity model

<figure><img src="image_34.png" /><figcaption>Figure 19: Partial Dependence Plot (PDP) in EBM and XGB for the claim severity model. We also show the Individual conditional expectation (ICE) curves for some observations.</figcaption></figure>

## Appendix D Comparison of runtime performance

Runtime for the competing models is displayed in Table 10. To maintain the inherent interpretability by keeping the individual terms additive, EBM pays an additional training cost and is therefore slower than similar methods like XGB. However, the additive structure makes EBM one of the fastest models at prediction time (Nori et al., 2019). Furthermore, on smaller datasets such as the one used for claim severity, EBM is significantly faster than GAM splines. Note that for EBM, EBM² and XGB, runtime for the default model is shown. It is worth noting that the GLM, GAM, and CART methods are implemented in R, while the EBM, EBM², and XGB methods are implemented in Python.
<table><caption>Table 10: Training runtime for the claim frequency and severity models in seconds.</caption><thead><tr><th>Model</th><th>GLM</th><th>GAM</th><th>CART</th><th>EBM</th><th>EBM<sup>2</sup></th><th>XGB</th></tr></thead><tbody><tr><th>Frequency</th><td>34</td><td>5079</td><td>8599</td><td>4221</td><td>6014</td><td>139</td></tr><tr><th>Severity</th><td>1</td><td>169</td><td>10</td><td>13</td><td>18</td><td>2</td></tr></tbody></table>