2.1. Contract dataset

In the contract dataset, each row represents a unique insurance contract. Contracts typically last for one year, but there are instances where their duration may be shorter or longer. Each vehicle is observed over one or more contracts; therefore, one vehicle can be represented by one or more rows in this dataset. Based on risk factors, a premium must be computed for each contract. When using a cross-sectional data model, contracts from the same vehicle are assumed to be independent of each other. On the other hand, a longitudinal data model assumes dependence between contracts, allowing it to use information from previous contracts (including traditional risk factors, telematics data, past claims, etc.) to compute the premium. The contract dataset includes attributes commonly used in vehicle insurance pricing models. These traditional risk factors, displayed in Table 1, are recorded for 117,268 insurance contracts initiated between December $\text{15}^{\text{th}}$ , 2015 and December $\text{31}^{\text{st}}$ , 2018.

Table 1. Variables of the contract dataset

In cases where multiple drivers are associated with a particular contract, attributes of the principal driver are used. Additionally, the dataset includes the vehicle identification number (VIN), allowing us to identify the insured vehicle accurately, alongside the reported claim count. As our goal is to perform claim count regression on contracts, the claim count variable will serve as the response for our supervised learning algorithms, namely the log-linear and the CANN models. The 117,268 contracts are associated with 49,671 distinct vehicles, resulting in an average of approximately 2.36 contracts per vehicle. The histogram of the number of contracts per vehicle is shown in Figure 1.

Figure 1. Number of contracts per vehicle.

2.2. Telematics dataset

All 117,268 contracts have been logged using an on-board diagnostics device, capturing driving information. These data are stored as trip summaries in the telematics dataset, which comprises 117,566,259 trips. Each row in the dataset represents a specific trip, and every trip is described by four attributes: the departure and arrival date and time, the distance driven, and the maximum speed reached. Additionally, each trip is associated with a VIN, and with the date information, it is thus possible to link each trip with one of the 117,268 contracts. An extract from the telematics dataset is presented in Table 2.

2.3. Training, validation, and testing datasets

In supervised learning analysis, splitting the available data into training, validation, and testing sets is paramount to ensure the reliability and generalization ability of the learned model. The training set, which usually comprises the largest portion of the data, is used to train the model’s parameters and optimize its performance. However, relying solely on the training set for performance assessment can lead to overfitting, particularly when the model has a high capacity. To address this, the validation set is used during the modeling process to assess the model’s performance on unseen data. It plays an important role in tuning hyperparameters, selecting the optimal model architecture, and preventing overfitting. By assessing the model’s performance on the validation set, one can obtain an estimate of its generalization performance and make necessary adjustments to improve its ability to generalize well to new, unseen examples. However, it is important to note that the back-and-forth process of evaluating the model on the validation set and adjusting its hyperparameters can introduce information leakage from the validation set into the training set. This can create an illusion of better performance than the model would exhibit in real-world scenarios. As a result, the testing set is reserved for the final evaluation of the learned model. It serves as an unbiased assessment of how well the model will perform on completely unseen data. This final evaluation provides an estimate of the model’s true performance and helps determine its reliability in real-world scenarios. By keeping the testing set separate from the training and validation sets, we can ensure an unbiased evaluation and avoid any potential data leakage. We partition the data as outlined in Table 3 for our analysis. Approximately 60% of the vehicles are allocated for training, while approximately 20% is assigned to the validation and testing sets.

Table 2. Extract from the telematics dataset. Dates are displayed in the yyyy-mm-dd format. The actual VINs have been hidden for privacy purposes

Table 3. Data partitioning

3.1. Cross-sectional models

In addition to assuming independence between vehicles, cross-sectional models also assume independence between contracts from the same vehicle. Consequently, these models do not use the history of a vehicle to estimate its future risk. The PMF of the number of claims can thus be written as:

(3.2) \begin{align} \mathbb{P}\left(Y_{it} = y_{it}| \boldsymbol{y}_{i, (1:t-1)}, \boldsymbol{x}_{i, (1:t)}\right) = \mathbb{P}\left(Y_{it} = y_{it}| \boldsymbol{x}_{it}\right), \quad y_{it} \in \mathbb{N}.\end{align}

3.1.1. Poisson regression

The Poisson distribution is widely used in supervised learning analysis for claim count data due to its good properties and simplicity. Under the Poisson specification, the PMF of the claim count for the $t^{\text{th}}$ contract of vehicle i, denoted by $Y_{it}$ , given its predictor vector, denoted by $\boldsymbol{x}_{it}$ , is defined by

(3.3) \begin{align} \mathbb{P}(Y_{it} = y_{it}|\boldsymbol{x}_{it}) = \frac{e^{-\mu(\boldsymbol{x}_{it})}\mu(\boldsymbol{x}_{it})^{y_{it}}}{y_{it}!},\quad \text{for} \quad y_{it} \in \mathbb{N}, \end{align}

with $\mathbb{E}\! \left[ Y_{it}|\boldsymbol{X}_{it} = \boldsymbol{x}_{it} \right] = \text{Var}[{Y_{it}|\boldsymbol{X}_{it} = \boldsymbol{x}_{it}}] = \mu(\boldsymbol{x}_{it})$ . The mean parameter $\mu(\boldsymbol{x}_{it})$ denotes the conditional expectation (and conditional variance) of $Y_{it}$ . The regression function $\mu(\cdot)$ captures the relationship between the predictors $\boldsymbol{x}_{it}$ and the mean parameter in the Poisson distribution, indicating how the conditional expected count is influenced by the predictors. Choosing a functional form for $\mu(\cdot)$ defines a hypothesis function space $\mathcal{H}$ that includes all the candidate functions for modeling $\mu(\cdot)$ . Typically, one aims at minimizing the Poisson deviance lossFootnote ¹ (see subsection 4.1.2 and Table 4.1 of Wüthrich and Merz, Reference Wüthrich and Merz2023) over the training set:

(3.4) \begin{align} \widehat{\mu} = \textrm{argmin}_{\mu \in \mathcal{H}} \left\{ -\frac{1}{|\mathcal{T}_r|} \sum_{(i, t)\in\mathcal{T}_r} y_{it} \ln[\mu(\boldsymbol{x}_{it})] - \mu(\boldsymbol{x}_{it}) - \ln(y_{it}!) \right \}. \end{align}

Note that this is equivalent to maximizing the likelihood function. For some specifications of $\mu(\cdot)$ , notably the log-linear specification, the criterion in Equation (3.4) is convex, which enables various convex optimization techniques to be applied.

Log-linear Poisson regression

In the Poisson regression context, one notable specification for the regression function is the log-linear form, where the mean parameter is expressed as the exponential of a linear function of the predictors:

(3.5) \begin{align} \mu^{\text{LL}}(\boldsymbol{x};\; \boldsymbol{\beta}) = \exp\left\{\langle\boldsymbol{x}, \boldsymbol{\beta}\rangle\right\}\!, \end{align}

where $\boldsymbol{\beta}$ denotes a vector of parameters, and $\langle\boldsymbol{x}, \boldsymbol{\beta}\rangle$ stands for the inner product between the predictor vector $\boldsymbol{x}$ and the coefficient vector $\boldsymbol{\beta}$ . Log-linear Poisson regression has favorable properties, notably its interpretability stemming from the quasi-linearity of the link function $\mu(\cdot)$ . Moreover, when maximum likelihood is used for parameter estimation, this regression model falls within the framework of GLMs. GLMs provide valuable properties, such as the asymptotic Gaussian distribution of the parameters $\boldsymbol{\beta}$ , allowing for the estimation of standard errors, hypothesis testing, and construction of confidence intervals.

However, log-linear regression does have a significant drawback – its regression function, being linear, lacks flexibility. To address this limitation, various techniques can be employed. In fact, any supervised learning technique could be used for the specification of $\mu(\cdot)$ , such as GAMs and tree-based models. Alternatively, one could incorporate non-linearity by computing polynomial terms (e.g., quadratic or cubic terms) of the predictors.

CANN Poisson regression

In some cases, the supervised learning problem may require even more flexibility, and neural networks are particularly useful in such scenarios. Neural networks are formidable function approximation machines, well-known for their ability to estimate a wide range of highly non-linear multivariate functions. One of the key advantages of neural networks is their ability to handle raw and unstructured data effectively. Because we deal with detailed telematics data, this capability forms the basis for adopting the CANN approach of Wüthrich and Merz (Reference Wüthrich and Merz2019), which embeds a classical actuarial model into a neural network architecture. A CANN model consists of two distinct components: the classical regression model component and the neural network component. This architecture offers great flexibility, allowing for seamless integration of any classical model whose regression function is compatible with a neural network architecture. Likewise, the neural network component can employ various types of supervised architectures, such as CNNs, recurrent neural networks (RNNs), and other architectures tailored to the specific problem at hand. The classical regression model provides good initial estimations and serves as a guide for the neural network component. It offers a starting point for the network’s optimization process, enabling faster convergence. The neural network component, in turn, refines the initial estimations, capturing additional signals and uncovering patterns that may have been missed by the classical model alone.

In our specific case, we use log-linear count regression as the classical model and a MLP as the neural network component in the CANN model. As a result, we have the following specification for the regression function:

(3.6) \begin{align} \mu^{\text{CANN}}(\boldsymbol{x};\; \boldsymbol{\beta}, \boldsymbol{\theta}) &= \mu^{\text{LL}}(\boldsymbol{x};\; \boldsymbol{\beta}) \times \mu^{\text{MLP}}(\boldsymbol{x};\; \boldsymbol{\theta}), \end{align}

where $\mu^{\text{MLP}}(\cdot)$ is the regression function of a MLP parametrized with $\boldsymbol{\theta}$ . In a nutshell, an MLP consists of interconnected layers, including an input layer, hidden layers, and an output layer. Each layer applies an affine transformation to the inputs it receives, followed by a non-linear activation function $\phi$ . This combination of linear transformations and non-linear activations allows MLPs to model complex non-linear relationships in the data. To delve into the mathematical description of an MLP, we can break down its structure starting from the input layer and progressing toward the output layer:

1. 1. Input layer ( $l = 0$ ): The input layer consists of $n_0$ nodes representing the input variables $\boldsymbol{x} = [x_1, x_2, \dots, x_{n_0}] = \boldsymbol{z}^{(0)}$ .

2. 2. Hidden layers ( $l = 1, 2, \dots, L - 2$ ): Each hidden layer l consists of $n_l$ nodes, connected to the nodes from the previous layer ( $l - 1$ ) and the nodes in the following layer ( $l + 1$ ). The computations in the hidden layers involve an affine transformation of the inputs $\boldsymbol{z}^{(l-1)}$ followed by the application of the non-linear activation function $\phi$ . Let us denote the weight matrix between layers $l - 1$ and l as $\boldsymbol{W}^{(l)}$ with dimensions $(n_l, n_{l-1})$ and the bias vector as $\boldsymbol{b}^{(l)}$ with dimensions $(n_l, 1)$ . The computations in each hidden layer l can then be expressed as:

   (3.7) \begin{align} \boldsymbol{a}^{(l)} = \boldsymbol{W}^{(l)} \boldsymbol{z}^{(l - 1)} + \boldsymbol{b}^{(l)}, \quad \boldsymbol{z}^{(l)} = \phi\left( \boldsymbol{a}^{(l)}\right), \end{align}
   where $\boldsymbol{a}^{(l)}$ represents the preactivation values in the $l^{\text{th}}$ hidden layer, and $\boldsymbol{z}^{(l)}$ represents the post-activation values.

3. 3. Output layer ( $l = L - 1$ ): The output layer consists of $n_{L-1}$ nodes, representing the final output(s) of the MLP. The computations in the output layer consist of an affine transformation followed by an activation function g. We denote the weight matrix between layers $L - 2$ (last hidden layer) and $L-1$ as $\boldsymbol{W}^{(L-1)}$ with dimensions $(n_{L-1}, n_{L-2})$ and the bias vector as $\boldsymbol{b}^{(L-1)}$ with dimensions $(n_{L-1}, 1)$ . The computations within the output layer can be expressed as:

   (3.8) \begin{align} \boldsymbol{a}^{(L - 1)} = \boldsymbol{W}^{(L - 1)} \boldsymbol{z}^{(L - 2)} + \boldsymbol{b}^{(L - 1)}, \quad \boldsymbol{z}^{(L - 1)} = g\left( \boldsymbol{a}^{(L - 1)}\right). \end{align}

Note that the number of output neurons $n_{L-1}$ should match the number of modeled distribution parameters. In the context of Poisson regression, where we are modeling a single parameter $\mu$ , only one output neuron is necessary. The choice of the output activation function $g(\cdot)$ is important and should be aligned with the specific problem being tackled since it determines the range and properties of the output values. In our case, we need to ensure that the parameter $\mu$ , which represents the expected count, is always positive. While the exponential function is a natural choice to enforce positivity, it can sometimes lead to numerical instability, especially for large input values. As a better alternative, we choose to use the softplus function as the activation function for the output layer, defined as $\zeta(x) = \log(1 + \exp(x))$ . The softplus function is well-behaved even for large input values, mitigating the issue of numerical instability that can arise with the exponential function. The parameters of the generic MLP described above, consisting of weight matrices and bias vectors, can be denoted as:

(3.9) \begin{align} \boldsymbol{\theta} = \left\{ \boldsymbol{W}^{(1)}, \boldsymbol{b}^{(1)}, \boldsymbol{W}^{(2)}, \boldsymbol{b}^{(2)}, \ldots, \boldsymbol{W}^{(L-1)}, \boldsymbol{b}^{(L-1)} \right\}.\end{align}

Naturally, these parameters must be estimated. However, the criterion in Equation (3.4) is typically not convex, making it challenging to find the global minimum of the empirical risk. In practice, the goal is to find a “good enough” local minimum that yields satisfactory performance on the task. Gradient descent algorithms, such as stochastic gradient descent and its variants, are commonly employed to update iteratively the parameters $\boldsymbol{\theta}$ in the direction of the steepest descent. The backpropagation algorithm efficiently computes the gradients and propagates them through the network, enabling parameter updates.

With the introduced notation for the MLP, we can now express the specification in (3.6) as

(3.10) \begin{align} \mu^{\text{CANN}}(\boldsymbol{x};\; \boldsymbol{\beta}, \boldsymbol{\theta}) &= \zeta\left\{\langle\boldsymbol{x}, \boldsymbol{\beta}\rangle + \boldsymbol{a}^{(L - 1)}(\boldsymbol{x};\; \boldsymbol{\theta})\right\}. \end{align}

The corresponding computational graph for $L = 5$ (i.e., 3 hidden layers) is shown in Figure 2. Like a standard MLP, the network parameters $\boldsymbol{\beta}$ and $\boldsymbol{\theta}$ can be estimated using gradient descent. A simplified pseudo-algorithm for the training of the Poisson CANN model is provided in Algorithm 3.1.1. The learning rate $\eta$ is a hyperparameter that determines the step size taken every time a gradient descent step is performed. In other words, it controls how quickly or slowly the network parameters are updated during training. A higher learning rate allows for larger steps, which can lead to faster convergence. However, an excessively high learning rate may cause the optimization process to overshoot or oscillate around the minimum, hindering convergence. Conversely, a very low learning rate might result in slow convergence, requiring more iterations to reach an acceptable solution. Finding the right learning rate is important and is typically an empirical process that requires experimentation and tuning. In practice, mini-batch gradient descent is commonly used for training neural networks. It works by dividing the training data into smaller subsets, called mini-batches, and computing the gradients and parameter updates based on these mini-batches. This approach offers computational efficiency and improved generalization compared to regular gradient descent, making it a preferred choice in practice. For a comprehensive understanding of neural networks, we refer to the excellent book (Goodfellow et al., Reference Goodfellow, Bengio and Courville2016).

Figure 2. CANN architecture for the Poisson specification. The MLP’s preactivation output value $a_1^{(4)}$ is added to the log-linear model’s preactivation output value $\langle \boldsymbol{x}, \boldsymbol{\beta} \rangle$ before being transformed with the softplus function $\zeta(\cdot)$ . The resulting $\mu$ value is compared to the ground truth y using Poisson deviance loss. The architecture shown employs a 3-hidden-layer MLP but can be customized with any number of layers.

Algorithm 1: Parameter estimation procedure – Poisson CANN model

3.1.2. Negative binomial regression

One issue with the Poisson distribution is its equidispersion assumption. Indeed, we have that $\mu(\boldsymbol{x}) = \mathbb{E}\! \left[ Y_{it}|\boldsymbol{X} = \boldsymbol{x} \right] = \text{Var}[{Y_{it}|\boldsymbol{X} = \boldsymbol{x}}]$ . In practice, claim count data often exhibit overdispersion, where the observed variance of the claim count is greater than the mean. To address this limitation, alternative distributions allowing for overdispersion can be used. Among them, the negative binomial distribution (see, for instance, Denuit et al., Reference Denuit, Maréchal, Pitrebois and Walhin2007 and Cameron and Trivedi, Reference Cameron and Trivedi2013) stands out as a common choice. Under the negative binomial specification, the PMF of the claim count for the $t^{\text{th}}$ contract of vehicle i ( $Y_{it}$ ), given its predictor vector ( $\boldsymbol{x}_{it}$ ), can be written as

(3.11) \begin{align} \mathbb{P}(Y_{it} = y_{it}|\boldsymbol{x}_{it}) = \frac{\Gamma(y_{it} + \phi)}{y_{it}! \Gamma(\phi)} \left(\frac{\phi}{\phi + \mu(\boldsymbol{x}_{it})}\right)^\phi \left( \frac{\mu(\boldsymbol{x}_{it})}{\mu(\boldsymbol{x}_{it}) + \phi}\right)^{y_{it}},\quad \text{for} \quad y_{it} \in \mathbb{N}, \end{align}

where $\phi > 0$ is a dispersion parameter. This can be seen as a generalization of the Poisson distribution. Indeed, the Poisson distribution is recovered when $\frac{1}{\phi} \rightarrow 0$ . The first two centered moments are given by:

(3.12) \begin{align} \mathbb{E}\! \left[ Y_{it}|\boldsymbol{X}_{it} = \boldsymbol{x}_{it} \right] = \mu(\boldsymbol{x}_{it}) \quad \text{and} \quad \text{Var}[Y_{it}|\boldsymbol{X}_{it} = \boldsymbol{x}_{it}] = \mu(\boldsymbol{x}_{it}) + \frac{\mu(\boldsymbol{x}_{it})^2}{\phi}.\end{align}

As can be seen, the negative binomial specification assumes overdispersion since $\text{Var}[Y_{it}|\boldsymbol{X}_{it} = \boldsymbol{x}_{it}] > \mathbb{E}\! \left[ Y_{it}|\boldsymbol{X}_{it} = \boldsymbol{x}_{it} \right]$ . One can estimate the parameters of the regression function $\mu(\cdot)$ along with the dispersion parameter $\phi$ by minimizing the negative binomial deviance loss over the training set:

(3.13) \begin{align} \{\widehat{\mu}, \widehat{\phi}\} = \textrm{argmin}_{\mu \in \mathcal{H}, \phi > 0} \left\{ -\frac{1}{|\mathcal{T}_r|}\sum_{(i,t) \in \mathcal{T}_r} \ln \left[\frac{\Gamma(y_{it} + \phi)}{y_{it}!\Gamma(\phi)}\right] + \phi \ln\left[\frac{\phi}{\phi + \mu(\boldsymbol{x}_{it})}\right] + y_{it} \ln\left[\frac{\mu(\boldsymbol{x}_{it})}{\mu(\boldsymbol{x}_{it}) + \phi}\right]\right\}. \end{align}

Log-linear negative binomial regression

As in the Poisson case, a common specification for $\mu(\cdot)$ is the log-linear form, defined in Equation (3.5). In this case, the criterion in (3.13) is convex, and convex optimization can be used to estimate $\boldsymbol{\beta}$ and $\phi$ .

CANN negative binomial regression

Similar to the approach used for the Poisson case, a CANN architecture can be used to model the mean parameter in the negative binomial distribution. The specification for the regression function $\mu(\cdot)$ remains identical to the Poisson case, as defined in Equation (3.10). In order to incorporate the extra distribution parameter $\phi$ , an additional output neuron is introduced in the network. This output neuron is connected to a neural network weight $w_\phi \in \mathbb{R}$ through the softplus function, ensuring that $\phi$ remains positive, that is, $\phi = \zeta(w_\phi)$ . It is important to highlight that the distribution parameter $\phi$ is not directly connected to the input variables $\boldsymbol{x}$ . As a result, no heterogeneity is incorporated into this parameter, and a common estimated value $\widehat{\phi}$ is used for all observations. The exact architecture for the negative binomial CANN model is depicted in Figure 3. The network parameters $\boldsymbol{\beta}$ , $\boldsymbol{\theta}$ , and $w_\phi$ can be learned by minimizing the criterion in Equation (3.13) using a procedure analogous to Algorithm 3.1.1.

Figure 3. CANN architecture for the negative binomial specification. The MLP’s preactivation output value $a_1^{(4)}$ is added to the log-linear model’s preactivation output value $\langle \boldsymbol{x}, \boldsymbol{\beta} \rangle$ before being transformed with the softplus function $\zeta(\cdot)$ to obtain the $\mu$ value of the negative binomial distribution. The $\phi$ value is obtained by transforming a real-valued parameter $w_\phi$ through the softplus function. The resulting parameters $\mu$ and $\phi$ are then compared to the ground truth y using negative binomial deviance loss. The architecture shown employs a 3-hidden-layer MLP but can be customized with any number of layers.

3.2. Longitudinal models

Cross-sectional models assume independence between all contracts. However, in our case, the data exhibit clustering due to contracts being grouped by vehicle. While it is reasonable to assume independence between contracts from distinct vehicles, this assumption is less valid for contracts from the same vehicle. In reality, the claim counts of contracts within the same vehicle may be influenced by shared vehicle-specific characteristics, unobserved risk factors, or policy-level effects, resulting in dependence between observations within each vehicle cluster. To appropriately address this dependence, we transition from cross-sectional to longitudinal models, enabling the introduction of within-vehicle dependence. In the case of claim count data, a longitudinal model can efficiently leverage the history of the vehicles to refine the risk estimation for future contracts.

While various models are available to analyze longitudinal data, such as random effects models, fixed effects models, generalized estimating equations, and autoregressive models, among others, empirical evidence in the context of claim count regression supports the effectiveness of random (or mixed) effects models (see Boucher et al., Reference Boucher, Denuit and Guillén2008). In these models, a random effect, which is a random variable, is introduced in the specified distribution. For instance, in the case of count data, the specified distribution could be the Poisson distribution. The random effect is assumed to follow a certain distribution, such as a normal, gamma, or another appropriate distribution. The inclusion of the random effect allows for capturing the unobserved heterogeneity or individual-specific effects that cannot be accounted for by the observed covariates. It introduces additional variability into the model and accounts for the dependence within clusters. In longitudinal analysis, we need, for each vehicle i, to model the random vector of claim counts $\boldsymbol{Y}_{i, (1:T_i)} = (Y_{i1}, \dots, Y_{i,T_i})$ . The joint PMF can be expressed with

(3.14) \begin{align} \mathbb{P}\left(\boldsymbol{Y}_{i, (1:T_i)} = \boldsymbol{y}_{i, (1:T_i)}| \boldsymbol{x}_{i, (1:T_i)}\right) = \int_{-\infty}^\infty \left(\prod_{t=1}^{T_i}\mathbb{P}(Y_{it} = y_{it}|\boldsymbol{x}_{i, (1:T_i)}, \theta_i)\right) f(\theta_i) d\theta_i, \end{align}

where $f(\theta_i)$ is the PDF of the random effect.

3.2.1. MVNB regression

A MVNB regression model is obtained by introducing a gamma-distributed random effect in the mean parameter of the Poisson distribution. Specifically, we assume that the conditional distribution of $Y_{it}$ , given $\Psi_i = \psi_i$ , follows a Poisson distribution with mean $\mu(\boldsymbol{x}_{it})\psi_i$ , where $\Psi_i$ is a gamma-distributed random variable with mean 1 and variance $1/\phi$ . The density of $\Psi_i$ is given by

(3.15) \begin{align}f_{\Psi_i}(\psi_i) = \frac{\phi^\phi}{\Gamma(\phi)} \psi_i^{\phi - 1} e^{-\phi\psi_i}, \quad \psi_i > 0.\end{align}

By using Equation (3.14), one can derive the joint distribution for the vector of claim counts:

(3.16) \begin{align} \mathbb{P}\left(\boldsymbol{Y}_{i, (1:T_i)} = \boldsymbol{y}_{i, (1:T_i)}| \boldsymbol{x}_{i, (1:T_i)}\right) = \prod_{t=1}^{T_i} \left(\frac{\mu(\boldsymbol{x}_{it})^{y_{it}}}{y_{it!}}\right) \frac{\Gamma(y_{i\bullet} + \phi)}{\Gamma(\phi)} \left(\frac{\phi}{\mu_{i\bullet} + \phi}\right)^\phi \left(\frac{1}{\mu_{i\bullet} + \phi}\right)^{y_{it}},\end{align}

where $\mu_{i\bullet} = \sum_{t = 1}^{T_i} \mu_{it}$ and $y_{i\bullet} = \sum_{t = 1}^{T_i} y_{it}$ . This joint distribution is commonly referred to as MVNB or negative multinomial distribution. Note that the Poisson distribution is retrieved when $\frac{1}{\phi} \rightarrow 0$ . Furthermore, given the past claim history denoted as $\boldsymbol{y}_{i, (1:t-1)} = (y_{i1}, \dots, y_{i,t-1})$ as well as current and past covariate vectors denoted as $\boldsymbol{x}_{i, (1:t)} = (\boldsymbol{x}_{i1}, \dots, \boldsymbol{x}_{it})$ , one can show that the number of claims at time (or contract) t follows a negative binomial distribution. The probability of observing $y_{it}$ claims at time t, given the past claim history as well as past and current covariate vectors, is thus expressed with

(3.17) \begin{align} \mathbb{P}(Y_{it} = y_{it}| \boldsymbol{y}_{i, (1:t-1)}, \boldsymbol{x}_{i, 1:t}) = \frac{\Gamma(y_{it} + \alpha_{it})}{y_{it}! \Gamma(\alpha_{it})} \left(\frac{\gamma_{it}}{\gamma_{it} + \mu(\boldsymbol{x}_{it})}\right)^{\alpha_{it}} \left( \frac{\mu(\boldsymbol{x}_{it})}{\mu(\boldsymbol{x}_{it}) + \gamma_{it}}\right)^{y_{it}},\quad t = 1, 2, \dots, T_i, \end{align}

where $\alpha_{it} = \phi + \Sigma_{it}^{(y)}$ and $\gamma_{it} = \phi + \Sigma_{it}^{(\mu)}$ . $\Sigma_{it}^{(y)} = \sum_{t' = 1}^{t-1} y_{it'}$ and $\Sigma_{it}^{(\mu)} = \sum_{t' = 1}^{t-1} \mu(\boldsymbol{x}_{it'})$ represent the number of past claims and the sum of past $\mu$ values for contract (i, t), respectively. In the special case when $t=1$ , there is no past history and we set $\Sigma_{it}^{(y)} = \Sigma_{it}^{(\mu)} = 0$ , which yields $\alpha_{i1} = \gamma_{i1} = \phi$ . The expected claim count, given the past history, is given by:

(3.18) \begin{align} \mathbb{E}\! \left[ Y_{it}| \boldsymbol{y}_{i, (1:t-1)}, \boldsymbol{x}_{i, 1:t} \right] &= \mu(\boldsymbol{x}_{it})\left(\frac{\phi + \Sigma_{it}^{(y)}}{\phi + \Sigma_{it}^{(\mu)}}\right)\end{align}

(3.19) \begin{align} &= \mu(\boldsymbol{x}_{it})\left(\frac{\alpha_{it}}{\gamma_{it}}\right). \end{align}

Fitting an MVNB model, therefore, amounts to fitting a negative binomial model, where the parameters $\alpha_{it}$ and $\gamma_{it}$ depend on the vehicle’s history. The parameter $\phi$ and the regression function $\mu(\cdot)$ can be estimated by minimizing the deviance loss over the training set. This can be achieved through the following optimization problem:

(3.20) \begin{align} \{\widehat{\mu}, \widehat{\phi}\} = \textrm{argmin}_{\mu \in \mathcal{H}, \phi > 0} \left\{ -\frac{1}{|\mathcal{T}_r|}\sum_{(i,t) \in \mathcal{T}_r} \ln \left[\frac{\Gamma(y_{it} + \alpha_{it})}{y_{it}!\Gamma(\alpha_{it})}\right] + \alpha_{it} \ln\left[\frac{\gamma_{it}}{\gamma_{it} + \mu(\boldsymbol{x}_{it})}\right] + y_{it} \ln\left[\frac{\mu(\boldsymbol{x}_{it})}{\mu(\boldsymbol{x}_{it}) + \gamma_{it}}\right]\right\}. \end{align}

Log-linear MVNB regression

If the specification for $\mu(\cdot)$ is the log-linear form, defined in Equation (3.5), the criterion in (3.20) is convex, and convex optimization can be used to estimate $\boldsymbol{\beta}$ and $\phi$ .

CANN MVNB regression

In the MVNB case, the CANN architecture, as defined in Equation (3.10), can also be used as a specification for the regression function $\mu(\cdot)$ . To incorporate the additional distribution parameters $\alpha_{it}$ and $\gamma_{it}$ , two additional output neurons are introduced in the network, as depicted in Figure 4. The distribution parameters $\alpha_{it}$ and $\gamma_{it}$ stem from a common parameter $\phi > 0$ , and for the $t^{\text{th}}$ contract of vehicle i, we have $\alpha_{it} = \phi + \Sigma_{it}^{(y)}$ and $\gamma_{it} = \phi + \Sigma_{it}^{(\mu)}$ . The neuron representing $\phi$ is connected to a network weight $w_\phi \in \mathbb{R}$ through the softplus function, i.e., $\phi = \zeta(w_\phi)$ . An MVNB CANN model can be trained with backpropagation and gradient descent, as outlined in Algorithm 3.2.1. Notice that for a vehicle i, the parameter $\gamma_{it}$ depends on the $\mu$ parameter values for its past contracts. As the training procedure of the CANN model is iterative, the estimated $\mu$ values change at each iteration. Hence, it is crucial to update $\Sigma_{it}^{(\mu)}$ for each contract (i, t) at every iteration. This updating procedure is carried out in step 2 of Algorithm 3.2.1.

Figure 4. CANN architecture for the MVNB specification. The MLP’s preactivation output value $a_1^{(4)}$ is added to the log-linear model’s preactivation output value $\langle \boldsymbol{x}, \boldsymbol{\beta} \rangle$ before being transformed with the softplus function $\zeta(\cdot)$ to obtain the $\mu$ value of the negative binomial distribution of Equation (3.17). The $\phi$ value is obtained by transforming a real-valued parameter $w_\phi$ through the softplus function. To obtain $\alpha$ , the sum of past claims $\Sigma^{(y)}$ is added to the $\phi$ parameter, while for $\gamma$ , the sum of past $\mu$ values $\Sigma^{(\mu)}$ is added to the same $\phi$ parameter. The resulting distribution parameters $\mu$ , $\alpha,$ and $\gamma$ are then compared to the ground truth y using negative binomial deviance loss. The architecture shown employs a 3-hidden-layer MLP but can be customized with any number of layers.

Algorithm 2: Parameter estimation procedure – MVNB CANN model

4.1. Log-linear models

The Poisson, negative binomial, and MVNB log-linear models are benchmarks for the Poisson, negative binomial, and MVNB CANN models. These models incorporate all 11 traditional risk factors from Table 1, including the real distance driven (although not strictly classified as a traditional risk factor). It is important to note that the expo variable is used as a covariate rather than as an offset. For each contract (i, t), these traditional risk factors are denoted by the vector $\boldsymbol{x}_{it}^{(\text{trad})}$ . Notice that among the 11 traditional risk factors, 4 are categorical: gender, marital_status, pmt_plan, and veh_use. For these risk factors, the approach involves initially grouping all rare categories, defined as those representing 5% or less of the total number of observations, and labeling them as “others.” We then encode them numerically using dummy encoding. All the resulting traditional covariates are then centered and scaled. Moreover, commute_distance contains missing values, which we fill in using median imputation.

Unlike neural networks, log-linear models have a hard time learning features directly from raw data. To put them on an equal footing with the CANN models, we manually engineer features from the telematics data used by these models. These 13 telematics features, described in Table 4, were specifically engineered from the telematics dataset as risk factors potentially correlated with the claiming risk. For each contract (i, t), these numerical handcrafted telematics features are denoted by the vector $\boldsymbol{x}_{it}^{(\text{hand})}$ . Note that these handcrafted telematics features are also centered and scaled prior to being input into the log-linear models. The regression function for the $\mu$ parameter can thus be written as

(4.1) \begin{align} \mu\left(\boldsymbol{x}^{(\text{trad, hand})};\;\boldsymbol{\beta}\right) = \exp\left(\langle \boldsymbol{x}^{(\text{trad, hand})}, \boldsymbol{\beta} \rangle\right),\end{align}

where $\boldsymbol{x}^{(\text{trad, hand})}$ is the concatenation of $\boldsymbol{x}^{(\text{trad})}$ and $\boldsymbol{x}^{(\text{hand})}$ . The parameters estimated on the training set are shown in a table that can be found on the project’s GitHub repository (see (see Duval, Reference Duval2023). Notably, when using telematics information, the estimated $\phi$ parameter in the MVNB log-linear model is higher. A higher $\phi$ value brings the correcting factor in Equation (3.19) closer to one, indicating reduced importance on past experience when telematics features are used. This underscores the relevance of the engineered telematics features.

Table 4. Handcrafted telematics features extracted from the telematics dataset.

$^1$ 20h–0h.

$^2$ 0h–6h.

$^3$ 11h–14h.

$^4$ 17h–20h Monday to Friday.

$^5$ 7h–9h Monday to Friday.

4.2. CANN models

For the CANN regression models, we extract low-level descriptor vectors that are specifically designed to accurately describe the driving patterns within a particular contract, at least with the dataset we have. We expect the MLP component within the CANN models to learn meaningful high-level features from these low-level vectors. The hope is that the learned features in the hidden layers will be more relevant than the handcrafted features of Table 4. Each contract (i, t) is described by the following descriptor vectors, which provide a summary of its telematics information:

\begin{align*} \boldsymbol{x}_{it}^{(h)} &= \left(x_{it, 1}^{(h)}, \dots, x_{it, 24}^{(h)}\right) \in \mathbb{R}^{24},\\ \boldsymbol{x}_{it}^{(d)} &= \left(x_{it, 1}^{(d)}, \dots, x_{it, 7}^{(d)}\right) \in \mathbb{R}^{7}, \end{align*}

\begin{align*} \boldsymbol{x}_{it}^{(a)} &= \left(x_{it, 1}^{(a)}, \dots, x_{it, 14}^{(a)}\right) \in \mathbb{R}^{14},\\ \boldsymbol{x}_{it}^{(m)} &= \left(x_{it, 1}^{(m)}, \dots, x_{it, 16}^{(m)}\right) \in \mathbb{R}^{16},\\ \boldsymbol{x}_{it}^{(k)} &= \left(x_{it, 1}^{(k)}, \dots, x_{it, 10}^{(k)}\right) \in \mathbb{R}^{10}. \end{align*}

• • The elements in vector $\boldsymbol{x}_{it}^{(h)}$ represent the fraction of driving during each of the 24 h of the day. Therefore, $x_{it, j}^{(h)}$ is the fraction of driving during the $j^{\text{th}}$ h of the day for contract (i, t).

• • The elements in vector $\boldsymbol{x}_{it}^{(d)}$ represent the fraction of driving during each of the 7 days of the week. Therefore, $x_{it, j}^{(d)}$ is the fraction of driving during the $j^{\text{th}}$ day of the week for contract (i, t). Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, and Sunday are denoted by $j = 1, 2, 3, 4, 5, 6, 7$ , respectively.

• • The elements in vector $\boldsymbol{x}_{it}^{(a)}$ represent the fraction of trips made in different average speed slots. For instance, $x_{it, j}^{(a)}$ denotes the fraction of trips made at an average speed between $10(j-1)$ and 10j kilometers per hour.

• • The elements in vector $\boldsymbol{x}_{it}^{(m)}$ represent the fraction of trips made in different maximum speed slots. For instance, $x_{it, j}^{(m)}$ denotes the fraction of trips made where the maximum speed reached falls between $10(j-1)$ and 10j kilometers per hour.

• • The elements in vector $\boldsymbol{x}_{it}^{(k)}$ represent the fraction of trips made in different distance slots. For instance, $x_{it, j}^{(k)}$ denotes the fraction of trips between $5(j-1)$ and 5j kilometers.

These descriptor vectors capture specific aspects of the driving patterns, such as hourly, weekly, average speed, and maximum speed distribution, providing valuable information for the MLPs. Since MLPs can only accept vectors as input, we concatenate these five vectors into a global telematics vector:

(4.2) \begin{align} \boldsymbol{x}_{it}^{(tele)} = \left(\boldsymbol{x}_{it}^{(h)}, \boldsymbol{x}_{it}^{(d)}, \boldsymbol{x}_{it}^{(a)}, \boldsymbol{x}_{it}^{(m)}, \boldsymbol{x}_{it}^{(k)}\right). \end{align}

We incorporate this telematics vector into the MLP component of the CANN models, together with the traditional risk factors $\boldsymbol{x}^{(trad)}$ , enabling interactions between telematics and traditional inputs. In contrast, the log-linear part of the CANN models only includes the traditional risk factors due to the difficulty of processing low-level information. The regression function for the $\mu$ parameter can thus be written as

(4.3) \begin{align}\mu^{\text{CANN}}\left(\boldsymbol{x}^{(\text{trad, tele})};\; \boldsymbol{\beta}, \boldsymbol{\theta}\right) &= \zeta\left\{\langle\boldsymbol{x}^{\text{trad}}, \boldsymbol{\beta}\rangle + \boldsymbol{a}^{(L - 1)}(\boldsymbol{x}^{\text{trad,tele}}; \;\boldsymbol{\theta})\right\}.\end{align}

where $\boldsymbol{x}^{(\text{trad, tele})}$ is the concatenation of $\boldsymbol{x}^{(\text{trad})}$ and $\boldsymbol{x}^{(\text{tele})}$ .

The CANN models are trained using the torch library in the R programming language, using mini-batch gradient descent with 256 observations per batch. The optimizer we use to perform gradient descent is the Adam optimizer, which is a fairly popular choice for training neural networks. Additionally, we use the reduce-on-plateau learning rate scheduler, which dynamically adjusts the learning rate based on the model’s performance, automatically reducing it when the improvement plateaus, allowing for better optimization and convergence during training. For the MLP component of our CANN models, we opt for 3 hidden layers $(L = 5)$ with 128, 64, and 32 hidden units, respectively $(n_1 = 128, n_2 = 64, n_3 = 32)$ . Please note that the number of hidden layers and units per layer are hyperparameters that typically require tuning. However, due to our computational constraints, we limit our investigation to this specific architecture, which we believe is a reasonable choice given our input vector’s size. We choose the rectified linear unit as the activation function $\phi(\cdot)$ used in the hidden layers as it is a widely accepted standard choice. Additionally, we add batch normalization and dropout layers in between fully connected layers, which contribute to network stabilization and regularization. Note that a batch normalization layer is present directly after the input layer, which means that the inputs are normalized before entering the network. Batch normalization applies a normalization transformation to the input of a layer by subtracting the mini-batch mean and dividing by the mini-batch standard deviation. By maintaining a stable mean and variance throughout the network, it can mitigate the vanishing or exploding gradients problem, enabling more effective and efficient training. The dropout layers, on the other hand, serve as a regularization technique that helps prevent overfitting. During training, dropout randomly sets a fraction of the hidden units of a given hidden layer to zero at each iteration, which forces the network to learn redundant representations and reduces the reliance on specific features. This regularization technique improves the model’s ability to generalize well to unseen data.

4.3. CANN hyperparameter tuning

To maximize the performance of our CANN models, we use grid search for hyperparameter tuning, with the average loss observed on the validation dataset $\mathcal{V}_a$ as our optimization criterion. Additionally, we incorporate a regularization technique known as “early stopping” to determine the best number of epochs. This approach allows us to prevent overfitting and select the optimal number of epochs based on the lowest average loss achieved during training. We focus on three key hyperparameters: l_start, denoting the initial learning rate used in the reduce-on-plateau learning rate scheduler; p, which represents the probability of dropout in the dropout layers; and factor, indicating the factor by which the learning rate is multiplied upon reaching a plateau. A plateau is the point where there is no observed improvement in the validation loss for two consecutive epochs. Overall, we compute the average validation loss for all 45 combinations derived from the following hyperparameter values:

• • l_start: $0.00001, 0.00005, 0.0001, 0.0005, 0.001$

• • p: $0.2, 0.3, 0.4$ .

• • factor: $0.3, 0.4, 0.5$

Remember that the network parameters in the classical components of the CANN models are initialized with the maximum likelihood estimators of the corresponding log-linear model, which is why we use relatively small learning rates. At the initialization stage, the network already produces reasonable predictions, reducing the need for large gradient descent steps. Considering the relatively small size of our networks, which are less prone to overfitting, we use rather small dropout probabilities p, all below or equal to 0.4. For the factor hyperparameter, three rather standard values (0.3, 0.4, and 0.5) are tested. The validation loss for each of the 45 combinations and the three specifications is presented in Table S1. It is worth noting that each model is trained for 30 epochs, and as early stopping is employed, the displayed average validation loss is based on the optimal number of epochs, which can be less than 30.

As can be seen, for all learning rates higher than 0.00001, the minimum average validation loss is achieved after a very small number of epochs, indicating that the network learns too quickly. Although the negative binomial and MVNB models perform best at a learning rate of 0.0001, we believe that with more epochs, we could achieve a lower average loss with a learning rate of 0.00001. This is particularly true since the average losses are quite similar for lr_start = 0.00001 and lr_start = 0.0001. When examining the first nine rows of Table S1, it becomes apparent that the factor hyperparameter has a negligible effect on the validation loss. On the other hand, the p hyperparameter only seems to have an impact on the validation loss for the Poisson model, performing best when p = 0.4. Although the dropout rate does not significantly affect the performance for both the negative binomial and MVNB models, we also choose p = 0.4 for these two models since the best performance is achieved at a high number of epochs (29 and 30 epochs, respectively). This suggests that with more epochs, there is potential for further performance improvement. Therefore, we select lr_start = 0.00001, factor = 0.3, and p = 0.4 as the hyperparameters for all three specifications. We train the models again on the training set, this time for 100 epochs. The performance of the three models on the validation set is displayed in Table 5. As can be seen, all 3 specifications require 35 epochs to minimize the average validation loss.

Table 5. Optimal CANN models’ performance on the validation set

5.1. Performance assessment on the testing set

After carefully tuning the hyperparameters of our CANN models, we have at hand promising claim count models that are now nearing implementation. The next crucial step is to estimate their generalization capabilities accurately. We assess the models’ generalization performance using the testing set $\mathcal{T}_e$ , which has remained untouched until now. Using this independent dataset, we can estimate the models’ predictive performance on unseen data points and determine their suitability for real-world applications.

We perform a comparative analysis between the CANN and the log-linear benchmark models, which allows us to evaluate our CANN models’ relative performance and effectiveness against established approaches. In addition to log-linear models using telematics handcrafted features, we also assess the performance of log-linear models using telematics information in the form of the global telematics vector of Equation (4.2). These models have access to exactly the same information as the CANN models and therefore allow an alternative insightful comparison.

In order to fully grasp the value of telematics data, we also evaluate the performance of all 6 models (Poisson, negative binomial and MVNB log-linear and CANN models) without telematics, using only the 11 traditional risk factors as inputs. In the CANN models, the MLP component therefore only comprises the 11 traditional risk factors. This analysis helps us understand the contribution of telematics information in improving the predictive power of the models.

Table 6. Performance of the baseline model on the testing set

Finally, we also explore a variant of the training procedure for CANN models in which the parameters in the log-linear component $\boldsymbol{\beta}$ are fixed, meaning they are not trainable. While this approach may result in reduced performance, it offers the benefit of preserving interpretability while still capturing signals that may have been overlooked by the log-linear component. For these three models (Poisson, negative binomial, and MVNB), we choose to initialize the parameters $\boldsymbol{\beta}$ with coefficients from the MVNB log-linear model, as these are estimated considering time dependence between contracts.

All 18 models are trained on the learning set, and their performance is evaluated on the testing set. To assess the performance, we employ three different scoring rules, namely the Poisson deviance, the logarithmic score, and the squared error. For each scoring rule, we compute the average value on the testing set. To assess the magnitude of the achieved performance, we begin by calculating the average scoring rule values for a baseline model. This baseline model is defined as a homogeneous Poisson log-linear model, where the estimation of the mean (and variance) parameter $\mu$ is estimated by the average number of claims per contract observed in the learning set:

(5.1) \begin{align} \hat{\mu} = \frac{1}{|\{\mathcal{T}_r, \mathcal{V}_a\}|}\sum_{(i, t)\in \{\mathcal{T}_r, \mathcal{V}_a\}} y_{it}.\end{align}

The average scoring rule values for this baseline model on the testing set are reported in Table 6. We can then evaluate the performance of each of the six models in terms of percentage improvement over the baseline model, as shown in Table 7. As can be seen, our CANN models consistently outperform their respective log-linear benchmark models across all scoring rules, indicating that the MLP component of the CANN architecture effectively captures signals that have been missed by the log-linear part. The improvement in prediction performance offered by CANN models obviously comes with computational costs. For the Poisson specification, the log-linear model using handcrafted telematics features had a training duration of 20 s, while the CANN model required about 24 min to complete 35 epochs. Furthermore, longitudinal models, including the MVNB CANN model, provide a substantial improvement over cross-sectional models, indicating that traditional and telematics data alone cannot fully capture the extent of the dependence between observations. This supports the use of longitudinal specifications, as they can capture dependence induced by unobserved variables. CANN models trained with fixed $\boldsymbol{\beta}$ parameters perform similarly to their trainable counterparts. Because CANN models using fixed $\boldsymbol{\beta}$ values are more interpretable and faster to train, we recommend using them over their counterpart with trainable parameters in the log-linear component. Lastly, it is worth noting that log-linear models using the global telematics vector are worse than those using handcrafted telematics features, indicating that they have a hard time learning from such detailed low-level data.

Table 7. Performance comparison of the CANN models and their corresponding log-linear benchmark model on the testing set. The percentages denote the improvement in the scoring rule compared to the baseline model

Table 8. Comparison of the sum of predicted values and the sum of actual values for CANN models using telematics

As highlighted in the works of Wüthrich (Reference Wüthrich2022) and (2020), one limitation of neural networks is their lack of the balance property. This property, which GLMs satisfy (provided that the canonical link function is used), ensures that the sum of predicted values matches the sum of actual values. This is indeed highly desirable in pricing models, as it ensures that premiums adequately cover the expected loss costs. Table 8 displays the sum of predicted and actual values on both learning (in-sample) and testing (out-of-sample) sets for CANN models using telematics. As can be seen, the models overshoot the premiums by about 7% for in-sample predictions and about 3.5% for out-of-sample predictions. This is an issue for pricing because premiums that are too high may lead to adverse selection. While our focus is primarily on comparing models based on the quality of individual predictions, it is essential for pricing models that are actually used in practice to meet the balance property. Fortunately, several techniques (see Wüthrich, Reference Wüthrich2020, Reference Wüthrich2022 and Denuit et al., Reference Denuit, Charpentier and Trufin2021) are available to address this balance issue in neural networks, including in CANN models.

5.2. Permutation feature importance and PDPs

One substantial drawback of neural networks is their difficulty of interpretation. However, researchers have developed tools to shed light on the inner workings of these black box algorithms. Two particularly useful tools in this context are permutation feature importance and PDPs.

Figure 5. Importance scores of the 20 most important variables obtained for the MVNB CANN model.

Permutation feature importance is a model-agnostic technique that computes an importance score for each input (or variable) in a supervised learning algorithm. It achieves this by randomly permuting the values of a specific input while holding the other inputs constant and observing the resulting effect on the model’s performance. By comparing the original model’s performance with the permuted performance, we can determine the variable’s importance relative to the chosen performance metric. Suppose we have a trained model and a holdout sample for evaluation purposes. We can initially score the model on this sample and measure its performance using a chosen metric, such as the average loss. Let us denote the average loss obtained with the original holdout sample as $\ell_{\text{original}}$ . To assess the importance of input j in the prediction process, we randomly shuffle the values of input j in the holdout sample and rescore the model. This process yields a new average loss, denoted as $\ell_{\text{permuted}}^{(j)}$ , where the superscript j indicates that input j has been permuted. If input j is indeed important for the model’s prediction, the permuted average loss $\ell_{\text{permuted}}^{(j)}$ is expected to be greater than the original average loss $\ell_{\text{original}}$ . This suggests that permuting the values of input j has a detrimental effect on the model’s performance. To obtain an importance score for input j, we can compute the difference between the permuted average loss and the original average loss, resulting in the feature importance score $\text{FI}_j$ :

\begin{align*} \text{FI}_j = \ell_{permuted}^{(j)} - \ell_{original}.\end{align*}

To obtain a more reliable estimate of the importance score, this procedure can be repeated a certain number of times for input j, creating a distribution of the increase or decrease in the average loss. The whole procedure can then be repeated for all inputs. In Figure 5, the importance scores of the 20 most important variables for our best model, the MVNB CANN, are visualized using boxplots. Please note that the names used for the telematics inputs in Figure 5 differ from the introduced notation. However, a translation table is provided in Table A.1 of Appendix A to clarify the correspondence between the names used and the introduced notation. Each boxplot represents the distribution of the 100 importance scores assigned to a specific input obtained by shuffling and assessing the model 100 times. The performance metric used is the deviance loss (up to a factor 2). The analysis reveals interesting findings regarding the claim count model. As can be seen, the top 5 most important variables are from our set of 11 traditional risk factors. Notably, veh_age, distance, and expo play a significant role in the model’s performance. When it comes to telematics inputs, those related to maximum speed demonstrate a substantial impact on the model’s performance. Particularly, vma_16, representing the fraction of trips made at a maximum speed exceeding 150 kilometers per hour, stands out as the most important input. In general, the fraction of trips made at high maximum speeds, such as vma_14, vma_15, and vma_16, proves to be valuable for predicting claims. Additionally, it is interesting to observe that h_22 and h_2, which represent the fraction of driving during night hours, contribute substantially to the assessment of risk. Importantly, the gender variable, often used by insurers as a risk factor, is rendered useless in the presence of telematics inputs. It ranks as the $\text{70}^{\text{th}}$ most important variable (not showed in Figure 5), indicating its insignificance in the model’s predictive power.

PDPs are valuable tools for understanding the relationship between a specific input variable and the output of a supervised learning model. PDPs are also model-agnostic, meaning they can be applied to different types of models. They provide insights into how changes in a particular input variable influence the model’s predictions while keeping all other variables at fixed values. In other words, they illustrate the marginal effect of an input variable on the predicted outcome. To compute a PDP for a specific input variable j, the process involves the following steps. First, a grid of values is defined to cover the entire or plausible range of the variable’s values. Next, while holding all other variables fixed, the input vector in the holdout dataset is sequentially replaced with each value from the defined grid. Subsequently, predictions are obtained using the trained model on the modified holdout dataset for each value. By plotting the input variable values on the x-axis and the corresponding average prediction on the y-axis, the resulting PDP visually showcases the relationship between the input variable and the model’s predictions. Figure S1 displays the PDPs of the eight most important telematics inputs in the MVNB CANN model. The plots reveal that the risk, expressed as the expected number of claims, appears to increase in a linear fashion with the proportion of trips made at high maximum speeds, indicated by the input variables vma_14, vma_15, and vma_16. Additionally, there appears to be a positive linear association between the expected number of claims and the proportion of driving taking place during nighttime hours, specifically between 9 p.m. and 10 p.m. (h_22) and between 1 a.m. and 2 a.m. (h_2). It is important to emphasize that when interpreting PDPs, caution must be exercised, as the procedure assumes that the input variables are independent of each other. In particular, the interpretation of the PDPs related to the fraction of driving on Tuesdays (p_2) is challenging due to the correlation between the proportions of driving on different days of the week. For instance, if an insured individual drives in smaller proportions on Tuesdays, they will systematically drive in larger proportions on other days of the week.