Mixture composite regression models with multi-type feature selection

Mixture composite regression models with multi-type feature selection

(M.V. W\"{u}thrich).} .} \and George Tzougas.}\\ \and Mario V. W\"{u}thrich\footnotemark[1] }

\abstract{The aim of this paper is to present a mixture composite regression model for claim severity modeling. Claim severity modeling poses several challenges such as multimodality, tail-heaviness and systematic effects in data. We tackle this modeling problem by studying a mixture composite regression model for simultaneous modeling of attritional and large claims, and for considering systematic effects in both the mixture components as well as the mixing probabilities. For model fitting, we present a group-fused regularization approach that allows us for selecting the explanatory variables which significantly impact the mixing probabilities and the different mixture components, respectively. We develop an asymptotic theory for this regularized estimation approach, and fitting is performed using a novel Generalized Expectation-Maximization algorithm. We exemplify our approach on a real motor insurance data set.}

Keywords: Splicing; Generalized Expectation-Maximization algorithm; Variable selection; Asymptotic normal theory; Multimodal and heavy-tailed claim losses

Introduction

Insurance claim severity modeling is a very challenging problem in actuarial science. Motivated by a Greek Motor Third Party Liability (MTPL) insurance data set which is further described in Section (ref), we observe that insurance claim severity data sets often exhibit several peculiar characteristics: Firstly, claim severity distributions are often multimodal, coming from the fact that there are systematic effects in the data due to unobserved heterogeneity and latent factors such as different claim types. Secondly, a claim severity distribution ranges over several magnitudes, from small attritional claims to large claim events, which often exhibits a heavy-tailed nature and a mismatch between body and tail behavior. Thirdly, insurance data are often accompanied by multiple types of policyholder attributes, including several continuous variables (e.g. driver's age), ordered categorical variables (e.g. sum insured categories) and nominal categorical variables (e.g. car brand). These variables may have different explanatory powers to different parts of the severity distribution.

For insurance pricing, reserving and risk management, it is crucial to have accurate descriptions of claim severity distributions, and to understand clearly the influence of policy attributes to the claim distribution. Therefore, it is essential to devise a claim severity modeling framework which possesses all of the following features to address the aforementioned modeling challenges:

There are several actuarial research works in addressing each of the aforementioned modeling requirements. To capture distributional multimodality and tail-heaviness (points 1 and 2), there are three main streams of claim severity modeling approaches described as follows.

To understand how the claim severity distribution is influenced by certain risk factors, covariates influence (point 3) has been extensively explored in actuarial literature using various types of severity regression models (nelder1972generalized). We refer readers to frees2009regression for an extensive summary.

Considering variable selection techniques (point 4), a popular approach is the use of penalty functions, such as LASSO, see tibshirani1996regression, or SCAD, see fan2001variable, to shrink unimportant regression coefficients to zero. In actuarial literature, jeong2020non used non-convex regularization methods in order to obtain stable estimation of loss development factors in insurance claims reserving. In multi-type variable setting, devriendt2020sparse is currently the only paper which considers multi-type feature selection under a Poisson regression framework for claim frequency modeling.

While the existing literature address some of the aforementioned claim severity modeling needs, we are still lacking a universal modeling framework, which not only provides versatility to fit a multimodal heavy-tailed severity distribution, but also explains the covariates' influence on multiple distributional parts with variable selections. As a result, the goal of this paper is to integrate, adapt and extend the existing modeling techniques, and devise a universal insurance claim severity modeling framework which simultaneously address all of the four modeling needs mentioned above. To this end, we make the following contributions:

Firstly, we introduce a mixture composite regression model for approximating claim severities based on the use of available covariate information. This extends the setup of reynkens2017modelling, who used a finite mixture distribution for the body and a Pareto-type distribution for the tail without using covariates, by incorporating covariates impacts on all three parts of the severity distribution: clustering probabilities, body part and tail part.

Secondly, we propose a group-fused regularization approach for variable selection. This approach allows us to select three different sets of variables which significantly impact the previously mentioned three parts of the claim severity distribution respectively. The set of variables chosen is homogeneous across all mixture components to preserve model interpretability. Furthermore, this approach enables regularization under multi-type variable settings.

Thirdly, we develop an asymptotic theory for the regularization approach. The following two main results theoretically justify the appropriateness of the proposed method: (i) The proposed method is consistent in terms of feature selection, in particular, as sample size goes to infinity, the proposed method will correctly merge and shrink regression coefficients across the various modeling parts. (ii) The parameters of the reduced model, after merging and shrinking the regression coefficients, are asymptotically normal with zero mean, and their variances are the same as the parameter uncertainties obtained by fitting the same mixture composite regression model to the reduced model (e.g. mean claim severity). The implication of the above two main results is that we can construct Wald-type confidence intervals and Efron percentile bootstrap confidence intervals of model parameters and other quantities of interest.

Finally, we present a novel Generalized Expectation-Maximization (GEM) algorithm for estimating the parameters of the proposed model with parameter regularization. The GEM algorithm is demonstrated to perform satisfactorily when the mixture-Gamma Lomax composite regression model is fitted to a Greek MTPL dataset which inherits all the previously described features.

The remainder of this article proceeds as follows. In Section 2, we introduce the framework of the mixture composite regression model. Section 3 presents the feature selection approach which can be used for selecting important variables for explaining the claim severity distribution in the presence multi-type covariates. In Section 4 we provide the theoretical foundations, such as consistency and asymptotic normality, upon which the proposed feature selection approach is based for merging and shrinking parameters correctly with high probability when the sample size is large. Furthermore, we develop Wald type and bootstrap two-sided confidence intervals for the parameters. The maximum likelihood estimation (MLE) procedure for our proposed model via the GEM algorithm is presented in Section 5. In Section 6, we describe the MTPL dataset that we use for our empirical analysis, and provide estimation and model comparison for various benchmark distributions. In Section 7, we fit the proposed mixture composite distribution and subsequently the mixture composite regression model with feature regularizations. Concluding remarks are given in Section 8, and other miscellaneous details are included in the Supplementary materials.

Modeling framework

This section summarizes the features that are incorporated in a regression modeling framework for addressing the challenges encountered in claim severity datasets in general insurance. In particular, motivated by the characteristics of the multimodal and heavy-tailed Greek MTPL insurance dataset studied below, we propose the following mixture composite regression model.

Let $Y\in\mathbb{R}^{+}$ be the claim severity random variable, and let $\bm{x}\in\mathbb{R}^{D}$ be the vector of covariate information\footnote{ Note also that all vectors are assumed to be column vectors.}. The density of the mixture composite regression model is given by

where $\pi_j(\bm{x};\bm{\alpha})$, $1\le j \le g+1$, are covariate-dependent component weights given by

with $\bm{\alpha}_{g+1}=\bm{0}$ for model identifiability, and $\bm{\alpha}=(\bm{\alpha}_1,\ldots,\bm{\alpha}_{g+1})\in\mathbb{R}^{D\times (g+1)}$. $f$ and $h$ are the body and tail density functions respectively, such that the first $g$ mixture components are specialized in capturing small to moderate claim amounts while the last component focuses on extreme claims. $F$ and $H$ are the corresponding cdfs.

In this paper, we specify $f$ and $h$ as Gamma (body) and Lomax (tail, also called Pareto type II) density functions given by, respectively,

and

The choice of Gamma density is motivated by its light-tailed and uni-modal characteristics to capture small to moderate claims. Also, mixture of Gammas provides sufficient flexibility to capture complex distributional structures like multimodality, thanks to the deneness property of Gamma mixture. The choice of the Lomax density for the tail is motivated by its polynomial tail characteristics with the tail index $\eta$ describing the tail-heaviness of the distribution. The analytical form of the truncated Lomax distribution also makes the model estimation procedures computationally desirable. Note however that one may choose other plausible model specifications as long as $f$ is a unimodal light-tailed distribution while $h$ is a heavy-tailed distribution. To avoid distorting the focus of this paper and given that the fitting results of the real dataset (Section (ref)) are satisfactory, we leave the comparisons among various model specifications as a future research direction.

Furthermore, $\bm{\beta}=(\bm{\beta}_1,\ldots,\bm{\beta}_g)\in\mathbb{R}^{D\times g}$ and $\bm{\nu}\in\mathbb{R}^{D}$ are the regression coefficients for the body and tail distributions, respectively. The proposed distribution is characterized by a splicing threshold $\tau>0$, which is predetermined using expert opinion via performing e.g. extreme value analysis instead of treated it as a parameter estimated by a likelihood approach; this is mainly motivated by estimation stability and is adopted by e.g. reynkens2017modelling.

The mean of $Y|\bm{x}$ is given by

The composite model in Equation ((ref)) can alternatively be regarded as a mixture of $g$ right truncated Gamma distributions for the body and a left truncated Lomax distribution for the tail. Each claim is classified to one of the $g+1$ subgroups ($g$ subgroups for body and one subgroup for tail) with probabilities $\{\pi_j(\bm{x};\bm{\alpha})\}_{j=1,\ldots,g+1}$, where each subgroup may correspond to a different claim sub-type. Regression coefficients $\bm{\alpha}$ explains the heterogeneity of the assignment probabilities across different claims, while $\bm{\beta}$ explain the systematic effects of the claims within a given subgroup. The regression coefficients $\bm{\nu}$ for the tail distribution, on the other hand, capture the effect of covariates to the tail-heaviness of claims.

The motivation of introducing a composite model in Equation ((ref)) instead of using a mixture-Gamma Lomax model is that there are no overlapping density regions between the body and tail distributions under the proposed framework. We will show in our motivating application in Section (ref) that this results in a more robust and stable estimation of the tail index, since it is not distorted by attritional claims from the body of the distribution. One should however note that the mixture probabilities connect tail and body regression parameter estimation, i.e., the proposed composite regression model does not decouple into independent estimation parts.

Feature selection method

In this section, we propose the group fused penalty approach to select the variables to describe the systematic effects in claim severities under a multi-type covariates setting. We will select three potentially different sets of variables that may influence, respectively, the subgroup probabilities, body and tail of the distribution. For the sake of model interpretability, we select the same set of variables for all mixture components of the body of the data and the mixing probabilities.

Suppose there are $n$ independent claims $\bm{Y}=(Y_1,\ldots,Y_n)^T$, and denote their realizations by $\bm{y}=(y_1,\ldots,y_n)^T$. For each claim $i=1,\ldots,n$, we have a covariates vector $\bm{x}_i=(x_{i1},\ldots,x_{iD})^T\in \mathbb{R}^D$ with $x_{i1}=1$ (for the intercept component). Define $\bm{X}=(\bm{x}_1,\ldots,\bm{x}_n)^T\in \mathbb{R}^{n\times D}$ as the design matrix containing the covariate information of all $n$ observations. The observed data log-likelihood is given by

where $\bm{\Phi}$ contains all model parameters. To incorporate variable selection, we propose a group fused regularization approach, where the penalty function for the regression parameters is as follows

with $P_{\bm{\lambda}_1,n}(\bm{\alpha})$, $P_{\bm{\lambda}_2,n}(\bm{\beta})$ and $P_{\bm{\lambda}_3,n}(\bm{\nu})$ being the penalty functions on the regression parameters $\bm{\alpha}\in\mathbb{R}^{D\times (g+1)}$, $\bm{\beta}\in\mathbb{R}^{D\times g}$ and $\bm{\nu}\in\mathbb{R}^{D}$. These are given by

where $\|\cdot\|_2$ is the $L^2$-norm, $\lambda_{1kn}$, $\lambda_{2kn}$ and $\lambda_{3kn}$ are penalty tuning parameters, $p_{1n}$, $p_{2n}$ and $p_{3n}$ are concave non-decreasing penalty functions (which will be chosen proportional to the sample size $n$) and $K_1$, $K_2$ and $K_3$ correspond to the numbers of penalization terms. Finally, $\{\bm{c}_{lk}\}_{l=1,2,3}$ are predetermined vectors of penalty coefficients which allow for different types of penalties, including standard LASSO to shrink continuous variables, fused LASSO to merge regression coefficients of various ordinal categorical variables, and generalized fused LASSO to merge regression coefficients for nominal categorical variables. For a full description on constructing predetermined vectors for each type of variables (continuous, ordinal categorical and nominal categorical), we refer the reader to oelker2017uniform, in the statistics literature, and devriendt2020sparse in the actuarial literature.

The aim is to maximize the following objective function (penalized log-likelihood)

We will use the following two commonly used penalty functions (for $l\in\{1,2,3\}$ and $\psi\geq 0$) to illustrate the usefulness of the proposed feature selection method:

Here, we set $n_1=n$ for the total number of observations, $n_2:=n_b=\sum_{i=1}^n1\{y_i\leq\tau\}$ for the number of observations in the body, and $n_3:=n_t=\sum_{i=1}^n1\{y_i>\tau\}$ for the number of observations in the tail.

Each of the parameters $\bm{\alpha}$ and $\bm{\beta}$ contains $g$ sets of regressors (one for each mixture component of the body, we initialize $\bm{\alpha}_{g+1}=\bm{0}$). For the sake of interpretability, the proposed group regularization method shrinks and merges regression coefficients of any variable uniformly across all mixture components, and the sets of variables does not vary across mixture components. Therefore, the proposed method allows us for choosing three sets of variables which significantly impact each of the three modeling parts -- the subgroup probabilities, the body part and the tail part of the severity distribution.

Asymptotic properties

This section presents two asymptotic theorems regarding to the proposed mixture composite model with the feature selection method. Our motivation is two-fold: First, we want to theoretically justify the ability of the proposed feature selection approach in correctly merging and shrinking regression coefficients. Second, the theorems provide guidance to estimate model uncertainty. We will only present the key results and discuss their implications in this section. All construction details, including the assumptions and proof details, are postponed to Appendix (ref). Suppose $Y_i$, given $\bm{x}_i$, is generated by the model of Equation ((ref)) with true model parameter $\bm{\Phi}_0=(\bm{\alpha}_0,\bm{\beta}_0,\bm{\phi}_0,\theta_0,\bm{\nu}_0)$. We first have the following theorem:

As discussed in Appendix (ref), both LASSO and SCAD penalty functions can be constructed to satisfy all assumptions H1 to H4. Therefore, the above theorem says that the estimated parameters $\hat{\Phi}_n$ under the proposed model setup will converge to the true model parameters as $n\rightarrow\infty$.

Apart from consistency, it is important to show sparsity of the proposed feature selection method, enabling consistent variable selection. To do so, we need to linearly transform the parameter space and formulate the asymptotic properties in the transformed space. We first define $\bm{C}_l=(\bm{c}_{l1},\ldots,\bm{c}_{lK_l})$ as a design matrix of penalty coefficients. Further, denote $\mathcal{Z}_1=\{k:\|\bm{c}_{1k}^T\bm{\alpha}_0\|_2=0\}$, $\mathcal{Z}_2=\{k:\|\bm{c}_{2k}^T\bm{\beta}_0\|_2= 0\}$ and $\mathcal{Z}_3=\{k:|\bm{c}_{3k}^T\bm{\nu}_0|= 0\}$, representing the regression coefficients to be merged or shrinked. W.l.o.g., we hereafter assume that under the true model, $\mathcal{Z}_l=\{1,2,\ldots,s_l\}$ for $l=1,2,3$, and we construct a reduced matrix $\bar{\bm{C}}_{\text{red},l}:=(\bm{c}_{l1},\ldots,\bm{c}_{lm_l})$ of the linearly independent vectors which span the space of the vectors $\{\bm{c}_{l1},\ldots,\bm{c}_{ls_l}\}$. Note that we always have $m_l\leq\min(s_l,P)$. Further, we construct linearly independent vectors $\bar{\bm{C}}_{\text{ind},l}:=(\bm{c}^*_{l,m_l+1},\ldots,\bm{c}^*_{l,D})$ which are also linearly independent of all vectors in $\bar{\bm{C}}_{\text{red},l}$. Then, define $\bar{\bm{C}}_l=(\bar{\bm{C}}_{\text{red},l},\bar{\bm{C}}_{\text{ind},l})$ which is a $D\times D$ full rank matrix, and define the transformed parameters $\bm{\alpha}^*=\bar{\bm{C}}_1^T\bm{\alpha}$, $\bm{\beta}^*=\bar{\bm{C}}_2^T\bm{\beta}$ and $\bm{\nu}^*=\bar{\bm{C}}_3^T\bm{\nu}$. Note that the transformed parameters can also be decomposed as $\bm{\alpha}^*=({\bm{\alpha}^{*}_{\text{red}}}^T,{\bm{\alpha}^{*}_{\text{ind}}}^T)^T$, $\bm{\beta}^*=({\bm{\beta}^{*}_{\text{red}}}^T,{\bm{\beta}^{*}_{\text{ind}}}^T)^T$ and $\bm{\nu}^*=({\bm{\nu}^{*}_{\text{red}}}^T,{\bm{\nu}^{*}_{\text{ind}}}^T)^T$, where $\bm{\alpha}^{*}_{\text{red}}=\bar{\bm{C}}_{\text{red},1}^T\bm{\alpha}$, $\bm{\alpha}^{*}_{\text{ind}}=\bar{\bm{C}}_{\text{ind},1}^T\bm{\alpha}$, $\bm{\beta}^{*}_{\text{red}}=\bar{\bm{C}}_{\text{red},2}^T\bm{\beta}$, $\bm{\beta}^{*}_{\text{ind}}=\bar{\bm{C}}_{\text{ind},2}^T\bm{\beta}$, $\bm{\nu}^{*}_{\text{red}}=\bar{\bm{C}}_{\text{red},3}^T\bm{\nu}$ and $\bm{\nu}^{*}_{\text{ind}}=\bar{\bm{C}}_{\text{ind},3}^T\bm{\nu}$. The above mathematical constructions allow us to re-write the penalized log-likelihood as a function of the transformed parameters $\bm{\Phi}^*:=(\bm{\alpha}^*,\bm{\beta}^*,\bm{\phi},\theta,\bm{\nu}^*)$ as follows:

with log-likelihood $\mathcal{L}_n^*(\bm{\Phi}^*)=\mathcal{L}_n(\bm{\Phi})$ and penalty term

where $\tilde{\bm{c}}_{lk}=\bar{\bm{C}}_l^{-1}\bm{c}_{lk}$, for $l=1,2,3$. Suppose that the true model parameters are given by $\bm{\Phi}^{*}_0:=(\bm{\alpha}^*_0,\bm{\beta}^*_0,\bm{\phi}_0,\theta_0,\bm{\nu}^*_0)$. This can be decomposed $\bm{\Phi}^{*}_0:=(\bm{\Phi}^{*}_{\text{red},0},\bm{\Phi}^{*}_{\text{ind},0})$, with $\bm{\Phi}^{*}_{\text{red},0}:=(\bm{\alpha}^*_{\text{red},0},\bm{\beta}^*_{\text{red},0},\bm{\nu}^*_{\text{red},0})$ and $\bm{\Phi}^*_{\text{ind},0}:=(\bm{\alpha}^*_{\text{ind},0},\bm{\beta}^*_{\text{ind},0},\bm{\phi}_0,\theta_0,\bm{\nu}^*_{\text{ind},0})$. Note that by construction $\bm{\Phi}^{*}_{\text{red},0}=\bm{0}$. Finally, denote $\hat{\bm{\Phi}}^*_n:=(\hat{\bm{\Phi}}^*_{\text{red},n},\hat{\bm{\Phi}}^*_{\text{ind},n})$ as the corresponding estimator of model parameters. We have the following theorem, which is an extension of the oracle property given by fan2001variable:

The above theorem shows that when the sample size is large, the proposed feature selection method merges and shrinks parameters correctly with high probability. Moreover, the estimated parameters of the reduced model are asymptotically normal. While detailed discussions are leveraged to Remark (ref) of Appendix (ref), the adjustment and bias terms ${\mathcal{P}^*}''(\bm{\Phi}^*_{\text{ind},0})$ and ${\mathcal{P}^*}'(\bm{\Phi}^*_{\text{ind},0})$ are both asymptotically negligible when the penalty function is LASSO with an adaptive approach (which will be discussed in Section (ref)) or SCAD. For large sample sizes, the estimated parameters are approximately unbiased and we may approximate the variance of the estimated transformed parameters as

where $\widehat{\mathcal{I}}_{\text{ind}}^{*}(\hat{\bm{\Phi}}^*_{\text{ind},n})$ is the sample Fisher information of the reduced model. In other words, parameter uncertainty under the proposed feature selection method is equivalent to that under the reduced model after selecting the variables. With this regards, the construction of confidence intervals (CI) of parameters is straightforward:

Model estimation

Direct optimization of the penalized log-likelihood in Equation ((ref)) is difficult. Firstly, the log-likelihood $\log h_{Y}(y_i;\bm{\alpha},\bm{\beta},\bm{\phi}, \theta,\bm{\nu},\bm{x}_i)$ is the logarithm of a sum of $(g+1)$ mixture terms. Secondly, observe that $\log h_{Y}(y_i;\bm{\alpha},\bm{\beta},\bm{\phi}, \theta,\bm{\nu},\bm{x}_i)$ contains distribution function $F(\tau;\exp\{\bm{\beta}_j^T\bm{x}_i\},\phi_j)$, which is not available in closed form; this is not the case for the Lomax distribution since $H(\tau;\theta,\exp\{\bm{\nu}^T\bm{x}_i\})$ has an analytical form. Thirdly, the penalty functions are not continuously differentiable.

Model estimation of an ordinary splicing model is typically simple, thanks to the non-overlapping density parts between the body and tail, so that one may factor out the likelihood function and separately calibrate the three parts of distribution -- subgroup probability, body and tail. Nonetheless, under the regression framework outlined in Equation ((ref)) with variable selection techniques embedded, the weight regression parameters $\bm{\alpha}$ share and interact across all $g$ body components and one tail component. Therefore, there is no straight-forward way to segregate the likelihood function and simplify the estimation procedure.

Motivated by the aforementioned computational challenges, this section presents the strategy to estimate the parameters and select important variables under the proposed modeling framework.

Construction of complete data

We first construct a hypothetical complete dataset. We present a modified version of the method introduced by FUNG2020MoECensTrun. Define the complete data

with three extra elements defined as follows:

We assume that the cases $(Y_i,\bm{Z}_i,\bm{K}_i,\{\bm{Y}_{ij}'\}_{j=1,\ldots,g})$ are independent in $1\le i \le n$. Moreover, we assume independence between $\bm{Z}_i$, $\bm{K}_i$, $Y_i$ and $\bm{Y}'_{ij}$, that $Y_{ij1}',\ldots,Y_{ijk_{ij}}'$ are i.i.d. given the covariates $\bm{x}_i$, and that $Y_i$ and $\bm{Y}'_{ij}$ are subgroup conditionally identically distributed. Furthermore, the components $K_{ij}$ of $\bm{K}_i$ are independent and follow the geometric distribution

The complete data log-likelihood function is given by

This is easier to evaluate and optimize compared to Equation ((ref)) given that $H$ has an analytical form, which is the case for the Lomax distribution. The complete data penalized log-likelihood is given by

The GEM algorithm

Parameter estimation is conducted using a Generalized Expectation-Maximazation (GEM) algorithm, where in the M-step is a modified version of the penalized iteratively re-weighted least squares (PIRLS) method proposed by oelker2017uniform, where this method is also technically justified. Following oelker2017uniform, due to the non-differentiability of the penalty function, we perturb the penalty function in Equation ((ref)) as follows

where the $\epsilon$-perturbed penalty functions $P_{\bm{\lambda}_1,n,\epsilon}(\bm{\alpha})$, $P_{\bm{\lambda}_2,n,\epsilon}(\bm{\beta})$ and $P_{\bm{\lambda}_3,n,\epsilon}(\bm{\nu})$ are given by

with $\|\bm{w}\|_{2,\epsilon}=\left(\bm{w}\bm{w}^T+\epsilon\right)^{1/2}$ and $|w|_{\epsilon}=(w^2+\epsilon)^{1/2}$ for any vector $\bm{w}$ and scalar $w$. In the following, instead of maximizing Equation ((ref)), we maximize the $\epsilon$-perturbed complete data penalized log-likelihood given by

Now, $\mathcal{F}^{\text{com}}_{n,\epsilon}(\bm{\Phi})$ is continuously differentiable w.r.t. any parameter and, hence, it is more computationally tractable. Furthermore, note that $\mathcal{P}_{n,\epsilon}(\bm{\Phi})\rightarrow\mathcal{P}_n(\bm{\Phi})$ and hence $\mathcal{F}_{n,\epsilon}^{\text{com}}(\bm{\Phi})\rightarrow\mathcal{F}^{\text{com}}_n(\bm{\Phi})$ as $\epsilon\rightarrow 0$, so choosing a very small $\epsilon>0$, the perturbation of the estimated parameters $\bm{\Phi}$ will be negligible. In simulation studies and real data analysis, we find that the choice of $\epsilon=10^{-10}$ works well.

E-step

In the $l^{\text{th}}$ iteration, the expectation of the complete data $\epsilon$-perturbed penalized log-likelihood is computed as follows:

where the updated quantities $z_{ij}^{(l)}$, $k_{ij}^{(l)}$, $\log \tilde{f}(\widehat{y_{ij}'}^{(l)},\widehat{\log y_{ij}'}^{(l)};\exp\{\bm{\beta}_j^T\bm{x}_i\},\phi_j)$, $\widehat{y_{ij}'}^{(l)}$ and $\widehat{\log y_{ij}'}^{(l)}$ are displayed in Equations (2.2) to (2.6) in the supplementary materials.

Note that $\widehat{\log y_{ij}'}^{(l)}$ is represented by a numerical integral, and hence, in general, it does not have an analytical solution. Here, we adopt a Stochastic EM approach, where for each $i=1,\ldots,n$ and $j=1,\ldots,g$ we simulate $\log Y_{ij}'^{(l)}$ from the conditional density of $Y_{ij}'^{(l)}$ given by

M-step

In this step, we attempt to find a parameter update $\bm{\Phi}^{(l)}$ in such that we will receive a monotonicity $Q_{\epsilon}(\bm{\Phi}^{(l)};\bm{y},\bm{X},\bm{\Phi}^{(l-1)})\geq Q_{\epsilon}(\bm{\Phi}^{(l-1)};\bm{y},\bm{X},\bm{\Phi}^{(l-1)})$. While $Q_{\epsilon}(\bm{\Phi};\bm{y},\bm{X},\bm{\Phi}^{(l-1)})$ is now differentiable w.r.t. any parameter, direct implementation of an iteratively re-weighted least squares (IRLS) algorithm is challenging due to concavity of the penalty functions. In this section, we propose the use of convex quadratic approximation to the penalty functions, analogously to fan2001variable and oelker2017uniform, such that the implementation of an IRLS algorithm is feasible. We approximate $\mathcal{P}_{n,\epsilon}(\bm{\Phi})$ by

where

The properties below justify the use of such approximations:

Next, we decompose $\tilde{Q}_{\epsilon}(\bm{\Phi};\bm{y},\bm{X},\bm{\Phi}^{(l-1)})$ into the following terms:

where

Update of $\bm{\alpha}^{(l-1)}$ to $\bm{\alpha}^{(l)}$ can be done by sequentially adopting the IRLS approach for $j=1,\ldots,g$:

where the derivatives are presented in Equations (2.7) and (2.8) of the supplementary material which are expressed in analytical forms.

Similarly, update of $\bm{\beta}^{(l-1)}$ to $\bm{\beta}^{(l)}$ using IRLS involves

with the analytical forms of derivatives given by Equations (2.9) and (2.10) of the supplementary material.

After updating $\bm{\beta}$, we update $\phi_j^{(l-1)}$ to $\phi_j^{(l)}$ directly using function optimize in R, which is found to involve little computational burden compared to the IRLS procedures above:

After that, the same IRLS procedure leads to an update of $\bm{\nu}^{(l-1)}$ to $\bm{\nu}^{(l)}$:

with the analytical forms of derivatives given by Equations (2.11) and (2.12) of the supplementary material.

Finally, $\theta$ can be updated directly using the optimize function or the Newton-Raphson method, aiming to achieve

Because of Corollary (ref), the M-step ensures $Q_{\epsilon}(\bm{\Phi}^{(l)};\bm{y},\bm{X},\bm{\Phi}^{(l-1)})\geq Q_{\epsilon}(\bm{\Phi}^{(l-1)};\bm{y},\bm{X},\bm{\Phi}^{(l-1)})$ given a very small $\epsilon>0$. The GEM algorithm is iterated until the observed data $\epsilon$-perturbed penalized log-likelihood is improved by less than a threshold $10^{-2}$ or the maximum number of iterations of 200 is reached.

Initialization procedures

Initialization of parameters $\bm{\Phi}^{(0)}$ can be done using the clusterized method of moments (CMM) approach proposed by gui2018fitting. It requires a $K$-means clustering method to assign observations $y_i$ with $y_i\leq\tau$ to one of the $g$ subgroups for the body, and observations $y_i$ with $y_i>\tau$ to the tail component. Then, we set initial parameters which match the first two moments for each mixture component, after initially fixing all the regression parameters be zero except for the intercepts. Refer to e.g. Section 3.3.3 of FUNG2020MoECensTrun for more details.

Choice of the number of mixture components

Usually, the choice of the number of mixture components $g$ of the body can be determined based on standard specification criteria, including Akaike's Information Criterion (AIC) and the Bayesian Information Criterion (BIC). However, for our motivating dataset, which will be described in Section 6, below, the AIC and BIC criteria would both lead to an excessively large number of components which would significantly impede the model interpretability. The main reason for obtaining a large $g$ is that the claim severity distribution has many small nodes for smaller claim amounts (i.e. less than $10,000$), as evidenced by Figure (ref) in Section 6. Excessive fitting and modeling of such smaller claim amounts does not bring much insight from an insurance ratemaking perspective because such smaller claims could even be modelled by an empirical distribution. As a result, for this particular dataset, we adopt a qualitative method, which chooses $g$ as the minimum number of components required for the proposed model to capture all nodes above a claim severity threshold of $10,000$.

Selection of variables

The proposed GEM algorithm with group fused penalty functions shrinks some regression coefficients to zero and merges some coefficients across different levels of a categorical variable. Being a variant of the PIRLS approach, the proposed algorithm also caters for a wide range of concave penalty functions. However, as pointed out by devriendt2020sparse, the parameters obtained by the proposed algorithm are not exact. Therefore, in order to select the variables and reduce model complexity, after every model fit we need to perform an automatic adjustment algorithm to remove parameters very close to zero and merge the parameters when their values are very close to each other.

Denote the fitted model parameter as $\hat{\bm{\Phi}}=(\hat{\bm{\alpha}},\hat{\bm{\beta}},\hat{\bm{\phi}},\hat{\theta},\hat{\bm{\nu}})$. Further, with slight abuse of notation, denote $\hat{\bm{\alpha}}_{p}$ as the $p^{\text{th}}$ row vector of $\hat{\bm{\alpha}}$, as opposed to $\hat{\bm{\alpha}}_{j}$ as the $j^{\text{th}}$ column vector of $\hat{\bm{\alpha}}$. Similarly, denote $\hat{\bm{\beta}}_{p}$ as the $p^{\text{th}}$ row vector of $\hat{\bm{\beta}}$. Also, let $\hat{z}_{ij}$, $\hat{k}_{ij}$, $\widehat{y_{ijk}'}$ and $\widehat{\log y_{ijk}'}$ be the $z_{ij}^{(l)}$, $k_{ij}^{(l)}$, $\widehat{y_{ijk}'}^{(l)}$ and $\widehat{\log y_{ijk}'}^{(l)}$ obtained by the E-step using the fitted parameters. Define the partial log-likelihood functions $S(\bm{\alpha})$, $T(\bm{\beta},\bm{\phi})$ and $V(\theta,\bm{\nu})$, respectively, for the mixing probabilities, body distributions and tail distribution analogously to Equations ((ref)) to ((ref)) as follows:

The general principle of the automatic adjustment algorithm is to fine tune the regression parameters, so the regression parameters (which are close to zero or very close to each other) are shrinked or merged in exact. Since fine tuning of parameters would lead to another source of error, the automatic adjustment algorithm needs to ensure that the likelihood-based quantities displayed above would not change significantly due to fine-tuning. We leverage the step-by-step algorithm to Section 2.2 of the supplementary materials.

Tuning of hyperparameters

The remaining problem is to select appropriate tuning parameters $\bm{\lambda}:=(\bm{\lambda}_1,\bm{\lambda}_2,\bm{\lambda}_3)$ which control the model complexity and hence select variables useful for explaining different parts of the claim severity distributions. The current theory in Section (ref) only provides guidance on the order of $\bm{\lambda}$, but in application it is obvious that grid search on $\bm{\lambda}$ is computationally prohibitive because of the curse of dimensionality. As a result, we adopt an adaptive-standardization approach similar to devriendt2020sparse, where we restrict $\lambda_{1kn}=w_{1k}\lambda_1$, $\lambda_{2kn}=w_{2k}\lambda_2$ and $\lambda_{3kn}=w_{3k}\lambda_3$. Here,

where $w_{1k}^{(\text{ad})}=\|\bm{c}_{1k}^T\hat{\bm{\alpha}}\|_2^{-1}$, $w_{2k}^{(\text{ad})}=\|\bm{c}_{2k}^T\hat{\bm{\beta}}\|_2^{-1}$ and $w_{3k}^{(\text{ad})}=|\bm{c}_{3k}^T\hat{\bm{\nu}}|^{-1}$ are the adaptive terms, and $w_k^{(\text{st})}$ is the standardization term. Note here that the estimated parameters $\hat{\bm{\alpha}}$, $\hat{\bm{\beta}}$ and $\hat{\bm{\nu}}$ are obtained on the fitting procedures obtained in Section (ref) starting with very small tuning parameters $\bm{\lambda}$ (or even $\bm{\lambda}=0$). $p_1$ and $p_2$ are the two categories that the $k^{\text{th}}$ penalty term is attempting to merge for categorical variables, and $(n_{p_1},n_{p_2})$ are the number of observations being classified to those respective categories. $p_{\bm{G}}$ is the number of categories for the respective explanatory variable, while $r_{\bm{G}}$ is the number of penalty terms for the respective explanatory variable. Note that $r_{\bm{G}}=p_{\bm{G}}-1$ for ordinal variables and $r_{\bm{G}}=p_{\bm{G}}(p_{\bm{G}}-1)/2$ for nominal variables. For continuous variables, we set $w_k^{(\text{st})}=1$ instead. The adaptive weights facilitate more efficient shrinkage or merger of regression coefficients, achieving the oracle property presented by zou2006adaptive. The standardization weights, on the other hand, address the issues of level imbalances and imbalances of numbers of terms on an explanatory variable involved in the penalty functions.

After specifying the weights, we can perform a grid search on $(\lambda_1,\lambda_2,\lambda_3)$ to find optimal tuning parameters. Motivated by the likelihood-based deviance approach by khalili2010new, we propose the following method, which allows us doing the separate grid search for $\lambda_1$, $\lambda_2$ and $\lambda_3$.

After obtaining the fitted model parameters $\bm{\Phi}$ starting with a small $\bm{\lambda}$, we compute the estimated latent variables $\hat{z}_{ij}$, $\hat{k}_{ij}$, $\widehat{y_{ijk}'}$ and $\widehat{\log y_{ijk}'}$ outlined in Section (ref) and assume that they are fixed during the whole process of grid searching. Then, for each $\lambda_1$, $\lambda_2$ and $\lambda_3$ within separate specified (one-dimensional) grids, we refit the models by maximizing the (unpenalized) partial log-likelihood functions $S_0(\bm{\alpha})$ in Equation ((ref)) (which only requires iterating Equations ((ref))), $T_0(\bm{\beta},\bm{\phi})$ in Equation ((ref)) (iterating Equations ((ref)) and ((ref))) and $V_0(\theta,\bm{\nu})$ in Equation ((ref)) (iterating Equations ((ref)) and ((ref))). We denote the resulting fitted parameters as $\hat{\bm{\alpha}}(\lambda_1)$, $(\hat{\bm{\beta}}(\lambda_2), \hat{\bm{\phi}}(\lambda_2))$ and $(\hat{\theta}(\lambda_3),\hat{\bm{\nu}}(\lambda_3))$. This avoids repeating the whole GEM procedure over a multidimensional grid of $(\lambda_1,\lambda_2,\lambda_3)$ which is computationally prohibitive. Optimal $(\lambda_1,\lambda_2,\lambda_3)$ can be determined by various choices of criteria, where in this paper we will present partial AIC (pAIC), partial BIC (pBIC) and $K$-fold cross-validation (CV). For pAIC or pBIC approach, we define

where $\mathcal{N}_1(\lambda_1)$, $\mathcal{N}_2(\lambda_2)$ and $\mathcal{N}_3(\lambda_3)$ are the effective number of parameters (i.e. the number excluding zeroes and redundant parameter values) for $\hat{\bm{\alpha}}(\lambda_1)$, $(\hat{\bm{\beta}}(\lambda_2),\hat{\bm{\phi}}(\lambda_2))$ and $(\hat{\theta}(\lambda_3),\hat{\bm{\nu}}(\lambda_3))$, respectively. Recall that $n_b$ and $n_t$ are the number of observations allocated to body and tail components, respectively. Now, $\lambda_1$, $\lambda_2$ and $\lambda_3$ can be chosen once at a time through minimizing the pAICs or pBICs.

For $K$-fold CV, we partition the data into $K$ disjoint folds and measure the performance on each fold after training the remaining $K-1$ folds. The performance metric we use in this paper is called “partial deviance" given by

For $l=1,2,3$, the optimal $\hat{\lambda}_l$ is the largest one such that the corresponding partial deviance is within one standard deviation of its minimum.

As the regularization functions make the estimated parameters of the fitted model biased towards zero, it is important to collapse the regression parameters $(\hat{\bm{\alpha}}(\hat{\lambda}_1),\hat{\bm{\beta}}(\hat{\lambda}_2),\hat{\bm{\nu}}(\hat{\lambda}_3))$ and covariate matrix $\bm{X}$, and re-estimate the model where the penalties are excluded (i.e. $\bm{\lambda}=\bm{0}$), using the full GEM algorithm outlined by Section (ref). This will also update the estimated latent variables for better accuracy. A related approach can be found by devriendt2020sparse.

Motivating dataset: Greece MTPL claim amounts

In this section, we present a dataset which motivates the proposed modeling and feature selection framework described in above. The characteristics of the dataset are first described in Section (ref). Then, we fit some state-of-the-art models to the dataset in Section (ref) to show the necessity of adopting the proposed modeling framework.

Data description

The dataset for our study was kindly provided by a major insurance company operating in Greece. It consists of 64,923 motor third-party liability (MTPL) insurance policies with non-zero property claims for underwriting years 2013 to 2017. The sample comprised of policyholders with complete records; i.e., with the availability of all explanatory variables under consideration, and with at least one reported accident over the five underwriting years. These explanatory variables are summarized in Table (ref).

An exploratory analysis was carried out in order to identify the challenges that need to be surmounted for efficiently modeling these property damage claim costs based on the subset of explanatory variables with the highest predictive power. Firstly, as we observe from Figures (ref) and (ref), the empirical claim severity distribution is multimodal and heavy-tailed. In particular, the empirical density plot of the claim amounts in the left panel of Figure (ref) shows that there are at least three major nodes or clusters in the empirical density function: one for small claim severities of $<$10,000, one for claim severities of about 30,000 and one for claim severities of about 80,000--100,000. Additionally, the density for the log claim amounts in the right panel of Figure (ref) reveals even more complex distribution characteristics illustrated by many small peaks of the density function, especially for small claim sizes. Furthermore, regarding the heavy-tailed nature of the data, the log-log plot in the left panel of Figure (ref) seems asymptotically linear (an asymptotic red straight line is fitted) with slope of roughly $-1.3$ (this represents the tail index $\alpha$ of the empirical distribution). Furthermore, the mean excess plot which is depicted in the right panel of Figure (ref) appears linear when the claim size exceeds a threshold of around 270,000 (black vertical line), with asymptotic slope of 3.35 which also suggests $\alpha\approx 1.3$. Secondly, as far as the explanatory power of the variables is concerned, we studied the influence of the explanatory variables to the claim amounts through plotting the empirical density plots across each level of each variables. The results for the variables Driver's age, Insurance duration, Payment way and Policy type are displayed in Figure (ref). From Figure (ref), we see that some variables have some apparent effects on the peak (or probability) of each cluster instead of the position of each cluster. For example, from the bottom right panel, the policy with “expensive type" has a higher probability assigned to the tail cluster and lower probabilities assigned to the remaining body clusters. Finally, it should be noted that considering all explanatory variables to be categorical, the 10 explanatory variables lead to 137 covariates in total. Therefore, since, as was previously mentioned, the impact of covariates on the claim severity distribution could be multi-fold (e.g. covariates may affect cluster assignment probabilities, average claim severity given a particular cluster and/or tail-heaviness), an appropriate regression model for this Greek MTPL dataset should contain multiple regressors. Obviously, this will lead to a large number of parameters without parameter regularizations which can potentially result in an over-fitting problem and impede model interpretations. Therefore, these issues outline the importance of variable selection.

Preliminary model fitting

In this subsection we first explore probability distributions which may appropriately fit the distribution of claim amounts, ignoring the effects of covariates. As we observed multiple nodes and heavy-tailed characteristics of the claim amount distribution, it is natural to consider a finite mixture model with both light- and heavy-tailed mixture components to capture such characteristics. Motivated by BLOSTEIN201935 who propose a finite mixture of various classes of distributions, a plausible benchmark is a mixture-Gamma Lomax model where the claim severity $Y$ is modelled by a density function

where $\bm{\pi}=(\pi_1,\ldots,\pi_{g+1})$ are the mixture probabilities, $\bm{\mu}=(\mu_1,\ldots,\mu_{g})$ and $\bm{\phi}=(\phi_1,\ldots,\phi_g)$ are the mean and dispersion parameters. $f$ is the Gamma density function for modeling the body, and $h$ is the Lomax density function for modeling the tail given by Equations ((ref)) and ((ref)) respectively.

Observing three major nodes in the density shown in the left panel of Figure (ref), we first start with the above model with $g=3$ components for the body. Summary statistics are shown in Table (ref), which compares the above model to several other classical unimodal models, including the Gamma (GA), Weibull (WEI), Weibull type three (WEI3), Generalized Gamma (GG) and Generalized Pareto distributions (GP) and a nonparametric maximum likelihood estimation (NPMLE) of a mixing distribution for mixtures of Exponential distributions. The probability distribution functions of the WEI, GG, GP and the NPMLE for Exponential mixtures are given by Equations (1.1) to (1.5) of the supplementary material.

The results show that the the NPMLE and in particular the three component mixture-Gamma Lomax model fit much better than all other preliminary models except for the mixture-Gamma Lomax case, revealing that a mixture-based model to capture distributional multimodality is necessary and important. However, one drawback of the latter mixture model is the instability of the estimation of the implied tail index $\alpha_0:=\exp\{\nu_0\}$. Model fitting has been tested across various numbers of Gamma components $g$ and we examine how robust the estimates of tail heaviness across different $g$'s is. The results are shown in Table (ref). We see that the implied tail index $\alpha_0$ fluctuates greatly from smaller than $1.5$ to greater than $1.8$ across $g$, which does not make sense in practice because $g$ should control the body part of severity distribution only and bring very little impact on the estimated tail index. The main reason of seeing such an undesirable phenomenon is that the Lomax distribution, which is designed to capture the tail distribution, also calibrates to the body of the distribution and the MLE approach is found not very stable in estimating the tail parameter. This motivates the use of the composite model proposed in Section (ref), where the tail component only interacts with the body via the mixture probability. Also, note that while both AIC and BIC suggest a bigger number of components for the body (the optimal $g$ goes way beyond 15), this mainly reflects improvements of fitting small claims below 3,000. This should not be over-weighted because exact prediction of these small claims is less relevant in pricing, while excessive model complexity may impede interpretability. As a result, AIC and BIC may be less appropriate in determining the number of mixture components under this dataset.

Fitting results

In this section, we analyze the performance under the proposed mixture composite model with multi-type feature regularization for the covariates.

Distributional fitting

As in the preliminary analysis, we first fit the distribution of claim amounts under the proposed modeling framework, without considering covariates. Notice that in the mean excess plot of claim amounts (right panel of Figure (ref)) under the preliminary analysis, the plot becomes linear beyond claim severity of $270,000$ indicated by the vertical line of the plot. As a result, a reasonable choice of the splicing threshold is $\tau=270,000$. After fitting the proposed model across various choices of the number of body components $g$, we find that $g=5$ is the minimum number of components required to capture all the density nodes above a claim severity of $10,000$. The summary statistics of the fitted model is presented in Table (ref), the fitted versus empirical density plots are shown in Figure (ref), and the Q-Q and log-log plots of claim sizes are illustrated in Figure (ref). The model estimated tail index of 1.3817 roughly resembles that estimated by the asymptotic slope of the log-log plot which is 1.3 (left panel of Figure (ref)). Also, as expected we find that the model estimated tail index is robust across various choices of $g$. The density plots indicate that the fitted distribution captures all nodes representing a larger amount of claims, with multiple small nodes for smaller claims explained smoothly by one single component (to be precise, by the subgroup $j=1$ indicated by Figure (ref)). From the Q-Q plot, we see that the fitting performance is satisfactory except for very small claims ($y<100$) which are less relevant from an insurance pricing perspective. The fitted versus empirical log-log plot also indicates satisfactory fitting performance for the tail part. The fitted log-likelihood is $-719,309.1$, with $\text{AIC}=1,438,652$ and $\text{BIC}=1,438,640$, which is even slightly superior compared to the 5-Gamma Lomax distribution illustrated in Table (ref).

Effects of the covariates

We now include all variables described in Table (ref) and fit our proposed mixture composite regression model with LASSO and SCAD regularizations. Since all variables are included as categorical covariates, there is a total of $D=138$ parameters for each set of regressors. The grid searches are performed on $\lambda_1\in\{0.1n,0.2n,\ldots,6553.6n\}$, $\lambda_2\in\{0.1n_b,0.2n_b,\ldots,6553.6n_b\}$ and $\lambda_3\in\{0.1n_t,0.2n_t,\ldots,6553.6n_t\}$ to find optimal tuning parameters. The fitting performances for different model settings (without regression vs. with regression), penalty settings (without regularization vs. with regularization) and model selection criteria (pAIC, pBIC or CV with one standard deviation rule) are summarized in Table (ref). With a large number of covariates, we first note from the table that regularization of regression coefficients is a must, or else some parameters would diverge to very large values (due to overfitting), causing the algorithm to collapse eventually because of numerical instability (spurious solutions). As a result, for a full model as a benchmark for comparison, we need to apply a weak LASSO penalty, which sets very small $\lambda_l>0$ ($l=1,2,3$) such that no covariates are removed or merged. We next investigate the effect of the model selection criteria to the resulting fitted model. For both LASSO and SCAD penalties chosen as regularization function, pAIC results to very large models with a total of $\mathcal{N}=809$ parameters for LASSO and $\mathcal{N}=613$ for SCAD, indicating that many variables have predictive power on explaining all parts (body, tail and subgroup probabilities) of the claim severity distribution. The large number of predictors, however, makes the fitted models very difficult to interpret. Also, the selected model severities can vary greatly across various choices of initializations or grids for tuning parameters, because we find that model sizes within a range of about 150 to 1,000 parameters all have very similar AICs. In contrast, pBIC heavily penalizes the regression parameters and leads to a very small fitted model which chooses very few or even no variables useful to describe any parts of the distribution.

On the other hand, using CV with a one standard deviation rule provides fitted models with more reasonable complexity ($\mathcal{N}=112$ under LASSO or $\mathcal{N}=197$ under SCAD). Both LASSO and SCAD penalties suggest that there are not any systematic effects in the tails that are explained by the available variables. On the other hand, both penalty functions reveal similar sets of variables important to explain the body and subgroup probability parts. The higher model complexity under SCAD is mainly due to more granular mergers among different levels of some variables (such as driver's age). Under LASSO, the resulting AIC under the CV approach is close to that under the corresponding pAIC approach, while the BIC is just slightly inferior to the pBIC approach.

Table (ref) also shows the performance of the LASSO and SCAD CV-selected models after the model refit procedure. Recall from Section (ref) that the refitting procedure involves re-estimation of parameters for the shrinked model with regularization terms excluded to reduce biasedness in the estimated parameters. For the LASSO penalty, the improvements of the log-likelihood, AIC and BIC are all expected after refitting. For the SCAD penalty, since the concavity of SCAD penalty function already mitigates the biasedness of estimated parameters (fan2001variable), there is no apparent improvement of the fitting performance after performing the refitting procedure. After refitting, the LASSO penalty approach results to superior fitting performance compared to the SCAD approach, as evidenced by lower AICs and BICs. As a result for conciseness concern, we focus solely on the CV approach with LASSO penalty as model selection criterion in the following analysis.

The final refitted model suggests that the subgroup probabilities $\pi_j(\bm{x};\bm{\alpha})$, $1\le j \le g+1$, are influenced by the variables as follows.

These results are also presented by plots in Figure (ref), which display the probability being classified to tail component versus various variables. The points indicated as triangle ($\triangle$) and square ($\square$) correspond to the fitted and empirical probabilities, respectively. Details on visualizing covariate influences through non-parametric approaches are discussed by FUNG2019MoEApplication. As we can see from the figure, conditioned on any categories/ levels of any explanatory variables, the fitted and empirical probabilities match very well, reflecting the ability of the proposed regression model to capture well the covariates influence. The green dotted line is the overall empirical tail probability across all observations. The blue and red intervals are respectively the 95% Wald-type and Efron bootstrap CIs presented in Section (ref). The CIs generated by the two approaches reconcile well.

The model chooses a smaller set of variables which are important in explaining the body distributions $f(y_i;\exp\{\bm{\beta}_j^T\bm{x}_i,\phi_j)\})$, reflecting more heterogeneity among subgroup probabilities than within-subgroup average claim sizes:

Finally, the fitted model suggests that none of the explanatory variables are significantly influential to the tail distribution $h(y_i;\theta,\exp\{\bm{\nu}^T\bm{x}_i\})$. Overall, the effects on various variables to the average claim severity are demonstrated in Figure (ref).

Summary findings

In this real data analysis, we get a deeper understanding on the influence of policyholder attributes to the claim severity distribution with highly complex structure including multimodality and tail-heaviness. Using the proposed mixture composite modeling framework embedded with a variable selection approach, we find that the explanatory variables most prominently impact the subgroup probabilities of the severity distribution, explaining the unobserved heterogeneity of policyholder risk profiles and/or claim types. Fewer variables explain well the body part of the distribution, reflecting relatively homogeneous claim severity distributions conditioned on the subgroups where each claim is belonging to. This finding is in contrast to many traditional regression models widely adopted in actuarial practice, including GLM and GAM, where regression links are set to capture the systematic effects in distributions instead of the subgroup heterogeneity.

Further, we do not find any variables significantly influencing the tail-heaviness of the claim severity distribution, which may be the result of scarcity of large claims (only around 2,400 claims exceed the splicing threshold $\tau$) to allow for statistically significant covariates influence to the tail part. This empirically verifies the legitimacy of actuarial practice where covariates influence is often excluded in modeling large claims. In actuarial literature, we refer to laudage2019severity who also refrains from incorporating regression in the tail part of their severity distribution.

Discussions

In this article, we considered a mixture composite regression model for addressing several challenges when modeling claim severities such as multimodality and tail-heaviness of claims, extending the framework of reynkens2017modelling who considered the case without covariates. For variables selection, we proposed a group-fused regularization approach. Our covariates may influence the mixture probabilities, the body and the tail of the claim size distribution, in such way that model interpretability is preserved. This approach enables regularization under multi-type variable settings. For this setup, we developed an asymptotic estimation theory which justified the efficiency of the proposed method. In particular, we showed that the method we presented is: (i) consistent in terms of covariate selection since, when the sample size goes to infinity, it will merge and shrink correctly regression coefficients across all modeling parts, and (ii) the parameters of the reduced model are asymptotically normal. The implementation was illustrated by a real data application which involved fitting claim size data from a Greek automobile insurance company. Maximum likelihood estimation of the model parameters was achieved through a novel Generalized Expectation-Maximization algorithm that was demonstrated to perform well.

Furthermore, it is worth noting that instead of following a data driven approach for selecting the number of mixture components in the body area based on specification criteria, as is done herein, an interesting direction of further research would be to extend the framework to a non-parametric maximum likelihood estimation approach which can be utilized for automated selection of the number of mixture components.

Finally, it is worth noting that while the proposed composite model mitigates instabilities of tail index estimations inherited by finite mixture models, selection of the splicing threshold is often subjective. Therefore, it would be worth to explore alternative approaches for robust estimation of the tail index. One possible way is to modify the maximum likelihood approach for parameter estimation such that an observation with a larger claim severity has a higher relative importance in determining the model parameters. Another possible way is to explore models which bridge the gap between finite mixture models and composite models, and share the advantages of both model classes.