# Zero-Inflated Tweedie Boosted Trees with CatBoost for Insurance Loss
Analytics

Banghee So Corresponding author: Department of Mathematics, Towson University,
7800 York Rd, Towson, 21204, MD, USA. Email<EMAIL_ADDRESS>Emiliano A.
Valdez Department of Mathematics, University of Connecticut, 341 Mansfield
Road, Storrs, CT, 06269-1009, USA. Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we explore advanced modifications to the Tweedie regression
model in order to address its limitations in modeling aggregate claims for
various types of insurance such as automobile, health, and liability.
Traditional Tweedie models, while effective in capturing the probability and
magnitude of claims, usually fall short in accurately representing the large
incidence of zero claims. Our recommended approach involves a refined modeling
of the zero-claim process, together with the integration of boosting methods
in order to help leverage an iterative process to enhance predictive accuracy.
Despite the inherent slowdown in learning algorithms due to this iteration,
several efficient implementation techniques that also help precise tuning of
parameter like XGBoost, LightGBM, and CatBoost have emerged. Nonetheless, we
chose to utilize CatBoost, a efficient boosting approach that effectively
handles categorical and other special types of data. The core contribution of
our paper is the assembly of separate modeling for zero claims and the
application of tree-based boosting ensemble methods within a CatBoost
framework, assuming that the inflated probability of zero is a function of the
mean parameter. The efficacy of our enhanced Tweedie model is demonstrated
through the application of an insurance telematics dataset, which presents the
additional complexity of compositional feature variables. Our modeling results
reveal a marked improvement in model performance, showcasing its potential to
deliver more accurate predictions suitable for insurance claim analytics.

Keywords: Tweedie, zero-inflated models, tree-based boosting, gradient
boosting, CatBoost, compositional features, aggregate claims models.
JEL classification: C46, C53, G22

##  1 Introduction

Insurance loss analytics, which is the backbone of modern insurance risk
management strategies, involves a systematic collection, processing, and
analysis of observed insurance loss data to extract insights and produce
predictive models to help quantify and manage risk exposure. Actuaries and
risk managers have historically endured the use of data analytics to uncover
patterns, identify emerging risks, and develop predictive models for assessing
claim probabilities, estimating aggregate claim amounts, and deriving
necessary optimal level of reserves to set aside to cover these future
estimated claims. See Klugman et al., (2019)

The two-part frequency-severity models have historically been the norm used in
insurance loss data analytics for modeling aggregate claims for pure premium.
These models recognize that the frequency of claims (how often claims occur)
and the severity of claims (the amount paid per claim) are often influenced by
different factors and follow different distributions. The claim frequency
model focuses on predicting the number of claims that occur within a specified
period, such as a year or a month, and often follows a count distribution
model, such as Poisson or negative binomial distribution. It considers factors
that influence claim occurrence, such as policyholder characteristics (e.g.,
age, location), policy features (e.g., coverage type, deductible), external
factors (e.g., weather conditions, economic trends), and any relevant
historical data. On the other hand, the severity model focuses on predicting
the monetary value of individual claims, i.e., how much each claim will cost
the insurer. It generally follows a continuous distribution, such as gamma or
log-normal, that captures the variability in claim amounts and their right-
skewed nature. It also considers factors that may influence claim amounts. See
Frees et al., 2014a .

Since the introduction of Tweedie et al., (1984) in the 1980s, the Tweedie
mixture distribution has gained popularity in insurance claims modeling as it
simplifies the modeling process by eliminating the need for separate frequency
and severity models. The distribution itself exhibits a point mass at zero
representing the probability of zero claims, followed by continuous density
function for positive claim values. See Frees et al., 2014a . Such integration
is a more accurate risk assessment as it also captures the inherent
interdependence between the frequency and severity components. Primarily due
to the distribution’s exceptional flexibility to accommodate a diverse range
of data, Tweedie regression models have a proven track record of providing
superior fits to historical insurance claims data when compared to the
conventional two-part models.

With the influx of machine learning, the abovementioned conventional methods
started to adopt powerful boosting algorithms. Gradient boosting has become
immensely popular in the field of machine learning and is increasingly being
adopted in insurance and actuarial science for many reasons, among which
include predictive accuracy and robustness. This technique was initially
introduced by Friedman, (2001) and is based around the idea of combining
multiple weak learners, such as small decision trees, to create a single,
strong, and robust predictive model. Suppose we are given a training dataset
$\mathcal{D}=\\{(\bm{x}_{i},y_{i})\\}^{n}_{1}$. Gradient boosting generates a
sequence of functions $W_{0},W_{1},\cdots,W_{T}$, by minimizing the expected
value of a specified loss function, $\ell(y_{i},W_{t})$. This loss function
takes two inputs: the $i$-th response $y_{i}$ and the $t$-th function $W_{t}$,
which estimates $y_{i}$. To improve on these estimates of $y_{i}$, we
iteratively search for another function $W_{t+1}=W_{t}+w_{t+1}(\bm{x})$ that
minimizes the expected value of the loss function, where $w_{t+1}$ is the weak
learner at the $(t+1)$-th iteration. That is,

$w_{t+1}=\operatorname*{arg\,min}_{w\in H}\mathbb{E}[\ell(y,W_{t+1})].$ (1)

Let $H$ be the set of candidate weak learners being evaluated, with the
objective of selecting one to incorporate into the model. Based on the
definition of $W_{t+1}$, we can express the expected value of the loss
function $\ell$ in terms of $W_{t}$ and $w_{t+1}$ as follows:

$\mathbb{E}[\ell(y,W_{t+1})]=\mathbb{E}[\ell(y,W_{t}+w_{t+1})].$ (2)

Our goal is to minimize the value of the loss function with respect to $y$ and
$W_{t}$, while considering the additional factor $w_{t+1}$. Given that the
loss function is continuous and differentiable, this process involves
adjusting $W_{t}$ towards the direction where the loss function decreases most
steeply. This direction corresponds to the rate of change of the loss
function. Thus, by aligning $w_{t+1}$ with the direction where the gradient of
the loss function with respect to $W_{t}$ declines most sharply, we
approximate a $w_{t+1}$ that closely approaches the minimum of
$\sum_{i=1}^{n}\ell(y_{i},W_{t+1}(\bm{x}_{i}))$. The direction of steepest
descent in this context is the negative gradient. Consequently, we define the
pseudo-residuals $r_{i,t+1}$ as follows:

$r_{i,t+1}=-\frac{\partial\ell(y_{i},W_{t}(\bm{x}_{i}))}{\partial
W_{t}(\bm{x}_{i})},$ (3)

and fit $w_{t+1}$ closely to these pseudo-residuals.

Since we are interested in modeling aggregate insurance claims, we design
boosting algorithms to optimize loss functions on the distribution of the
target variable, e.g., Tweedie loss function calculated as the negative log-
likelihood of the data. Insurance claims data frequently features a
significant number of zero values, especially for low-frequency events or
claims. Additionally, the data often includes a small number of high-value
claims, leading to a skewed distribution with extreme values. Given these
inherent characteristics of insurance claims data, choosing a loss function
associated with the Tweedie distribution is a promising strategy when
incorporating gradient boosting methods. In the case of the Tweedie loss
function, it incorporates three parameters $\mu$, $p$, and $\phi$ that all
control the shape of the function. We find from Figure 1(a) the lack of
flexibility of these parameters to control the function in handling the zeros
and right-skewness of the data.

Henceforth, we propose an enhancement to the Tweedie distribution through the
introduction of an inflation probability, denoted as $q$. This modification
leads to the formulation of what we refer to as the Zero-Inflated Tweedie
distribution. In this revised distribution, the probability of zero is
represented by the expression $(q+(1-q)\times f_{Tw}(y=0)$, where
$f_{Tw}(y=0)$ denotes the probability of a zero from the original Tweedie
distribution. In order to provide an even further flexibility in accommodating
the inherent characteristics of the data, we suggest introducing the inflation
probability $q$ as a function dependent on $\mu$, with an additional
parameterization $\gamma$. Figure 1(b) showcases a more flexible $q$ as a
function of $\mu$, addressing the limitations of the Tweedie distribution.

(a) Tweedie distribution

(b) Zero-Inflated Tweedie Distribution

Figure 1: Comparing Tweedie and Zero-Inflated Tweedie Distributions

To enhance the flexibility of modeling insurance losses, we made adjustments
to the conventional gradient boosted trees approach in the following manners.
First, we implemented a customized loss function based on the zero-inflated
Tweedie distribution described above. This decision was made because we
believe the traditional Tweedie distribution is insufficient to provide a more
accurate characterization of insurance losses. While there are some papers
(cite come) that discuss its potential, the use of zero-inflated Tweedie
distribution remains limited in the actuarial and statistical literature.
Second, we next reparameterized the zero-inflated Tweedie loss function to
enable the modeling of the inflation probability $q$ as a function of the mean
parameter $\mu$. Because of this implied relationship between $q$ and $m$,
this results in a singular model that becomes simplified for purposes of model
interpretation and model efficiency. Such approach becomes more attractive to
practitioners and hence, becomes highly likely to be adopted in practice. In
the data analysis section, we demonstrated various ways to interpret the
model’s estimates and outcomes. Third, for empirical demonstration, we
analyzed a synthetic telematics dataset drawn from So et al., (2021). This
data contains records of 100,000 automobile insurance policies and showcases a
notable zero-inflation pattern, with only 2,698 policies experiencing at least
one claim. We specifically picked this dataset to allow us to examine the need
for any specific allowance or special treatments related to compositional
features. It is known that when GLM (Generalized Linear Model) is used, we
need adjustments to the compositional feature variables to avoid bias.
However, we found that such adjustments become unnecessary when boosting
algorithms are employed. Finally, we implemented CatBoost, among other
boosting libraries, to train our Gradient Boosted Decision Tree (GBDT) models.
In So, (2024), So explored various boosting libraries, including XGBoost,
LightGBM, and CatBoost, on the same dataset for estimating various frequency
models and found that CatBoost is the most efficient. For our purposes,
CatBoost further offers enhanced efficiency in handling special types of
feature variables, particularly categorical features that are not uncommon in
insurance data.

By leveraging the zero-inflated Tweedie loss function and employing CatBoost,
we were able to create a more flexible and accurate model for predicting
insurance losses. These adjustments not only improved the model’s performance
but also enhanced its interpretability and ability to handle real-world
insurance data effectively.

The structure of this paper is outlined as follows. Section 2 provides an
overview of zero-inflated Tweedie boosted trees, discussing their relevance in
modeling zero-inflated data commonly encountered in various domains. Section 3
provides an in-depth discussion about CatBoost, an efficient gradient boosting
algorithm known for its robustness and performance in handling categorical
variables. Section 4 presents a comprehensive overview of the methodology
employed in the study, including data preprocessing steps, model training
techniques, and evaluation metrics used to assess model performance. Section 5
delves into the empirical analysis of the dataset, showcasing the application
of the proposed methods (zero-inflated Tweedie boosted trees and CatBoost)
using a synthesized data containing telematics information. We implemented
several ways of interpreting the model estimates and outcomes. Section 6
concludes this paper.

##  2 Zero-Inflated Tweedie Boosted Decision Trees

This section provides an overview for a comprehensive understanding of the
Tweedie model, particularly focusing on its application in modeling the
aggregate claim amount. It serves as a precursor to detailed discussions on
the zero-inflated Tweedie model and the zero-inflated Tweedie boosted tree
model. To begin, let us consider a Poisson random variable $N$, which
represents the number of claims, and a sequence of independent and identically
distributed (i.i.d.) gamma random variables $X_{i}$, for $i=1,2,\ldots,N$. By
assuming independence between $N$ and the $X_{i}$’s, we can define the random
variable $Y$ as follows:

$Y=\begin{cases}0,&\text{if }N=0,\\\\[4.30554pt]
X_{1}+X_{2}+\cdots+X_{N},&\text{if }N\geq 1.\end{cases}$ (4)

Here, $Y$ represents the Poisson sum of independent gamma random variables. In
the context of actuarial science, we interpret $Y$ as the aggregate claim
amount, where $N$ signifies the count of claims, and each $X_{i}$ represents
the amount associated with the $i$-th claim. This specific distribution of $Y$
is known as the compound Poisson-gamma distribution.

Smyth, (1996) highlighted that the compound Poisson-gamma distribution falls
within a distinctive class of exponential dispersion distributions known as
Tweedie distributions (Tweedie et al., (1984)). The Tweedie distributions
belong to the linear exponential family, characterized by a disperson
parameter $\phi$ and a power parameter $p$, with known mean and variance
functions as follows:

  * mean: $\mathbb{E}(Y)=\mu$

  * variance: $\text{Var}(Y)=\phi\mu^{p}$

When constructing the compound Poisson-gamma model for the insurance aggregate
claims, we restrict the power parameter $p$ to the interval $1<p<2$. The
Tweedie models, denoted as $\text{Tw}(\mu,\phi,p)$, are defined by the
following probability density function:

$f_{\text{Tw}}(y|\mu,\phi,p)=a(y,\phi,p)\exp\left(\frac{1}{\phi}\left(y\frac{\mu^{1-p}}{1-p}-\frac{\mu^{2-p}}{2-p}\right)\right),\;\;y\geq
0,$ (5)

where $a(\cdot)$ represents a normalizing function, $\mu>0$ denotes the
expected value of $Y$, and $\phi>0$ represents the dispersion parameter.
Throughout this paper, we aim to delve deeper into the Tweedie model described
above, with specific emphasis on expanding its applicability and
understanding.

The boosted tree model is an approach that assumes the logarithm of the
expected value of the target variable $Y$, given a set of features $\bm{x}$,
can be effectively modeled using the prediction scores generated by gradient
boosting. This relationship can be expressed as follows:

$\ln\mathbb{E}(Y\mid\bm{x})=\ln E+W_{T}(\bm{x}),$ (6)

where $\ln E$ is the offset term and $W_{T}(\bm{x})$ denotes the prediction
score produced by summing the prediction of individual trees in the gradient
boosting model, expressed as:

$W_{T}(\bm{x})=w_{1}(\bm{x})+w_{2}(\bm{x})+\cdots+w_{t}(\bm{x})+\cdots+w_{T}(\bm{x}).$

Here, $w_{t}(\bm{x})$ represents the prediction of the $t$-th tree in the
gradient boosting model. This framework allows for a flexible and powerful
modeling of complex relationships between features and the target variable.

In the context of the Tweedie distribution, the Tweedie boosted tree extends
this concept by assuming that the response variable $Y$ follows a Tweedie
model. Thus, the parameter $\mu$ in the Tweedie model aligns with the expected
value $\mathbb{E}(Y\mid\bm{x})$ in the boosted tree model. This integration of
the Tweedie distribution enhances the model’s ability to handle insurance-
related data, which often exhibits characteristics such as zero-inflated and
right-skewness.

While libraries like CatBoost offer built-in parameters for training Tweedie
boosted trees, challenges can arise in scenarios where observed zeros in
insurance claim data exceed that expected under the Tweedie model assumption.
To overcome this challenge, we adopted a “zero-inflated” approach, known as
the zero-inflated Tweedie (ZITw) model.

The ZITw model combines a point mass at zero with the Tweedie model, improving
the accuracy of estimating $\mu$ especially when dealing with zero-claim
observations. The probability density function in the ZITw model can be
formulated as follows:

$f_{\text{ZITw}}(y|\mu,\phi,p,q)=\begin{cases}q+(1-q)\cdot\exp\left(-\displaystyle\frac{1}{\phi}\displaystyle\frac{\mu^{2-p}}{2-p}\right),&\quad\text{if
}y=0\\\\[12.91663pt] (1-q)\cdot
a(y,\phi,p)\cdot\exp\left(\displaystyle\frac{1}{\phi}\left(y\cdot\frac{\mu^{1-p}}{1-p}-\frac{\mu^{2-p}}{2-p}\right)\right),&\quad\text{if
}y>0.\\\ \end{cases}$ (7)

In this formulation, $q$ represents the inflation probability, indicating the
degree of zero inflation. The expected value of $Y$ under the ZITw model is
determined by both $q$ and $\mu$, given by $(1-q)\mu$. Hence, estimating both
$\mu$ and $q$ accurately becomes crucial, and the gradient boosting framework
provides techniques for achieving this, which we further discuss in Section 4.

##  3 Categorical Boosting: CatBoost

Boosting algorithms have gained immense popularity among actuaries as they are
considered one of the most advanced learning algorithms developed in recent
years. This is because boosting is an iterative fitting process that combines
weak learners – those slightly better than random – into a strong learner,
resulting in more precise and accurate predictions. This process is often
referred to as a stage-wise additive model, where a new weak learner is added
at a time while existing weak learners in the model remain fixed and
unchanged. (need actuarial boosting papers: Michel Denuit, check out Montse
Guillen)

The first boosting algorithm, AdaBoost.M1, was developed by Freund and
Schapire, (1997). Initially designed for classification tasks, it has since
been extended to include regression problems (Hastie et al., (2009)).
Friedman, (2001) introduced gradient boosting, an algorithm tailored for
regression, which constructs a predictive model by leveraging numerical
optimization to minimize the model’s loss through a gradient descent process.
Gradient Boosting Decision Tree (GBDT) is a specific type of gradient boosting
that employs decision trees as weak learners. Notable software libraries for
GBDT implementation include XGBoost, LightGBM, and CatBoost.

CatBoost, developed by Yandex ((Prokhorenkova et al.,, 2018)), is recognized
for its effectiveness in handling heterogeneous datasets, a common scenario in
insurance data. Recent studies (So, (2024)) have demonstrated CatBoost’s
superior performance compared to its counterparts when processing insurance
data. This section will elaborate on CatBoost’s distinctive implementation
strategy.

CatBoost is increasing in popularity for its effective implementation of the
Gradient Boosting Decision Tree (GBDT) algorithm. It excels in handling
heterogeneous data that includes both continuous and categorical features,
while employing decision trees as weak learners. As explained in the preceding
section, the process of the tree construction in CatBoost is sequential, with
each new tree aiming to approximate the negative gradients $r$ of the loss
function $\ell$ at the predictions made by the current function $W$. CatBoost
provides a variety of built-in loss functions, such as Poisson and Tweedie,
that can be utilized with GBDT. Additionally, users have the option to
implement customized loss functions by writing functions that return both the
gradient and the Hessian, which are then used to find a tree that minimizes
the loss function. When incorporating a new tree $w$, CatBoost uses a score
function to assess potential trees. For a prospective tree $w$, the score
function is defined as follows:

$S(w(\bm{x}_{i}),r_{i})=-\sum_{i}(w(\bm{x}_{i})-r_{i})^{2}.$ (8)

The structure of the new tree is determined by the index $j$ of specific
features and a cutoff value $c$. The objective is to determine the optimal
values $j$ and $c$ that maximize the score function. Let $x_{i,j}$ denote the
value of the $j$-th feature on the $i$-th instance. For a tree $w$ with a
depth of one, for example, the function $w(\bm{x}_{i})$ and the score function
are given by:

$w(\bm{x}_{i})=\begin{cases}a_{\text{left}},&\text{if }x_{i,j}\leq
c\\\\[4.30554pt] a_{\text{right}},&\text{if }x_{i,j}>c,\end{cases}$

and

$S(w(\bm{x}_{i}),r_{i})=-\left(\sum_{i:x_{i,j}\leq
c}(a_{\text{left}}-r_{i})^{2}+\sum_{i:x_{i,j}>c}(a_{\text{right}}-r_{i})^{2}\right).$
(9)

Given $i$ and $j$, the optimal values for $a_{\text{left}}$ and
$a_{\text{right}}$ can be determined by employing a regularized loss function
$L$. To avoid overfitting, a regularized loss function can be used to limit
the complexity of the trees: the regularization term is added to the loss
function and penalizes the model for having too many leaves or a tree that is
too deep. At the $t$-th iteration, the regularized loss function is expressed
as:

$L=\sum_{i=1}^{n}\ell(y_{i},W_{t}(\bm{x}_{i}))+\sum_{j=1}^{t}u(w_{j}),$ (10)

where

$W_{t}(\bm{x}_{i})=W_{t-1}(\bm{x}_{i})+w_{t}(\bm{x}_{i}),$ (11)

and $u(w_{j})$ serves as the regularization term, penalizing the complexity of
the tree $w_{j}$.

At each iteration $t$, we train a single decision tree and integrate it into
the existing model, as depicted in equation (12).

$L=\sum_{i=1}^{n}\ell(y_{i},W_{t-1}(\bm{x}_{i})+w_{t}(\bm{x}_{i}))+u(w_{t})+\text{constant}.$
(12)

By expanding the regularized loss function up to the second order using Taylor
series and eliminating constants, the regularized loss function at the $t$-th
iteration can be simplified as:

$L=\sum_{i=1}^{n}\left[g_{i,t}w_{t}(\bm{x}_{i})+\frac{1}{2}h_{i,t}w_{t}(\bm{x}_{i})^{2}\right]+u(w_{t}),$
(13)

where $g_{i,t}=\partial_{W_{t-1}}\ell(y_{i},W_{t-1}(\bm{x}_{i}))$ and
$h_{i,t}=\partial^{2}_{W_{t-1}}\ell(y_{i},W_{t-1}(\bm{x}_{i}))$.

When incorporating L2 regularization with parameter $\lambda$ and
restructuring equation (13) by leaves, we obtain:

$L=\sum_{i\in\text{left}}\left[g_{i,t}a_{t,\text{left}}+\frac{1}{2}h_{i,t}a_{t,\text{left}}^{2}\right]+\frac{1}{2}\lambda
a_{t,\text{left}}^{2}+\sum_{i\in\text{right}}\left[g_{i,t}a_{t,\text{right}}+\frac{1}{2}h_{i,t}a_{t,\text{right}}^{2}\right]+\frac{1}{2}\lambda
a_{t,\text{right}}^{2}.$ (14)

We aim to find the optimal value for each leaf $a_{t,.}$ by minimizing the
regularized loss function specified in equation (14). Therefore, by taking
derivatives with respect to each leaf value $a_{t,.}$, we can identify the
optimal values:

$a_{t,\text{left}}=-\frac{\sum_{i\in\text{left}}g_{i,t}}{\sum_{i\in\text{left}}h_{i,t}+\lambda}\qquad\text{and}\qquad
a_{t,\text{right}}=-\frac{\sum_{i\in\text{right}}g_{i,t}}{\sum_{i\in\text{right}}h_{i,t}+\lambda}\;\;.$
(15)

By combining equations (15) and (9), we can determine the optimal values for
$j$ and $c$ that maximize the score function.

Catboost stands out among GBDT libraries due to its unique handling of
categorical variable. It is important to note that categorical features do not
inherently possess a natural ordering, and therefore, before making a split in
the decision tree, categorical features must be converted into numerical
values. CatBoost addresses this challenge by employing a technique known as
“Ordered Target Statistic.” This method involves randomly rearranging the
dataset and then performing target encoding for each instance, taking into
account only the instances that precede the current one. The process of
transforming categorical features into their numerical counterparts in
CatBoost involves the following steps:

  * •

The training instances are permutated in random order.

  * •

The target (or response) variable is converted from a floating point to an
integer, with the target values divided into buckets for regression problems.

  * •

Categorical features are encoded to numerical features.

After splitting target values into $K$ buckets, the target variable takes
integer values in the range from $0$ to $K-1$. For an instance having target
value $i$, the encoding of the categorical feature is expressed as:

$\text{encoding}^{i}_{u}=\frac{\text{countInClass}+\text{prior}}{\text{totalCount}+1}$

Here, $u$ denotes one of the classes for the categorical feature. The term
“countInClass” represents the frequency of instances with class $u$ whose
target values are equal to $i$. “totalCount” is the total number of instances
with class $u$. The “prior” parameter, which is predetermined, is added to the
encoding formula to ensure smoother results.

It is essential to ensure that the target encoding is solely based on
instances preceding the current one to prevent data leakage. Consequently,
both “countInClass” and “totalCount” are computed following this rule.

##  4 Methodology

In conventional zero-inflated models, training is usually conducted separately
for the mean $\mu$ and the inflation probability $q$. This approach requires
the zero-inflated Tweedie (ZITw) boosted tree to use twice the number of trees
compared to the Tweedie (Tw) model because each parameter is independently
estimated as shown below:

$\ln\mu=\ln E+\displaystyle{W^{mean}_{T}(\bm{x})},$ (16)
$\text{logit}(q)=\ln\frac{q}{1-q}=W^{prob}_{T}(\bm{x}).$ (17)

From a practical standpoint, this leads to a final model consisting of two
separate boosted trees, making it challenging to interpret feature effects and
analyze feature interactions. To address this issue, we propose a novel
algorithm that models $q$ as a function of $\mu$, which allows for the direct
estimation of $q$ from $\mu$. This approached is supported by So, (2024), who
demonstrated that a zero-inflated Poisson boosted tree, treating $q$ as a
function of $\mu$, outperformed traditional models when applied to auto
telematics data. The subsequent sections describe the operational mechanics of
the ZITw boosted tree in these contexts.

Moreover, we discuss strategies to address the specific challenges posed by
compositional data. Compositional data is characterized by feature values that
collectively sum to a constant value, typically 1. Typical examples of such
feature variables, which are commonly found in auto telematics data, are the
percentages of time driving for each day of the week. The importance of this
discussion is highlighted by our analysis of auto telematics data using the
proposed techniques. Auto telematics data serves as a prime example of
compositional data, making it a relevant setting for our investigation.

###  4.1 Scenario 1: Functionally unrelated

We refer to Scenario 1 as the methodology that pertains to the case where $q$
is not directly funtionally related to $\mu$. In this case, as indicated by
equations (16) and (17), it becomes crucial to train $W^{mean}_{T}(\bm{x})$
and $W^{prob}_{T}(\bm{x})$ separately. This separation is necessary because
$W^{mean}_{T}(\bm{x})$ represents the prediction score for the mean parameter
$\mu$, which is obtained through an exponential transformation, while
$W^{prob}_{T}(\bm{x})$ denotes the prediction score for the inflation
probability $q$, which is achieved via a sigmoid transformation. The loss
function used in this context is defined as follows:

$\ell(W^{mean}_{T}(\bm{x}),W^{prob}_{T}(\bm{x}))=\begin{cases}-\ln\left(q+(1-q)\cdot\exp\left(-\displaystyle\frac{1}{\phi}\displaystyle\frac{\mu^{2-p}}{2-p}\right)\right),&\quad\text{if
}y=0\\\\[12.91663pt] -\ln(1-q)-\ln
a(y,\phi,p)-\left(\displaystyle\frac{1}{\phi}\left(y\cdot\frac{\mu^{1-p}}{1-p}-\frac{\mu^{2-p}}{2-p}\right)\right),&\quad\text{if
}y>0,\\\ \end{cases}$ (18)

where $\mu=Ee^{W^{mean}_{T}}$ and $q=\text{logit}^{-1}(W^{prob}_{T}(\bm{x}))$.
The gradients and Hessians of the loss function with respect to $W^{mean}_{t}$
are given by:

$g^{mean}_{t+1}=\partial_{W^{mean}_{t}}\ell(y,W^{mean}_{t},W^{prob}_{t})=\begin{cases}\displaystyle\frac{\alpha\cdot\beta}{q+\alpha},&\quad\text{if
}y=0\\\\[12.91663pt]
-\displaystyle\frac{1}{\phi}y\mu^{1-p}+\beta,&\quad\text{if }y>0,\\\
\end{cases}$ (19)

$h^{mean}_{t+1}=\partial^{2}_{W^{mean}_{t}}\ell(y,W^{mean}_{t},W^{prob}_{t})=\begin{cases}\displaystyle\frac{\alpha\cdot\beta((2-p-\beta)\cdot(q+\alpha)+\alpha\cdot\beta)}{(q+\alpha)^{2}},&\quad\text{if
}y=0\\\\[12.91663pt]
-\displaystyle\frac{1}{\phi}y(1-p)\mu^{1-p}+(2-p)\cdot\beta,&\quad\text{if
}y>0,\\\ \end{cases}$ (20)

where
$\alpha=(1-q)\exp\left(-\displaystyle\frac{1}{\phi}\displaystyle\frac{\mu^{2-p}}{2-p}\right)$,
and $\beta=\displaystyle\frac{1}{\phi}\mu^{2-p}$.

Similarly, the gradients and the Hessians of the loss function with respect to
$W^{prob}_{t}$ are given by:

$g^{prob}_{t+1}=\partial_{W^{prob}_{t}}\ell(y,W^{mean}_{t+1},W^{prob}_{t})=\begin{cases}\displaystyle\frac{1}{1+e^{\delta+W^{logit}_{t}}}-\frac{1}{1+e^{W^{logit}_{t}}},&\quad\text{if
}y=0\\\\[12.91663pt]
\displaystyle\frac{e^{W^{logit}_{t}}}{1+e^{W^{logit}_{t}}},&\quad\text{if
}y>0,\\\ \end{cases}$ (21)

$h^{prob}_{t+1}=\partial^{2}_{W^{prob}_{t}}\ell(y,W^{mean}_{t+1},W^{prob}_{t})=\begin{cases}\displaystyle\frac{e^{W^{prob}_{t}}}{\left(1+e^{W^{prob}_{t}}\right)^{2}}-\frac{e^{\delta+W^{prob}_{t}}}{\left(1+e^{\delta+W^{prob}_{t}}\right)^{2}},&\quad\text{if
}y=0\\\\[12.91663pt]
\displaystyle\frac{e^{W^{prob}_{t}}}{\left(1+e^{W^{prob}_{t}}\right)^{2}},&\quad\text{if
}y>0,\\\ \end{cases}$ (22)

where $\mu=Ee^{W^{mean}_{t+1}}$ and
$\delta=\displaystyle\frac{1}{\phi}\displaystyle\frac{\mu^{2-p}}{2-p}$.

To minimize the loss function with two separate models, we adopt the training
technique proposed by Meng et al., (2022). This approach employs the principle
of coordinate descent optimization, where during each iteration, one parameter
(e.g., $\mu$ or $q$) is optimized while keeping the other fixed. This process
is repeated alternately for $\mu$ and $q$ until convergence, leading to
improve predictive accuracy for both parameters. More specifically, during the
$t$-th iteration, with a fixed inflation parameter $q$, we train
$w_{t}^{\text{mean}}$ (the $t$-th decision tree) to estimate the mean
parameter $\mu$. Then, with the estimated $\mu$ held constant, we train
$w_{t}^{\text{prob}}$ to estimate the inflation probability $q$. After
completion of an estimation cycle for both $\mu$ and $q$, the dispersion
parameter $\phi$ is estimated by minimizing the mean deviation described in
subsection 5.1.1. The power parameter is set to a predetermined constant,
typically 1.5. For a detailed explanation of this training algorithm, please
refer to Algorithm 1 in Appendix A.

###  4.2 Scenario 2: Functionally related

We introduce Scenario 2 as a novel methodology wherein the parameter $q$ is
functionally linked to $\mu$. This innovative approach, absent from existing
literature, acknowledges the relationship between the inflation parameter $q$
and the mean parameter $\mu$. Given our limited prior understanding of this
relationship, we propose a refinemen of the parameterization, initially
introduced by Lambert, (1992) for zero-inflated Poisson Generalized Linear
Models (GLM).

Our proposed parameterization is depicted by the following equations:

$\ln\mu=\ln E+W_{T}(\bm{x}),$ $\text{logit}(q)=\ln\frac{q}{1-q}=-\gamma(\ln
E+W_{T}(\bm{x})).$ (23)

From (23), we derive:

$q=\frac{1}{1+\mu^{\gamma}}\;.$ (24)

Equation (24) reveals that for $\gamma>0$, the inflation probability
diminishes as $\mu$ increases. Guided by this assumption, we define the loss
function as follows:

$\ell(y,W_{T}(\bm{x}))=\begin{cases}-\ln\left(1+\mu^{\gamma}\cdot\exp\left(-\displaystyle\frac{1}{\phi}\displaystyle\frac{\mu^{2-p}}{2-p}\right)\right)+\ln(1+\mu^{\gamma}),&\quad
y=0\\\\[12.91663pt] -\gamma\ln\mu+\ln(1+\mu^{\gamma})-\ln
a(y,\phi,p)-\displaystyle\frac{1}{\phi}\left(y\cdot\frac{\mu^{1-p}}{1-p}-\frac{\mu^{2-p}}{2-p}\right),&\quad
y>0,\\\ \end{cases}$ (25)

where $\mu=Ee^{W_{T}}$. The gradient and Hessian of the loss function with
respect to $W_{t}$ are:

$g_{t+1}=\partial_{W_{t}}\ell(y,W_{t})=\begin{cases}-\displaystyle\frac{\zeta\cdot\eta}{1+\zeta}+\kappa,&\quad
y=0\\\\[12.91663pt]
-\gamma+\kappa-\displaystyle\frac{1}{\phi}y\mu^{1-p}+\displaystyle\frac{1}{\phi}\mu^{2-p},&\quad
y>0,\\\ \end{cases}$ (26)

$h_{t+1}=\partial^{2}_{W_{t}}\ell(y,W_{t})=\begin{cases}\displaystyle\frac{\left(\displaystyle\frac{2-p}{\phi}\mu^{2-p}-\eta^{2}\right)\cdot\zeta\cdot(1+\zeta)+\zeta^{2}\cdot\eta^{2}}{(1+\zeta)^{2}}+\kappa^{2}\cdot\mu^{-\gamma},&\quad
y=0\\\\[12.91663pt]
\kappa^{2}\cdot\mu^{-\gamma}-\displaystyle\frac{(1-p)y}{\phi}\mu^{1-p}+\displaystyle\frac{(2-p)}{\phi}\mu^{2-p},&\quad
y>0,\\\ \end{cases}$ (27)

where
$\zeta=\mu^{\gamma}\cdot\exp\left(-\displaystyle\frac{1}{\phi}\frac{\mu^{2-p}}{2-p}\right)$,
$\eta=\gamma-\displaystyle\frac{1}{\phi}\mu^{2-p}$, and
$\kappa=\displaystyle\frac{\gamma\mu^{\gamma}}{1+\mu^{\gamma}}$.

After constructing an additional tree, we update the value of $\mu$, determine
the inflation parameter $q$, and subsequently, determine the dispersion
parameter $\phi$ by minimizing the mean deviance with the updated values of
$\mu$ and $q$. Finally, with $\mu$, $q$ and $\phi$ fixed, we estimate $\gamma$
by minimizing the mean deviance. The power parameter remains consistently set
to a predetermined constant, typically at 1.5. For a detailed description of
the training algorithm, please refer to Algorithm 2 in Appendix A.

###  4.3 Adjustment of Compositional Data

Compositional data is characterized by multiple non-negative features that sum
up to a constant, typically 100% or 1. Due to the inherent statistical
dependence among these features, transformations are often necessary to map
the data onto the real Euclidean space. This transformation facilitates the
application of traditional statistical methodologies (Aitchison, (1994)).
Among the notable transformations are the logratio methods, which include the
centered logratio transformation (CLR), the additive logratio transformation
(ALR), and the isometric logratio transformation (ILR).

When dealing with compositional data comprising $J$ features, denoted as
$\\{\bm{x}_{\cdot 1},\bm{x}_{\cdot 2},\ldots,\bm{x}_{\cdot J}\\}$, where the
features sum to 1, we refer to these features as a $J$-part composition. The
CLR involves taking the logarithm of a part, standardized relative to the
geometric mean of the logarithms of all parts. Mathematically, $\text{CLR}(j)$
is defined as:

$\text{CLR}(j)=\ln\left(\frac{\bm{x}_{\cdot j}}{\left(\prod_{i}\bm{x}_{\cdot
i}\right)^{1/J}}\right),\quad j=1,2,\ldots,J.$

On the other hand, the ALR contrasts a specific reference part with all other
parts. By selecting the $d$-th part as the reference, the ALR transformation
is represented as:

$\text{ALR}(j|d)=\ln\left(\frac{\bm{x}_{\cdot j}}{\bm{x}_{\cdot
d}}\right),\quad j\neq d.$

The ILR transformation, as proposed by Egozcue et al., (2003), is expressed
as:

$\text{ILR}(\bm{x})=R\cdot\text{CLR}(\bm{x}),$

where $\bm{x}$ is a $J\times n$ data matrix comprising $J$ features, and $R$
is a $(J-1)\times J$ matrix satisfying the condition:

$RR^{T}=I_{J-1}.$

The matrix $R$ is referred to as the contrast matrix (Pawlowsky-Glahn et al.,
(2015)). Although these logratio transformation methods are equivalent up to
linear transformations, the ALR does not preserve distances between data
points in the original space. While CLR retains such distances, it can lead to
singular covariance matrices. On the contrary, the ILR overcomes these issues
but results in transformed features that are complex combinations of the
original features, making interpretation additionally challenging.

The purpose of logratio transformations is to eliminate the inherent spurious
correlations among compositional features. Following transformation,
dimensionality reduction and decorrelation of these features become feasible.
However, researchers have explored applying exponential family principal
component analysis (EPCA) to compositional data (Gan and Valdez, (2021), Yin
et al., (2020), Avalos et al., (2018)). To apply EPCA to compositional data,
it is initially transformed into the real Euclidean space via a logratio
transformation, such as CLR. Probabilistic principal components analysis
(PPCA) is a subset of EPCA (Tipping and Bishop, (1999)). Yin et al., (2020)
emphasized the superiority of PPCA in compositional data analysis,
highlighting its ability to restore original features, which is not achievable
with EPCA alone. Therefore, we confine our study by employing PPCA subsequent
to CLR transformation for its analysis.

##  5 Empirical Investigation

In this section, we undertake a comprehensive comparative analysis to evaluate
the effectiveness of the models proposed in Section 4 compared to conventional
models like the Tweedie boosted tree and Tweedie generalized linear model
(GLM). The primary objective of this analysis is to determine the most
suitable aggregate claims models for managing auto claim datasets with
excessive zeros. Additionally, we aim to explore whether the performance of
the boosted model is influenced by the additional complexity introduced by the
compositional features in the dataset; this helps us assess the necessity for
specific treatments related to compositional features.

Our empirical analysis is based on a synthetic telematics dataset developed by
So et al., (2021). This dataset comprises of 100,000 policies and demonstrates
a zero-inflation characteristic, with only 2,698 policies experiencing at
least one claim. For this study, a total of eight different models were
trained:

  1. 1.

Zero-inflated Tweedie boosted tree with scenario 1 (ZITwBT1)

  2. 2.

Zero-inflated Tweedie boosted tree with scenario 2 (ZITwBT2)

  3. 3.

Tweedie boosted tree (TwBT)

  4. 4.

Tweedie GLM (TwGLM)

  5. 5.

ALR

  6. 6.

CLR

  7. 7.

ILR

  8. 8.

PPCA after CLR transformation

For the latter four models listed above, we applied the suitable logratio
transformations only to zero-inflated Tweedie boosted tree with scenario 2 to
avoid overwhelming the reader.

The Tweedie boosted tree (TwBT) and Tweedie GLM (TwGLM) models were
implemented using the CatBoost and statsmodels Python packages, respectively.
These models utilized the built-in loss function during training. On the other
hand, ZITwBT1 and ZITwBT2 were developed using CatBoost with a customized loss
function, as outlined in Algorithms 1 and 2 in Appendix A, to apply the
specified training methods.

For the boosted trees in our study, we carefully set hyperparameters to ensure
robust and comprehensive model training. Specifically, we established the
number of trained trees, denoted as $T$, to be 500. This choice was made to
strike a balance between model complexity and computational efficiency. Next,
we tackled the crucial hyperparameters of learning rate $\alpha$ and L2
penalty parameter $\lambda$. The learning rate, which controls the
contribution of each tree to the overall model, was selected from a grid of
values including 0.01, 0.05, and 0.10. This range allowed us to explore
different rates of learning while avoiding overly aggressive or sluggish
updates. Similarly, the L2 penalty parameter, influencing the regularization
of our model, was chosen from a grid spanning values from 0 to 100 with
intervals of 10. This grid search enabled us to strike a balance between model
complexity and regularization strength, preventing overfitting and improving
generalization.

The ZITwBT1 model required a different approach due to its unique structure.
To train this model effectively, we divided the total number of trained
decision trees $T$ equally between two distinct boostings – one for $\mu$ and
another for $q$. This division ensures fairness in comparison with other
models and in alignment with total iterations across all models. Each boosting
in the ZITwBT1 model underwent $T=250$ iterations, maintaining consistency in
training depth and performance evaluation.

Throughout our experiments, we maintained a fixed maximum tree depth of 10 and
a power parameter set consistently at 1.5 across all models. These choices
were made based on domain knowledge and empirical observations to ensure
stable and reliable model performance. Moving forward, we conducted rigorous
model comparisons based on four key metrics, elaborated in the subsequent
subsection. To ensure unbiased evaluations, we held 20% of the total dataset
as the test set, while leveraging a 3-fold cross-validation method on the
remaining training dataset to determine optimal hyperparameters. This approach
allowed us to assess model performance accurately and make informed decisions
regarding model selection and refinement.

###  5.1 Performance Evaluation

Comparing Tweedie models with zero-inflated Tweedie models presents a complex
challenge because of their distinct assumptions and structural differences.
This study examines a comprehensive evaluation and comparison of these models
based on a diverse set of assessment metrics. Among these metrics are
traditional measures like deviance and Mean Absolute Deviation (MAD), which
assess the overall goodness of fit and prediction accuracy. The study also
utilizes the Vuong test, a statistical tool specifically designed for
comparing non-nested models and determining their relative performance.
Furthermore, the analysis includes two variants of the Gini index, a widely-
used metric for evaluating predictive power and discrimination ability in
models. By employing this array of assessment metrics, this paper aims to
provide a rigorous comparison between Tweedie and zero-inflated Tweedie
models, which can provide insights on their respective strengths and
limitations in modeling insurance claim data and other similar applications.

####  5.1.1 Deviance and Mean Absolute Deviation

Deviance quantifies the difference between a statistical model and a saturated
model, which perfectly predicts the data. This is done by comparing the log-
likelihoods of the two models and is a key concept in GLMs. Although zero-
inflated boosted models are not exactly GLMs, Martin and Hall, (2016)
investigated whether deviance could be similarly defined for zero-inflated
models. Their research, which included simulations and real-world examples,
produced positive results which suggest that deviance can be defined for zero-
inflated models.

The unit deviance for the $i$ -th observation in zero-inflated Tweedie models
can be determined by the equation below:

$d(y_{i},(\hat{\mu}_{i},\hat{\phi},p,\hat{q}_{i}))=\begin{cases}-2\ln\left(\hat{q}_{i}+(1-\hat{q}_{i})\cdot\exp\left(-\displaystyle\frac{\hat{\mu}_{i}^{2-p}}{\hat{\phi}(2-p)}\right)\right),&\text{if
}y_{i}=0\\\
2\left(\displaystyle\frac{1}{\hat{\phi}}\left(\displaystyle\frac{y^{2-p}}{1-p}-\frac{y^{2-p}}{2-p}\right)-\ln(1-\hat{q}_{i})-\displaystyle\frac{1}{\hat{\phi}}\left(\displaystyle\frac{y_{i}\cdot\hat{\mu}_{i}^{1-p}}{1-p}-\frac{\hat{\mu}_{i}^{2-p}}{2-p}\right)\right),&\text{if
}y_{i}>0.\end{cases}$ (28)

The mean deviance represents the average of the unit deviances for a specific
model across all observations. A significant decrease in deviance suggests a
better alignment between the model and the data.

The mean absolute deviation (MAD) is a statistical measure used to quantify
the average absolute difference between the observed values and the predicted
values, defined as:

$\text{MAD}=\frac{1}{n}\sum_{i=1}^{n}|y_{i}-\hat{y}_{i}|.$

The MAD offers a simple way to assess the error in the predictions. A lower
MAD value suggests higher precision in the predictions because it indicates
that, on average, the predicted values are closer to the actual values.

####  5.1.2 Vuong Test

The Vuong test, proposed by Vuong, (1989), provides a reliable method for
comparing the zero-inflated Tweedie model with other non-nested models, such
as the Tweedie model. Let $f_{K}(y_{i}\bm{x}_{i})$ be the probability
distribution functions for the i-th observation from model $K$. When both
models exhibit similar data fit, their respective likelihood functions tend to
be nearly identical. Conversely, any disparities in the likelihood functions
can indicate which model provides a better fit. The Vuong test operates on
this principle. We define $m_{i}$ as:

$m_{i}=\ln\displaystyle\frac{f_{1}(y_{i}|\bm{x}_{i})}{f_{2}(y_{i}|\bm{x}_{i})}$

We will test the hypothesis that the expected value of $m_{i}$ equals zero,
with the Vuong statistic to gauge this:

$V=\displaystyle\frac{\sqrt{n}\left(\displaystyle\frac{1}{n}\sum_{i=1}^{n}m_{i}\right)}{\sqrt{\displaystyle\frac{1}{n}\sum_{i=1}^{n}(m_{i}-\bar{m})^{2}}}$

Under the assumption of the null hypothesis, the Vuong statistic follows an
asymptotic normal distribution. If the statistic surpasses 1.96, it indicates
that the first model is more likely to be correct; conversely, if it falls
below -1.96, the second model is deemed more accurate, at a 5% significance
level.

####  5.1.3 Gini Index ($\text{Gini}^{a}$ and $\text{Gini}^{b}$)

The Gini index is a well-established metric for assessing model prediction
performance. This paper will specifically examine two popular variants of the
Gini index, $\text{Gini}^{a}$ and $\text{Gini}^{b}$, frequently used in
insurance modeling.

$\text{Gini}^{a}$ It has the following expression:

$\text{Gini}^{a}=\displaystyle\frac{\displaystyle\frac{\sum_{i=1}^{n}y_{i}R(\hat{y}_{i})}{\sum_{i=1}^{n}y_{i}}-\sum_{i=1}^{n}\frac{n-i+1}{n}}{\displaystyle\frac{\sum_{i=1}^{n}y_{i}R(y_{i})}{\sum_{i=1}^{n}y_{i}}-\sum_{i=1}^{n}\frac{n-i+1}{n}}$

In this context, $R(y_{i})$ represents the rank of $y_{i}$. The
$\text{Gini}^{a}$ metric depends on the model’s prediction of the expected
claim amount. A higher $\text{Gini}^{a}$ value indicates an improved
discriminatory power of the model.

$\text{Gini}^{b}$ The $\text{Gini}^{b}$ metric, as proposed by Frees et al.,
2014b , offers a sophisticated approach for comparing a collection of models.

Consider two predictive aggregate claim amounts:

  * •

$\hat{y}^{B}_{i}$ from the reference model, and

  * •

$\hat{y}^{P}_{i}$ from the model under comparison.

Their relative difference can be represented as:

$R_{i}=\frac{\hat{y}^{P}_{i}}{\hat{y}^{B}_{i}}.$

Subsequent to determining $R_{i}$, the instances are sorted in ascending order
based on its values. The ordered Lorenz curve can then be drawn as:

$\left(\frac{\sum_{i=1}^{n}\hat{y}^{B}_{i}\cdot I[R_{i}\leq
s]}{\sum_{i=1}^{n}\hat{y}^{B}_{i}},\frac{\sum_{i=1}^{n}y_{i}\cdot I[R_{i}\leq
s]}{\sum_{i=1}^{n}y_{i}}\right)\quad\text{for }s\in[0,1],$

where $I[\cdot]$ is the indicator function.

The Ginib coefficient is computed by doubling the area between the ordered
Lorenz curve and the line $y=x$. When evaluating multiple models, the one with
the highest Ginib coefficient is deemed the most suitable compared to a
reference model. To determine the best model, each model is considered as the
base model in turn, and the maximum Ginib value for each base model is
recorded. The model with the lowest Ginib among these maximum values is then
selected as the optimal model. This selection method, known as the “min-max”
approach, helps identify and discriminate the most dependable model among the
alternatives (Frees et al., 2014b ).

###  5.2 Synthetic Telematics Data

The synthetic telematics dataset developed by So et al., (2021) used in this
study is publicly accessible at http://www2.math.uconn.edu/~valdez/data.html.
This dataset was meticulously generated using Synthetic Minority Over-sampling
Technique (SMOTE), coupled with a feedforward Neural Network (NN), and trained
based on authentic data from a Canadian insurance provider. It consists of
100,000 data samples, with only 2,698 policies reporting at least one claim.
This notable occurrence of zero-inflation within the dataset highlights a
distinct characteristic of the data. The dataset is structured around 52
distinct features categorized into Traditional data (e.g., exposure, insured
age, and gender), Telematics data (e.g., total miles driven, harsh braking,
and harsh acceleration), and Claims Data (encompassing claim counts and
aggregate claim amounts). Evidently, the aggregate claim amount models were
rigorously developed utilizing exposure and 49 features, with the response
variable being the claim amounts. Among these 49 features, five are
categorical variables, each having varying numbers of categories with 2, 2, 4,
2, and 55 categories, respectively. A noteworthy attribute of this dataset is
the presence of a few categorical attributes. This suggests the need for the
selection of models capable of effectively handling such features to ensure
optimal results and accurate data fitting. To this end, the CatBoost library
was employed to implement boosted tree algorithms, which offers the necessary
capabilities to handle these intricate features. For a detailed overview of
the dataset and its attributes, readers are directed to Table 7 in Appendix B,
where a comprehensive description of the variables in the dataset can be
found. Figure 2 provides an overview of the distribution of the aggregate
claim amounts for the synthetic telematics data. The graph shows a point mass
at zero and a histogram of the positive aggregate claim amounts.

Figure 2: Distribution of Aggregate Claim Amounts for the Synthetic Telematics
Data

####  5.2.1 Results of Model Performance

First, we examine how the zero-inflated Tweedie boosted tree models perform
relative to the conventional Tweedie GLM and Tweedie boosted tree models
ignoring the compositional features in the dataset. Table 1 shows that the
zero-inflated Tweedie boosted tree models (ZITwBTs) have similar MAD values,
which are notably lower than those of the Tweedie Generalized Linear Model
(TwGLM) and Tweedie boosted tree model (TwBT). This suggests that the zero-
inflated models provide more accurate predictions. Furthermore, the deviances
observed in the ZITwBT models are consistently lower compared to the TwGLM and
TwBT models, which suggests the suitability of the zero-inflated approach for
modeling our synthetic telematics data. In terms of $\text{Gini}^{a}$, ZITwBT2
has a significantly higher value than the TwGLM, TwBT models, and the model
from scenario 1 (ZITwBT1), which indicates that it performs better in training
our synthetic telematics data. Therefore, Table 1 suggests the strong promise
of the use of the ZITwBT2 model as it outperforms these other models in
training telematics data.

Table 1: Comparison of MAD, Deviance, and $\text{Gini}^{a}$ across 4 Models

| TwGLM | TwBT | ZITwBT1 | ZITwBT2
---|---|---|---|---
MAD | 257 | 193 | 145 | 146
Deviance | 0.065 | 0.062 | 0.0018 | 0.0017
$\text{Gini}^{a}$ | 0.634 | 0.736 | 0.194 | 0.748

Next, we investigate the impact of the compositional features on model
performance, specifically focusing on ZITwBT2 models. Table 2 presents the
outcomes of these models, comparing their performance with and without
adjustments made to compositional features. The ZITwBT2 model, when
implemented without any compositional feature adjustments, exhibits noteworthy
results, as evidenced by the highest observed $\text{Gini}^{a}$ value
alongside the lowest MAD and deviance scores. This particular configuration
emerges as the most reliable among the 8 models scrutinized, considering these
three critical evaluation metrics. This preference of the ZITwBT2 model
without compositional feature adjustments underscores its possible robustness
in capturing the underlying patterns and nuances present in the dataset. By
leveraging compositional features in their raw form, the model leads to more
accurate predictions and a better understanding of the underlying data
dynamics, which sheds light on how to handle compositional features.

Table 2: Comparison of MAD, Deviance, and $\text{Gini}^{a}$ in ZITwBT2 Models
with and without Compositional Data Adjustment

| ZITwBT2 | ZITwBT2
---|---|---
| alr | clr | ilr | clr+PPCA
MAD | 146 | 149 | 176 | 149 | 149
Deviance | 0.0017 | 0.0020 | 0.0026 | 0.0032 | 0.0025
$\text{Gini}^{a}$ | 0.748 | 0.740 | 0.712 | 0.748 | 0.736

Table 3 presents the Vuong statistics calculated for the four distinct models.
This statistical metric, which serves as a comparative tool, allows us to
assess the relative performance of these models. Specifically, when the Vuong
statistic exceeds 2.33 at a 1% significance level, it signifies that the first
model is superior to the second. Conversely, a value lower than -2.33
indicates that the second model outperforms the first. The corresponding
$p$-values, denoted in parentheses within the table, further support these
comparisons. Upon examining the results, a notable trend emerges when
considering the ZITwBT2 model as the first model and TwGLM and TwBT as the
second models. In this scenario, the $p$-values approach zero, signifying a
high level of statistical significance, while the Vuong statistic is positive.
These combined indicators strongly suggest that the ZITwBT2 model exhibits
superior performance compared to both TwGLM and TwBT. Furthermore, when
juxtaposing the ZITwBT2 model as the first model against ZITwBT1 as the second
model, a different pattern emerges. Here, the $p$-values surpass the 1%
threshold, indicating a lack of statistical significance in differentiating
between the performances of these models. In essence, the evidence suggests
that the performances of ZITwBT2 and ZITwBT1 are comparable, without a clear
distinction in superiority at the specified significance level.

Table 3: Vuong Statistics and Corresponding $p$-values across 4 models

| Second Model
---|---
|  | TwGLM | TwBT | ZITwBT1 | ZITwBT2
First Model | TwGLM |  | - | - | -
TwBT | 0.93 (0.35) |  | - | -
ZITwBT1 | 27.56 (0.00) | 17.03 (0.00) |  | -
ZITwBT2 | 27.65 (0.00) | 17.16 (0.00) | 0.08 (0.94) |

  * •

Numbers represent Vuong statistics, with p-values shown in parentheses.

The results depicted in Table 4 provide further intriguing finding - both
ZITwBT2 models, with and without compositional data adjustment, exhibit
equivalent performance levels. This equivalence suggests that the inclusion or
exclusion of compositional data adjustment does not significantly impact the
predictive power of accuracy of the models. This analysis based on Vuong
statistics indicates that the performance of ZITwBT2 models, both with and
without compositional data adjustment, is equivalent. This equivalence further
underscores the stability, consistency, and robustness of the ZITwBT2 models
as reliable predictive models.

Table 4: Vuong Statistics and Corresponding $p$-values in ZITwBT2 Models with
and without Compositional Data Adjustment

| Second Model
---|---
|  |  | ZITwBT2 | ZITwBT2
|  |  | alr | clr | ilr | clr+PPCA
First Model | ZITwBT2 |  | - | - | - | -
ZITwBT2 | alr | -0.30 (0.77) |  | - | - | -
clr | -0.70 (0.48) | -0.53 (0.60) |  | - | -
ilr | -0.88 (0.38) | -0.73 (0.47) | -0.30 (0.76) |  | -
clr+PPCA | -0.64 (0.52) | -0.46 (0.65) | 0.14 (0.89) | 0.42 (0.68) |

  * •

Numbers represent Vuong statistics, with p-values shown in parentheses.

According to Frees et al., 2014b , the model selected by the “min-max” method
stands out as a robust model selection approach, particularly in the context
of insurance data analysis. This method involves a strategic selection process
aimed at identifying the base model that exhibits the least susceptibility to
competing models. In essence, it focuses on identifying the model with the
smallest maximum $\text{Gini}^{b}$ indices across all competing models,
thereby highlighting its robustness and stability under various conditions.
Table 5 presents a comprehensive overview of the results obtained through the
“min-max” method. Specifically, the ZITwBT2 model emerges with the smallest
maximum $\text{Gini}^{b}$ index recorded at 0.127 among all competing models.
This finding further underscores the ZITwBT2 model’s robustness and
resilience, particularly when compared to alternative models in the analysis.

Table 5: $\text{Gini}^{b}$ across 4 Models

| Competing Model
---|---
|  | TwGLM | TwBT | ZITwBT1 | ZITwBT2
Base Model | TwGLM | - | 0.489 | 0.120 | 0.504
TwBT | -0.043 | - | -0.275 | 0.266
ZITwBT1 | 0.695 | 0.598 | - | 0.704
ZITwBT2 | 0.127 | 0.035 | -0.105 | -

Table 6 presents insightful results regarding the performance of the ZITwBT2
model, particularly in terms of the $\text{Gini}^{b}$ metric. Upon
examination, it becomes evident that the ZITwBT2 model, when implemented
without compositional data adjustment, yielded the lowest maximum
$\text{Gini}^{b}$ value among the tested scenarios, reaching a value of 0.128.
The recorded $\text{Gini}^{b}$ value of 0.128 suggests that, contrary to
expectations, adjusting for compositional features did not lead to a
substantial improvement in the ZITwBT2 model’s fit to the telematics data.
This finding challenges the assumption that incorporating compositional
adjustments inherently enhances model accuracy or predictive capabilities in
this context.

Table 6: $\text{Gini}^{b}$ in ZITwBT2 Models with and without Compositional
Data Adjustment

| Competing Model
---|---
|  |  | ZITwBT2 | ZITwBT2
|  |  | alr | clr | ilr | clr+PPCA
Base Model | ZITwBT2 | - | 0.128 | 0.122 | 0.124 | 0.011
ZITwBT2 | alr | 0.160 | - | 0.145 | 0.119 | 0.033
clr | 0.762 | 0.760 | - | 0.755 | 0.735
ilr | 0.152 | 0.167 | 0.138 | - | 0.059
clr+PPCA | 0.335 | 0.349 | 0.342 | 0.304 | -

To summarize, the results indicate that both the ZITwBT1 and ZITwBT2 models
demonstrate comparable accuracy levels and outperform the TwGLM and TwBT
models across metrics such as Mean Absolute Deviation (MAD), deviance, and
Vuong statistics. Examination of the $\text{Gini}^{a}$ and $\text{Gini}^{b}$
measures indicates that the ZITwBT2 model, without any logratio adjustments,
shows up as the most reliable choice. Specifically, incorporating methods to
adjust compositional data do not significantly enhance model performance. This
suggests that Gradient Boosted Decision Tree (GBDT) models exhibit robustness
against compositional features and thus offer notable advantages over
Generalized Linear Models (GLMs) in terms of convenience and effectiveness.
Unlike GLMs which require additional processing steps to handle compositional
features, GBDT models maintain high performance without the need for such
adjustments.

####  5.2.2 Interpretation

Using CatBoost for constructing an aggregate claim amount model offers
numerous advantages, especially due to its comprehensive set of model analysis
tools that facilitate the interpretation and assessment of diverse risk
features’ effects and interactions on claim amounts. One of the key methods
CatBoost provides is the computation of feature importance values, which plays
a crucial role in understanding the model’s decision-making process. Feature
importance values in CatBoost are calculated based on the impact each feature
has on the change in predictions. This impact is measured by examining how
distinctively the values of leaves split for a particular feature are
distributed. The resulting feature importance values are then normalized,
ensuring that the sum of all feature importances equals 100. This
normalization allows us for a clearer understanding of the relative importance
of each feature in influencing the model’s predictions. In Figure 3, we
present feature importance for the ZITwB2 model, which has been identified as
the most reliable model based on prior analyses conducted on synthetic
telematics data. In the figure, we ranked feature importance for telematics
variables first and later followed by the traditional variables. In
particular, the “Total.miles.driven” feature holds the highest rank in terms
of importance, which indicates its significant impact on predicting aggregate
claim amounts. Following closely are features such as “Accel.06miles,”
“Pct.drive.fri,” “Pct.drive.rush am,” and “Left.turn.indensity11,” all of
which are telematics-related features. This underscores the crucial role that
telematics features play in determining the overall claim amount, for which it
might indicate the importance of leveraging such data in insurance risk
assessment and pricing models.

Figure 3: Feature Importance in ZITwBT2 CatBoost

CatBoost provides insights into the strength and impact of relationships
between pairs of features. This analysis is particularly valuable for
understanding how different features influence each other and contribute to
the model’s predictive power. When CatBoost calculates interactions between
two features, it considers all splits of the feature pair across the trees in
the model. This process involves evaluating the magnitude of leaf value
changes when the splits of the two features are either aligned in the same
direction or differ. Essentially, this approach simulates the effect of fixing
one feature while observing how changes in the other feature affect the
model’s predictions. In practice, understanding feature interactions can
significantly enhance our understanding of complex relationships within the
data. For example, Figure 4 depicts the top 10 feature interactions identified
by the ZITwBT2 model, which was trained on synthetic telematics data. Among
these interactions, one of the most prominent ones is observed between the
features “Credit.score” and “Total.miles.driven,” which indicates a strong
relationship between an individual’s credit score and the total miles they
drive. The analysis further reveals other noteworthy interactions where
“Total.miles.driven” consistently appears paired with features such as
“Car.age,” “Avgdays.week,” and “Annual.pct.driven.” These findings provide
valuable insights into how driving behavior, vehicle age, and creditworthiness
may collectively influence certain outcomes or predictions made by the model.
To better understand these interaction effects and their contributions to
model predictions, techniques like SHAP (SHapley Additive exPlanations) values
can be leveraged. SHAP values help elucidate the impact of individual features
and their interactions on model predictions, and this allows for more
interpretable and actionable insights from the model.

Figure 4: Top 10 Features Interaction Strength in ZITwBT2 CatBoost

CatBoost, known for its efficiency in handling categorical variables and
robust performance in various machine learning tasks, can be enhanced further
by integrating it with the SHAP (SHapley Additive exPlanations) Python
library. SHAP values, with roots in game theory principles, provide for a
better understanding of feature importance by quantifying each feature’s
contribution to the model’s prediction for a specific observation Lundberg et
al., (2020). When interpreting SHAP values, a positive value for a feature
indicates that the feature’s presence is pushing the prediction above the
average prediction value for that particular observation. Conversely, a
negative SHAP value suggests that the feature is pulling the prediction below
the average.

Figure 5: SHAP Values of Top Importance Features in ZITwBT2 CatBoost

In Figure 5, the SHAP values for the top importance features are displayed,
which are also depicted in Figure 3. These SHAP values provide valuable
insights into how different features impact the model’s predictions:

  * •

The feature “Total.miles.driven” exhibits a positive correlation with expected
claim amounts. This finding implies that individuals who cover longer
distances are at a heightened risk of encountering substantial claims,
possibly due to increased exposure on the road.

  * •

Interestingly, increases in the values of “Pct.drive.fri” and “Pct.drive.sat”
are associated with lower expected claim amounts. This observation could be
attributed to safer driving practices adopted by individuals during leisure or
family-related trips on weekends.

  * •

Higher values of “Left.turn.intensity 11” and “Left.turn.intensity 12”
correspond to increased expected claim amounts. This suggests that drivers who
frequently execute left turns are more likely to incur significant claim
amounts. This highlights the elevated risk associated with maneuvers that
involve crossing oncoming traffic.

  * •

For features such as “Brake.09miles”, “Brake.06miles”, and “Accel.06miles,”
both extremely high and low values contribute to a greater expected claim
amount. This may indicate potential interactions with other features that
influence claim predictions.

To gain a deeper understanding of the interaction effects between different
features, we utilized scatter plots to visualize the SHAP (SHapley Additive
exPlanations) values associated with each feature across all observations.
These plots were created with points color-coded by another relevant feature,
as illustrated in Figure 6. The strength of interaction between features was
assessed and reflected in Figure 4, resulting in the generation of six
distinct plots. In these scatter plots, the x-axis and the right y-axis
represent normalized values of the primary feature and the interacting
feature, respectively. Meanwhile, the left y-axis corresponds to the SHAP
values attributed to the primary feature, providing a comprehensive view of
the feature interactions within the dataset.

Figure 6: Six Scatter Plots Describing Feature Interaction Through SHAP Values

The plot in the top-left corner of Figure 6 illustrates a clear positive
correlation between the total annual miles driven and the annualized
percentage of time spent on the road. This trend indicates that as both
variables increase, there is a corresponding escalation in the expected claim
amount. A similar positive correlation is observed between
“Total.miles.driven” and “Avgdays.week,” which provides additional support for
the idea that individuals who drive more frequently tend to have higher
expected claim amounts. Interestingly, the graph also reveals that if
individuals drive nearly every day, regardless of the number of miles driven,
the expected claim amount significantly decreases.

An analysis of credit scores reveals a clear correlation between lower scores
and an elevated risk of substantial insurance claims. Intriguingly, this risk
escalates among individuals who log higher mileage. Conversely, higher-than-
average credit scores not only mitigate the likelihood of major claims but
also indicate safer driving behaviors among those with such scores, which
leads to a decrease in claim frequency. Examining the age of vehicles, those
below the average age present a notable risk of larger claims, particularly
when coupled with extended driving distances. However, as vehicles surpass the
average age threshold, the expected claim amount tends to decrease, with a
more pronounced decline observed among drivers covering extensive distances.
This suggests that longer journeys in older vehicles may potentially
contribute to reducing the overall claim amount.

The bottom-left plots in Figure 6 reveal an intriguing insight: drivers who
frequently cover long distances, brake frequently, and execute numerous left
turns tend to have higher insurance claims. On the contrary, individuals who
drive extensively but brake infrequently and make fewer left turns are less
likely to file large claims. This observation underscores the complex
interplay between driving behavior and insurance claim likelihood,
particularly with the importance of understanding driving patterns for risk
assessment and pricing strategies in the insurance industry.

To better understand the impact of key feature variables on aggregate claim
amounts, we present 7 that shows a heat map of aggregate claim amounts based
on the interaction of “Total.miles.driven” and other feature variables. The
top graph shows that a low credit score combined with more miles driven can
lead to increased claim amounts. In contrast, a high credit score, even with
more miles driven, does not result in increased claim amounts. The other two
graphs in the figure show similar interactions between miles driven
interacting with “Car.age” and “Avgdays.week”. Note that the white cells in
the heat map indicate no observations.

Figure 7: Three Heatmaps Describing Feature Interaction with Aggregate Claim
Amount

##  6 Concluding Remarks

Broadly speaking, for short-term insurance contracts, it is not uncommon to
encounter claim portfolios that exhibit a large mass at zero and a right-
skewed distribution for positive claims. This distribution is characterized by
a prominent peak at zero, indicating a high likelihood of no claims being
made. The adoption of “zero-inflated” models, such as the zero-inflated
Poisson, negative binomial, and Tweedie regression models, has become more
prevalent in practice. While GLM has traditionally been a favored choice in
actuarial studies, its linear structure can pose limitations, especially when
confronted with nonlinearities and intricate interactions inherent in datasets
such as auto telematics data, which captures diverse driving behaviors.
Consequently, there is a growing preference for machine learning approaches to
better model these complex interactions. Among these approaches, gradient
boosting techniques have emerged as powerful tools for constructing predictive
models for insurance claims. See Zhou et al., (2022) and and Meng et al.,
(2022).

In this paper, we integrated the advantages of applying a zero-inflated
Tweedie loss function within a gradient boosting algorithm. To achieve this,
we employed several adjustments. Specifically, we reparameterized the zero-
inflated Tweedie loss function to enable us to express the inflation
probability $q$ as a function of $\mu$. This reparameterization resulted in a
singular model, which enabled us to optimally utilize the CatBoost libraries.
This approach enhances interpretation and facilitates better model comparison
by implementing various performance metrics. CatBoost excels at iteratively
transforming weak learners into strong learners. Among existing boosting
libraries, So, (2024) recommended the use of CatBoost for its efficient
implementation in insurance loss estimation.

Importantly, our research demonstrated that boosted models exhibit robustness
against the influence of compositional features (Aitchison,, 1994), in
contrast to GLM models. This robustness implies that boosted models do not
require additional adjustments for compositional data, underscoring their
superiority over GLMs. Additionally, our study identifies significant risk
factors impacting auto insurance claims and meticulously evaluates their
individual and combined effects on claim amounts. Leveraging the CatBoost
framework in model training has provided dual advantages: efficient handling
of categorical features and a comprehensive toolkit for exploring intricate
interaction effects embedded within the zero-inflated Tweedie Boosted model.
Ultimately, our research illuminates the complex interplay of various features
in determining aggregate claim amounts for short-term insurance portfolios.

### Appendix A. Detailed steps of ZITwBT1 and ZITwBT2 algorithms

Input: Training dataset $D=\\{\bm{x}_{i},y_{i}\\}_{1}^{n}$, total iterations
$T$, learning rate $\alpha$, regularization parameter $\lambda$, power
parameter $p$, dispersion parameter $\phi$

Output: Final model: $\hat{\mu}_{i}=w_{i}\exp(W_{T}^{mean}(\bm{x}_{i}))$,
$\hat{q}_{i}=\text{logit}^{-1}(W^{prob}_{T}(\bm{x}_{i}))$

1 Initialize trees: $w^{mean}_{0}=0$, $w^{prob}_{0}=0$, $\phi=1$ ;

2 for  _$t=1,\ldots,T$_ do

3 Compute first(19) and second(20) derivatives
$g^{mean}_{it}=\partial_{W^{mean}_{t-1}}L(y_{i},W^{mean}_{t-1}(\bm{x}_{i}),W^{prob}_{t-1}(\bm{x}_{i}))$
$h^{mean}_{it}=\partial^{2}_{W^{mean}_{t-1}}L(y_{i},W^{mean}_{t-1}(\bm{x}_{i}),W^{prob}_{t-1}(\bm{x}_{i}))$
for each $i=1,\ldots,n$ ;

4 Fit the tree $w^{mean}_{t}$ minimizing regularized loss function(10);

5 Update $\ W^{mean}_{t}=W^{mean}_{t-1}+\alpha w^{mean}_{t}$;

6

7 Compute first(21) and second(22) derivatives
$g^{prob}_{it}=\partial_{W^{prob}_{t-1}}L(y_{i},W^{mean}_{t}(\bm{x}_{i}),W^{prob}_{t-1}(\bm{x}_{i}))$
$h^{prob}_{it}=\partial^{2}_{W^{prob}_{t-1}}L(y_{i},W^{mean}_{t}(\bm{x}_{i}),W^{prob}_{t-1}(\bm{x}_{i}))$
for each $i=1,\ldots,n$ ;

8 Fit the tree $w^{prob}_{t}$ minimizing regularized loss function(10);

9 Update $\ W^{prob}_{t}=W^{prob}_{t-1}+\alpha w^{prob}_{t}$;

10 Update $\phi$ minimizing mean deviation(28) ;

11

12 end for

13Return final prediction scores: $\
W^{mean}_{T}(\bm{x}_{i})=\sum_{t=1}^{T}\alpha w^{mean}_{t}(\bm{x}_{i})$, $\
W^{prob}_{T}(\bm{x}_{i})=\sum_{t=1}^{T}\alpha w^{prob}_{t}(\bm{x}_{i})$;

Algorithm 1 Zero-Inflated Tweedie Boosted Tree Scenario 1 (ZITwBT1)

Input: Training dataset $D=\\{\bm{x}_{i},y_{i}\\}_{1}^{n}$, total iterations
$T$, learning rate $\alpha$, regularization parameter $\lambda$, inflation
parameter $\gamma$, power parameter $p$, dispersion parameter $\phi$

Output: Final model: $\hat{\mu_{i}}=w_{i}\exp(W_{T}(\bm{x}_{i}))$,
$\hat{q_{i}}=\displaystyle\frac{1}{1+\hat{\mu_{i}}^{\gamma}}$

1 Initialize tree: $w_{0}=0$, $\phi=1$ ;

2 for  _$t=1,\ldots,T$_ do

3 Compute the first(26) and second(27) derivatives
$g_{it}=\partial_{W_{t-1}}L(y_{i},W_{t-1}(\bm{x}_{i}))$
$h_{it}=\partial^{2}_{W_{t-1}}L(y_{i},W_{t-1}(\bm{x}_{i}))$ for each
$i=1,\ldots,n$ ;

4 Fit the tree $w_{t}$ minimizing regularized loss function(10);

5 Update $\ W_{t}=W_{t-1}+\alpha w_{t}$;

6 Update $\phi$ minimizing mean deviation(28) ;

7 Update $\gamma$ minimizing mean deviation(28) ;

8

9 end for

10Return final prediction score: $\ W_{T}(\bm{x}_{i})=\sum_{t=1}^{T}\alpha
w_{t}(\bm{x}_{i})$ ;

Algorithm 2 Zero-Inflated Tweedie Boosted Tree Scenario 2 (ZITwBT2)

### Appendix B. Descriptive details of datasets

Table 7: Variable Names and Descriptions for the Synthetic Telematics Dataset

Type | Variable | Description
---|---|---
Traditional | Duration | Total exposure in yearly units
| Insured.age | Age of insured driver
| Insured.sex † | Sex of insured driver: Male, Female
| Car.age | Age of vehicle (in years)
| Marital † | Marital status: Single, Married
| Car.use † | Use of vehicle: Private, Commute, Farmer, Commercial
| Credit.score | Credit score of insured driver
| Region † | Type of region where driver lives: Rural, Urban
| Annual.miles.drive | Annual miles expected to be driven declared by driver
| Years.noclaims | Number of years without any claims
| Territory † | Territorial location of vehicle: 55 labels in {11, 12, 13, . . ., 91}
Telematics | Annual.pct.driven | Annualized percentage of time on the road
| Total.miles.driven | Total distance driven in miles
| Pct.drive.xxx | Percent of driving day xxx of the week: mon/tue/…/sun
| Pct.drive.x hrs | Percent vehicle driven within x hrs: 2hrs/3hrs/4hrs
| Pct.drive.xxx | Percent vehicle driven during xxx: wkday/wkend
| Pct.drive.rush xx | Percent of driving during xx rush hours: am/pm
| Avgdays.week | Mean number of days used per week
| Accel.xxmiles | Number of sudden acceleration 6/8/9/…/14 mph/s per 1000miles
| Brake.xxmiles | Number of sudden brakes 6/8/9/…/14 mph/s per 1000miles
| Left.turn.intensityxx | Number of left turn per 1000miles with intensity xx: 08/09/10/11/12
| Right.turn.intensityxx | Number of right turn per 1000miles with intensity xx: 08/09/10/11/12
Response | NB_Claim | Number of claims on the given policy
| AMT_Claim | Amount of claims on the given policy

  * $\dagger$

Indicates categorical variable.

## References

  * Aitchison, (1994) Aitchison, J. (1994).  Principles of compositional data analysis.  Lecture Notes-Monograph Series, pages 73–81.
  * Avalos et al., (2018) Avalos, M., Nock, R., Ong, C. S., Rouar, J., and Sun, K. (2018).  Representation learning of compositional data.  Advances in Neural Information Processing Systems, 31.
  * Egozcue et al., (2003) Egozcue, J. J., Pawlowsky-Glahn, V., Mateu-Figueras, G., and Barcelo-Vidal, C. (2003).  Isometric logratio transformations for compositional data analysis.  Mathematical Geology, 35(3):279–300.
  * (4) Frees, E. W., Derrig, R. A., and Meyers, G. (2014a).  Predictive modeling applications in actuarial science, volume 1.  Cambridge University Press.
  * (5) Frees, E. W., Meyers, G., and Cummings, A. D. (2014b).  Insurance ratemaking and a gini index.  Journal of Risk and Insurance, 81(2):335–366.
  * Freund and Schapire, (1997) Freund, Y. and Schapire, R. E. (1997).  A decision-theoretic generalization of on-line learning and an application to boosting.  Journal of Computer and System Sciences, 55(1):119–139.
  * Friedman, (2001) Friedman, J. H. (2001).  Greedy function approximation: a gradient boosting machine.  Annals of Statistics, pages 1189–1232.
  * Gan and Valdez, (2021) Gan, G. and Valdez, E. A. (2021).  Compositional data regression in insurance with exponential family pca.  arXiv preprint arXiv:2112.14865.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer: New York.
  * Klugman et al., (2019) Klugman, S. A., Panjer, H. H., and Willmot, G. E. (2019).  Loss Models: From data to Decisions.  John Wiley & Sons: Hoboken, New Jersey.
  * Lambert, (1992) Lambert, D. (1992).  Zero-inflated poisson regression, with an application to defects in manufacturing.  Technometrics, 34(1):1–14.
  * Lundberg et al., (2020) Lundberg, S. M., Erion, G., Chen, H., DeGrave, A., Prutkin, J. M., Nair, B., Katz, R., Himmelfarb, J., Bansal, N., and Lee, S.-I. (2020).  From local explanations to global understanding with explainable ai for trees.  Nature Machine Intelligence, 2(1):56–67.
  * Martin and Hall, (2016) Martin, J. and Hall, D. B. (2016).  R 2 measures for zero-inflated regression models for count data with excess zeros.  Journal of Statistical Computation and Simulation, 86(18):3777–3790.
  * Meng et al., (2022) Meng, S., Gao, Y., and Huang, Y. (2022).  Actuarial intelligence in auto insurance: Claim frequency modeling with driving behavior features and improved boosted trees.  Insurance: Mathematics and Economics, 106:115–127.
  * Pawlowsky-Glahn et al., (2015) Pawlowsky-Glahn, V., Egozcue, J. J., and Tolosana-Delgado, R. (2015).  Modeling and analysis of compositional data.  John Wiley & Sons.
  * Prokhorenkova et al., (2018) Prokhorenkova, L., Gusev, G., Vorobev, A., Dorogush, A. V., and Gulin, A. (2018).  Catboost: unbiased boosting with categorical features.  Advances in neural information processing systems, 31.
  * Smyth, (1996) Smyth, G. K. (1996).  Regression analysis of quantity data with exact zeros.  In Proceedings of the Second Australia-Japan Workshop on Stochastic Models in Engineering, Technology and Management, Gold Coast, Australia, pages 17–19. Citeseer.
  * So, (2024) So, B. (2024).  Enhanced gradient boosting for zero-inflated insurance claims and comparative analysis of catboost, xgboost, and lightgbm.  Scandinavian Actuarial Journal, 0(0):1–23.
  * So et al., (2021) So, B., Boucher, J.-P., and Valdez, E. A. (2021).  Synthetic dataset generation of driver telematics.  Risks, 9(4):58.
  * Tipping and Bishop, (1999) Tipping, M. E. and Bishop, C. M. (1999).  Probabilistic principal component analysis.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 61(3):611–622.
  * Tweedie et al., (1984) Tweedie, M. C. et al. (1984).  An index which distinguishes between some important exponential families.  In Statistics: Applications and new directions: Proc. Indian statistical institute golden Jubilee International conference, volume 579, pages 579–604.
  * Vuong, (1989) Vuong, Q. H. (1989).  Likelihood ratio tests for model selection and non-nested hypotheses.  Econometrica: journal of the Econometric Society, pages 307–333.
  * Yin et al., (2020) Yin, S., Dey, D. K., Valdez, E. A., and Li, X. (2020).  Flexible modeling of hurdle Conway-Maxwell-Poisson distributions with application to mining injuries.  arXiv preprint arXiv:2008.05968.
  * Zhou et al., (2022) Zhou, H., Qian, W., and Yang, Y. (2022).  Tweedie gradient boosting for extremely unbalanced zero-inflated data.  Communications in Statistics - Simulation and Computation, 51(9):5507–5529.