Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
75,019 characters · 16 sections · 82 citation commands
A Posteriori Risk Classification and Ratemaking with Random Effects in the Mixture-of-Experts Model
A well-designed framework for risk classification and ratemaking in automobile insurance is key to insurers' profitability and risk management, while also ensuring that policyholders are charged a fair premium according to their risk profile. For a new policyholder, risk classification and ratemaking are usually done on an a priori basis, whereby the insurer only knows a set of the policyholder's covariates such as age, gender, vehicle specifications, etc. As time goes by, the insurer gains additional, up-to-date insights into the policyholder's risk profile from their claim history, including frequency and severity, which leads to a posteriori risk classification and ratemaking.
The use of claim history for a posteriori risk classification and ratemaking is a classical problem which has been studied in depth in the actuarial literature. Early works in credibility theory, such as buhlmann1967experience, norberg1979credibility and buhlmann2005course, assume some common parameters underlying the distribution of insurance losses. One uses the observed claim history to infer the posterior distribution of the parameters, which then yields the posterior distribution of future losses given the history. From a practical perspective, the Bonus-Malus System (BMS) is perhaps one of the most widely used approaches, see e.g. lemaire1995bonus and denuit2007actuarial. Based on the claim history (typically the number of claims in the year prior to policy renewal), policyholders are (re-)classified into one of a number of pre-specified risk classes according to certain transition rules, whereby each risk class corresponds to a premium relativity which reflects the level of risk. However, in their classical formulation, neither credibility theory nor BMS considers covariate information, which is usually deemed as important indicators of policyholders' risk characteristics. To this end, there has been an abundance of literature that aims to apply more sophisticated statistical models, which typically involve a regression component, to the problem of a posteriori risk classification and ratemaking. Most notably, random effects have been a popular choice for modelling the temporal dependence between past and future claim behaviour. For example, many authors have considered adding random effects in Generalized Linear Models (GLM), which results in Generalized Linear Mixed Models (GLMM), see e.g. dionne1989generalization, dionne1992automobile, pinquet1998designing, frangos2001design and boucher2006fixed, whereby the posterior distribution of random effects given claim history is used for prediction. Another important consideration is the dependence structure between multiple coverages which is common in automobile insurance, see e.g. pinquet1998designing, gomez2008univariate, boucher2009number, gomez2016bivariate and tzougas2021multivariate for using shared random effects to model such dependence. Besides, while some works mainly focus on claim frequency alone, many researchers have also attempted to incorporate claim severity and its dependence structure with frequency, for example, ni2014bonus, park2018does, oh2020bonus and oh2021designing. Furthermore, to overcome certain restrictive assumptions in GLM, finite mixture models have recently become popular in a posteriori risk classification and ratemaking for more flexible and accurate modelling of claim frequency and severity, as used in bermudez2012finite, tzougas2014optimal, tzougas2018bonus and tzougas2021multivariate.
In this paper, we propose to apply a flexible regression model, called the Mixed LRMoE, to the problem of a posteriori risk classification and ratemaking. Compared with existing approaches to this problem, our proposed method enjoys several distinct advantages, such as an intuitive and interpretable model structure (see (ref)), the flexibility to model any mixed effects model (see (ref)), and superior performance in goodness-of-fit and adequacy in a posteriori risk classification and ratemaking compared with benchmark models (see (ref)). The Mixed LRMoE as a general modelling framework has recently been introduced in Fung2022MixedLRMoETheory as an extension to the Logit-weighted Reduced Mixture-of-Experts (LRMoE) model. The latter was first developed in fung2019class and has subsequently been applied to various insurance modelling problems such as correlated claim frequencies and reporting delay (see (ref) for an overview). In order to adapt to the problem of a posteriori risk classification and ratemaking, we propose to add policyholder-level random effects in a multiyear portfolio which results in the Mixed LRMoE. Similar to many papers cited above, the addition of random effects introduces dependence between observations across multiple policy years of the same policyholder, from which the posterior distribution of random effects is inferred and then utilized for a posteriori risk classification and ratemaking. Our work also intersects with mixture model-based approaches such as tzougas2021multivariate, in that the Mixed LRMoE model allows for more flexible and accurate modelling of the loss distribution compared with classical regression models such as GLM. In the broader class of general mixture-of-experts (MoE) models, our work is closely related to yau2003finite, ng2007extension and ng2014mixture, where random effects are also incorporated to account for heterogeneity observed in real data. However, the Mixed LRMoE presented in this paper has an arguably simpler model structure. A detailed discussion on various properties of the Mixed LRMoE and a brief comparison between our work and existing literature are provided in (ref)
From a modelling perspective, Fung2022MixedLRMoETheory shows that the Mixed LRMoE is dense in the space of any mixed effects models subjected to mild regularity conditions. It means that the Mixed LRMoE is flexible enough to resemble any complex characteristics inherited from any mixed effects models, including the joint distribution, the regression pattern, the random intercept, and the random slope, to an arbitrary degree of accuracy. This theoretical result is an extension of fung2019class, whereby the LRMoE is shown to be dense in the space of regression models, justifying the versatility and parsimony of the LRMoE with a reduced model structure. The addition of random effects is also crucial for modelling the temporal dependence between observations across different policy years in a large, multiyear insurance portfolio. These desirable features of Mixed LRMoE make it a powerful tool for a posteriori risk classification and ratemaking, as demonstrated by our real data analysis in (ref), where we apply our proposed framework to an automobile insurance dataset. Our model is shown to outperform classical models in terms of goodness-of-fit to data, while offering fair and interpretable risk classification and ratemaking which accurately reflect policyholders’ claim history.
Besides methodologically applying the Mixed LRMoE model to a posteriori risk classification and ratemaking, our second major contribution is the development of a variational inference (VI) algorithm for parameter estimation and posterior inference of random effects. In general, when random effects are included in regression models, parameter estimation and inference may be challenging due to typically intractable likelihood functions. As a classical approach, one may consider applying the Best Linear Unbiased Predictor (BLUP) procedure for obtaining the realization of random effects, combined with Restricted/Residual Maximum Likelihood (REML) for estimating the model parameters, see e.g. henderson1973sire, henderson1975best, mclean1991unified for Linear Mixed Models, mcgilchrist1994estimation and mcgilchrist1995derivation for Generalized Linear Mixed Models, and yau2003finite and ng2007extension for MoE models. Alternatively, one may choose to estimate the parameters from the marginal likelihood by numerically integrating out the random effects using e.g. the Gauss-Hermite Quadrature (pinheiro1995approximations) or the Laplace approximation (e.g. breslow1993approximate and raudenbush2000maximum). One may also apply Markov Chain Monte Carlo (MCMC) methods (e.g. zeger1991generalized, booth1999maximizing and brooks2011handbook) for generating samples of random effects from their posterior distribution given the observed data, based on which the posterior of model parameters can also be obtained. A comparison of these methods for models with random effects can be found in browne2006comparison. However, the aforementioned methods may not be suitable for the application of a posteriori risk classification and ratemaking. For example, when working with large insurance portfolios, it desirable to develop an algorithm which scales with the number of random effects and the size of datasets, which may be difficult for numerical integration or MCMC methods. Also, it is desirable to obtain posteriori distributions, rather than point estimates, of certain quantities of interest (e.g. a posteriori premium based on different premium principles), which are not produced by either BLUP or numerical integration methods. Hence, in place of these classical methods, we opt to use VI primarily for its superior speed and scalability for large insurance portfolios. Besides estimating model parameters with computational efficiency, our VI algorithm also directly produces the approximated posterior distribution of random effects for each individual policyholder, which is key for a posteriori risk classification and ratemaking for future policy years. Further, while VI methods have been widely used in the machine learning community as an alternative to computationally more expensive methods such as MCMC (blei2017variational), there has been little application of VI in the actuarial literature (see e.g. kuo2020individual and gomes2021insurance). We hope our paper serves as another example to showcase the potentials of VI methods for analyzing the ever-growing amount of data available for insurance applications.
The remainder of this paper is organized as follows. (ref) reviews the LRMoE model and introduces the Mixed LRMoE. Then, (ref) develops a stochastic variational Expectation-Conditional-Maximization (ECM) algorithm for estimating model parameters and inferring the posterior distribution of random effects. Next, (ref) presents two simulation studies which aim to numerically illustrate and examine the proposed estimation algorithm, and (ref) contains an application of our proposed framework on a real insurance dataset. Finally, (ref) concludes with a brief discussion and outlook for future research directions.
In this section, we first give an overview of the LRMoE modelling framework, including model formulation, theoretical properties, implementation and application in actuarial contexts. Then, we extend the LRMoE model with random effects to account for the temporal dependence across different policy years. Finally, we provide some discussion on the Mixed LRMoE and a brief comparison between our work and existing literature.
The LRMoE model first introduced in fung2019class is formulated as follows. Let $\bm{x}_i$ denote a $P$-dimensional vector of covariates of policyholder $i$ such as demographic information and vehicle specification. Given $\bm{x}_i$, the policyholder is classified into one of $g$ latent risk classes by the logit gating function
where $\bm{\alpha}_{j}$ is a vector of regression coefficients for latent class $j$. Within each latent class $j$, a $D$-dimensional vector of response variable(s) $\bm{y}_i$ such as claim frequency and severity is modelled by an expert function $f_{j}(\bm{y}_i; \bm{\psi}_j)$, where $\bm{\psi}_j$ denotes the parameters of the expert function. Consequently, the likelihood function for a portfolio of $n$ policyholders is given by
where $\bm{\alpha} = (\bm{\alpha}_{1}^T, \bm{\alpha}_{2}^T, \dots, \bm{\alpha}_{g}^T)^T$ and $\bm{\Psi} = \{ \bm{\psi}_1, \bm{\psi}_2, \dots, \bm{\psi}_g \}$ are the model parameters to estimate given the observed data $(\bm{X}, \bm{Y}) = \{ (\bm{x}_i, \bm{y}_i): i=1, 2, \dots, n \}$. We assume conditional independence among all dimensions in $\bm{y_i}$ given the latent class $j$ such that $f_{j}(\bm{y}_i; \bm{\psi}_j) = \prod_{d=1}^{D}f_{jd}(y_{id}; \bm{\psi}_{jd})$ for $d=1, 2, \dots, D$, where $y_{id}$ is the $d$-th dimension in $\bm{y}_i$ and $f_{jd}$ is the expert function for $y_{id}$ with parameters $\bm{\psi}_{jd}$.
The LRMoE model can be viewed as a simplification of the general MoE model (see e.g. jordan1994hierarchical), whereby the gating function is restricted to multiple logistic functions and the regression on covariates in the expert functions is eliminated. It is shown in fung2019class that such simplification will not reduce modelling flexibility, provided the expert functions satisfy some mild conditions. In other words, the LRMoE model is capable of achieving the same level of goodness-of-fit as the general mixture-of experts with a much simpler model structure. In the meantime, the simplified model structure of LRMoE provides the following intuitive model interpretation in insurance contexts. Based on covariates $\bm{x}_i$ which are indicative of individual risk profiles, policyholders are classified into latent risk groups by a commonly used function for classification problems. Within the same latent group $j$, the individual risk profiles are naturally assumed to be homogeneous by sharing the same expert function $f_{j}(\bm{y}_i; \bm{\psi}_j)$ whose parameters are independent of policyholder information.
Thanks to its flexibility and interpretability, the LRMoE model has been applied to many actuarial modelling problems. fung2019classapplication used it for modelling correlated claim frequencies of two types of automobile insurance coverage, where the LRMoE mixture of Erlang Count experts is shown to outperform the negative binomial GLM (with and without zero inflation). fung2020fittingcensor discussed fitting LRMoE to censored and truncated data which are commonly encountered when modelling claim severity or reporting delays. The extended model is applied to insurance pricing with policy deductibles and prediction of incurred but not reported (IBNR) claims. In fung2021mixture, the LRMoE is further extended to include composite or slicing expert functions which account for multi-modal and heavy-tailed distributions. For implementation of LRMoE, software packages written in R (tseung2020lrmoeR) and in Julia (tseung2021lrmoejl) are readily available for use, which offer a wide selection of expert functions commonly used for actuarial modelling and utility functions for predictive analysis and model visualization.
As with many mixture models, parameter estimation for LRMoE is done using the Expectation-Conditional-Maximization (ECM) algorithm (see e.g. DEMPSTER1977EM and Mclachlan2004Finite). Details of the ECM algorithm for LRMoE can be found in the papers cited above. For Mixed LRMoE, we combine the same ECM algorithm with VI methods in order to deal with intractable marginal likelihood due to the presence of random effects, which will be presented in (ref).
\afterpage{
}
In the context of a posteriori risk classification and ratemaking, it is important to utilize information about policyholders' claim history to make predictions for the upcoming policy years. In effect, one takes advantage of the dependence structure in the claim history across different policy years generated by the same policyholder. Note that such dependence structure has not been accounted for by the LRMoE model, due to the assumption of independence between observations $(\bm{x}_i, \bm{y}_i)$ as indicated by the likelihood function in (ref). To incorporate dependence between observations across different policy years, we propose to add policyholder-level random effects to the LRMoE model, which results in the Mixed LRMoE model. In this subsection, we first formulate the Mixed LRMoE in a general setting following Fung2022MixedLRMoETheory, and then discuss the special case with only policyholder-level random effects.
Assume each observation $(\bm{x}_i, \bm{y}_i)$ is equipped with a vector of random effects $\bm{w}_i = (w_{i1}, w_{i2}, \dots, w_{iL})$, where $L$ is the total number of levels of different random effects. For the $l$-th level of random effect, $l=1, 2, \dots, L$, we assume there are in total $S_l$ factors $\{ w_{l}^{(s)}\}_{s=1, 2, \dots, S_l}$, and each observation $i$ is mapped into one of these factors by a known function $c_l(\cdot)$ such that $w_{il} = w_{i'l} = w_{l}^{(s)}$ if $c_l(i) = c_l(i') = s$ for $s = 1, 2, \dots, S_l$. Equivalently, the mapping function $c_l(i)$ can be represented by a $S_l$-vector $\bm{t}_{il}$ where exactly the $c_l(i)$-th element is one and the others are zero (see also (ref) for an example).
Let $\bm{w} = \{w_{l}^{(s)}\}_{l=1, 2, \dots, L; s=1, 2, \dots, S_l}$ denote the collection of random effects across all levels and all factors, which are assumed to be independent across $l$ and $s$. We also assume their distribution and density functions are pre-specified by $\bm{\Phi}(\cdot)$ and $\bm{\phi}(\cdot)$ with no extra parameters such that
where $\Phi_l(\cdot)$ and $\phi_l(\cdot)$ are, respectively, the distribution and density functions for the $l$-th level of random effects $\{ w_{l}^{(s)}\}_{s=1, 2, \dots, S_l}$ for $l=1, 2, \dots, L$. In general, one may specify a priori any distribution for $\bm{\Phi}(\cdot)$, but a common choice for random effects is the normal distribution. In this paper, we will set each $\Phi_l(\cdot)$ to be a standard normal distribution for $l=1, 2, \dots, L$. More discussions on the choice of $\bm{\Phi}(\cdot)$ are given in (ref).
Similar to the covariates $\bm{x}_i$, we assume the random effects $\bm{w}_i$ influences only the gating function. In addition, we assume there are coefficients $\bm{\beta}_j$, $j=1, 2, \dots, g$, multiplied to the random effects, which serve as scaling factors that also affect the gating functions and add to the modelling flexibility by compensating the lack of parameters in $\bm{\Phi}(\cdot)$. Consequently, the gating function in a Mixed LRMoE model is given by
Unlike the gating functions, the expert functions are assumed to be independent of both the covariates $\bm{x}_i$ and the random effects $\bm{w}_i$, as illustrated in (ref). Note this is the same assumption used in the LRMoE model without random effects. Consequently, given the realization of random effects $\bm{w}$, the likelihood function of Mixed LRMoE is
while the likelihood with random effects integrated out is given by
where $d\bm{w} = \prod_{l=1}^{L}\prod_{s=1}^{S_l} d w_l^{(s)}$ and the subscript of the expectation operator $\E$ indicates the expectation is calculated by integrating out $\bm{w}$ with respect to $\bm{\phi}(\cdot)$.
The most important property of the Mixed LRMoE is the denseness property, which justifies the flexibility of the proposed model in capturing a broad range of complex multilevel data characteristics. While the theoretical result has been rigorously developed by Fung2022MixedLRMoETheory, we hereby briefly describe and interpret the result without extensive mathematical treatments.
Let $F(\bm{Y};\bm{\alpha},\bm{\beta},\bm{\Psi}|\bm{X})$ be the joint distribution function of $\bm{Y}$ given $\bm{X}$ under the proposed Mixed LRMoE model, which is given by
where $F_{j}(\bm{y}_i; \bm{\psi}_j)$ is the distribution function of $f_{j}(\bm{y}_i; \bm{\psi}_j)$. Also, denote $H(\bm{Y}|\bm{X})$ as the joint distribution of $\bm{Y}$ given $\bm{X}$ under an arbitrary mixed effects model. Under some mild regularity conditions, Fung2022MixedLRMoETheory proves that for any target mixed effects model $H(\bm{Y}|\bm{X})$, there exists a sequence of model parameters $\{(\bm{\alpha}^{[s]},\bm{\beta}^{[s]},\bm{\Psi}^{[s]})\}_{s=1,2,\ldots}$ (note that the number of latent risk classes $g$ may increase as $s$ increases) such that $F(\bm{Y};\bm{\alpha}^{[s]},\bm{\beta}^{[s]},\bm{\Psi}^{[s]}|\bm{X})$ converges in distribution to $H(\bm{Y}|\bm{X})$ uniformly on $\bm{X}$ as $s\rightarrow\infty$. Note that the target mixed effects model $H(\bm{Y}|\bm{X})$ may carry very complicated model characteristics, including but not limited to the joint loss distribution (e.g., distributional multimodality and dependence across business lines), the regression link (e.g., non-linear or interactive influence of policyholder attributes to the losses), the random intercept (e.g., latent impacts to each policyholder), and the random slope (e.g., random effects interact with policyholder attributes). As a result, the denseness theorem justifies the versatility of the proposed Mixed LRMoE in simultaneously capturing all these features to an arbitrary degree of accuracy. Moreover, the denseness theorem only requires that $\bm{\Phi}(\cdot)$ is continuous. Hence, one has the freedom to choose any continuous distributions for the random effects without impeding the flexibility of the Mixed LRMoE. Motivated by the computational convenience (see (ref) below), we select $\Phi_l(\cdot)$ ((ref)) to be a standard normal distribution, such that $\bm{\Phi}(\cdot)$ follows a multivariate standard normal distribution.
\afterpage{
}
Before proceeding to parameter estimation, we make the following remarks on the model formulation of Mixed LRMoE and provide a brief comparison with existing literature.
Firstly, in (ref) we have given a general formulation of Mixed LRMoE with potentially multiple levels of random effects when $L>1$. For the application in a posteriori risk classification and ratemaking in this paper, we set $L=1$ to add only policyholder-level random effects. In this case, $n$ is the total number of policy year observations out of $N_0$ unique policyholders, such that each factor in $\{w_{1}^{(s)}\}_{s=1, 2, \dots, N_0}$ represents the individual risk of one unique policyholder. An illustration for one such policyholder is shown in (ref). Other than a posteriori risk classification and ratemaking, one may consider applying the Mixed LRMoE to other modelling problems with multiple levels of latent risks, such as modelling geographical risks with a nested structure for random effects with $L=2$ levels, where $l=1$ represents city-level random effects and $l=2$ represents the latent risks for specific neighbourhoods. For illustration purposes, we will leave the application of Mixed LRMoE with $L>1$ for future investigation, and only demonstrate a simulation study for $L=2$ in (ref).
Secondly, similar to many previous works such as those cited in (ref), our paper also utilizes random effects for modelling temporal dependence among different policy years of the same policyholder, but we have done so in a slightly different fashion. Many previous papers have proposed mixed models whereby the certain model parameters are shared across different observations. For example, one may assume the claim frequency $N_{it}$ of policyholder $i$ in the $t$-th year follows $\textrm{Poisson}(\theta_{it})$, and then uses the observed data $\{N_{it}: t=1, 2, \dots\}$ to infer the posterior of the intensity parameter. In contrast, our formulation of the Mixed LRMoE treats the random effects $\bm{w}$ in a similar way as the fixed effects $\bm{x}_i$, which essentially serve as a regressor in the gating function. Rather than imposing certain changing dynamics on model parameters, the formulation of Mixed LRMoE actually resembles, to a large extent, classical approaches of longitudinal data modelling with random effects, see e.g. diggle2002analysis and fitzmaurice2012applied.
Finally, the Mixed LRMoE model shares varying degrees of similarity with previous works which attempt to incorporate random effects in the general MoE framework. For example, yau2003finite proposes a two-component MoE with random effects in both the logit gating function and normal experts. ng2007extension considers a similar framework but uses Bernoulli experts for a classification problem, while ng2014mixture adds random effects only to the expert functions. In contrast, our present work focuses on a specific subclass of Mixed MoE model where random effects only influence the latent class probabilities through the gating function, while the expert functions are kept independent of covariates and random effects. Besides possessing the same level of modelling flexibility due to denseness (as discussed in (ref)), this simplified model structure leads to an easier implementation of parameter estimation. As will be evident in (ref), since the estimation procedures of gating and expert functions can be separated to some extent, the Mixed LRMoE model actually allows for more flexible choices and combinations of expert functions which are customized to different modelling problems (see also (ref)). By restricting the random effects to only the gating functions, we are able to develop a unified estimation algorithm which caters for different choices and combinations of expert functions.
In this section, we develop a stochastic variational ECM algorithm for estimating model parameters and for inferring the posterior distribution of random effects for Mixed LRMoE. We first present an overview of variational inference methods in general, and then provide details of the implementation for Mixed LRMoE with one single type of random effect. Discussion on model identifiability, model selection and generalization of this algorithm is given at the end of this section.
In this subsection, we first provide an overview and motivation of variational inference methods. We start with the exact posterior distribution of random effects $\bm{w}$
which may be complicated due to the dependence on both the model parameters $(\bm{\alpha}, \bm{\beta}, \bm{\Psi})$ and the observed data $(\bm{X}, \bm{Y})$. To circumvent this numerical challenge, we assume the exact posterior can be reasonably approximated by a variational distribution $\bm{q}(\bm{w} ; \bm{\Theta})$ where $\bm{\Theta}$ is the variational parameters, which are assumed to be independent of the model parameters and observed data. This produces a numerically more tractable lower bound of the marginal likelihood in (ref), also known as the Evidence Lower Bound (ELBO) in the variational inference literature. More specifically, by taking logarithm of (ref), utilizing the variational distribution, and applying Jensen's inequality, we obtain the following ELBO of the marginal loglikelihood.
where $\textrm{KL}\left[ \bm{q}(\bm{w} ; \bm{\Theta}) || \bm{\phi}(\bm{w}) \right]$ is the Kullback-Leibler (KL) divergence between the variational posterior $\bm{q}(\bm{w} ; \bm{\Theta})$ and the prior $\bm{\phi}(\bm{w})$ of random effects.
Instead of directly maximizing the marginal likelihood in (ref), we aim to maximize the ELBO $\underline{\ell}(\bm{\alpha}, \bm{\beta}, \bm{\Psi}, \bm{\Theta}; \bm{X}, \bm{Y})$ in (ref), hoping that the optimal parameters which maximize this lower bound are close to the true optimal parameters which maximize the actual loglikelihood. The main advantage is the tractability of the approximate posterior of random effects $\bm{w}$, which is essentially specified by parameters $\bm{\Theta}$ independent of all the other model parameters and observed data. As will be evident in the next subsection, sampling from the approximated posterior is easier and faster than MCMC methods, since the latter works with a more complex exact posterior and typically requires a burn-in period. This may offer significant numerical efficiency, especially in high-dimensional cases where there are many types of random effects and each type of random effect has many levels. Meanwhile, the obvious trade-off is obtaining only the approximated solutions to the estimated model parameters and the approximated posterior distributions of the random effects. While the goodness of approximation and convergence properties for variational inference remain an open problem (see e.g. blei2017variational), our numerical simulations in (ref) and real data analysis in (ref) show promising results. This may serve as an empirical evidence for applying variational inference methods to insurance problems where an approximated solution may be acceptable in the presence of large datasets.
For variational inference, one needs to specify a family of parametric distributions for the approximated posterior $\bm{q}(\bm{w} ; \bm{\Theta})$. In this paper, we follow standard practices and use the mean-field variational family, whereby the posterior of latent variables, i.e. random effects $\bm{w}$, is a factorized multivariate normal distribution. More specifically, we assume the posterior of $w_{l}^{(s)}$ is a normal distribution with mean $\mu_{l}^{(s)}$ and standard deviation $\sigma_{l}^{(s)}$ for $s=1, 2, \dots, S_l$ and $l=1, 2, \dots, L$, which are independent across all levels $l$ and all factors $s$. Mathematically,
For notational convenience, we write $\bm{\Theta} = \{(\bm{\mu}_l, \bm{\Sigma}_l)\}_{l=1, 2, \dots, L}$, where $\bm{\mu}_l = (\mu_{l}^{(1)}, \mu_{l}^{(2)}, \dots, \mu_{l}^{(S_l)})^T$ is the posterior mean vector and $\bm{\Sigma}_l = \textrm{diag}((\sigma_{l}^{(1)})^2, (\sigma_{l}^{(2)})^2, \dots, (\sigma_{l}^{(S_l)})^2)$ the diagonal covariance matrix for the $l$-th level of random effect.
When $L=1$, given the factorization of likelihood across $s=1, 2, \dots, S_1$, different factors of the same level of random effect are in fact independent, both in the prior and the posterior distribution. Hence, in our application of the Mixed LRMoE with only policyholder-level random effects, the only source of error of variational inference is the approximation of the exact posterior by a normal distribution. However, when there are multiple types of random effects (e.g. the multilevel example in (ref)), especially in the case of certain dependence structures (e.g. multiple crossed random effects), the independence assumption in the mean-field variational family may create an additional source of error of approximation.
With the approach of variational inference and the choice of the mean-field variational family $\bm{q}(\bm{w} ; \bm{\Theta})$, we now develop a stochastic variational ECM algorithm for estimating the model parameters $(\bm{\alpha}, \bm{\beta}, \bm{\Psi})$, as well as inferring the posterior of random effects $\bm{w}$ represented by the variational parameters $\bm{\Theta} = \{(\bm{\mu}_l, \bm{\Sigma}_l)\}_{l=1, 2, \dots, L}$.
On a high level, our estimation algorithm proceeds in an iterative manner which seeks to conditionally maximize the ELBO in (ref) with respect to one set of parameters while keeping others fixed. Consequently, the algorithm will ultimately arrive at a local optimum for the ELBO of the marginal loglikelihood. First, we initialize the model parameters $({\bm{\alpha}}, {\bm{\beta}}, {\bm{\Psi}})$ using the clusterized method of moments (CMM), similar to e.g. gui2018fitting. Meanwhile, the variational parameters $\bm{\Theta}$ can be initialized such that $\bm{\mu}_l = \bm{0}$ and $\bm{\Sigma}_l = \bm{I}$ for $l=1, 2, \dots, L$ (i.e. assuming a multivariate standard normal distribution), which is consistent with standard practices in the VI literature. Then, our algorithm iterates through the following steps until convergence.
E-Step: At iteration $t+1$, given the current model parameters $(\bm{\alpha}^{(t)}, \bm{\beta}^{(t)}, \bm{\Psi}^{(t)})$ and variational parameters $\bm{\Theta}^{(t)} = \{(\bm{\mu}_l^{(t)}, \bm{\Sigma}_l^{(t)})\}_{l=1, 2, \dots, L}$, we calculate the expectation of the complete-data ELBO, which results in the objective function $Q^{(t+1)}(\bm{\alpha}, \bm{\beta}, \bm{\Psi}, \bm{\Theta}; \bm{X}, \bm{Y})$.
CM-Steps:
Next, we describe each of these steps in more detail. In the E-Step, we augment the usual latent variables $\bm{Z} = \{Z_{ij}, i=1, 2, \dots, n$ and $j=1, 2, \dots, g \}$ such that $Z_{ij}=1$ indicates $\bm{y}_i$ is generated by the $j$-th latent class and $Z_{ij}=0$ otherwise. Consequently, the complete-data ELBO is given by
We then calculate the expected value of $\bm{Z}$ given the current values of model and variational parameters, which yields the following objective function $Q^{(t+1)}(\bm{\alpha}, \bm{\beta}, \bm{\Psi}, \bm{\Theta}; \bm{X}, \bm{Y})$ to be maximized in the CM-Steps.
where
and the change of order of integration is justified by $\sum_{i=1}^{n} \sum_{j=1}^{g} Z_{ij} \log \left[ \pi_{j}(\bm{x}_i, \bm{w}_i; \bm{\alpha}, \bm{\beta}) \right] \leq 0$. Given the realization of random effects $\bm{w}$, the conditional expectation on the right-hand-side of (ref) is evaluated as
Note that the unconditional expectation of $Z_{ij}$ by integrating out $\bm{w}$ admits no closed-form solution. However, the normality assumption on the posterior of $\bm{w}$ allows for the following numerical evaluation through Monte Carlo simulation which entails little computational burden.
where $\bm{w}^{[m]}$ denotes the $m$-th sample of random effects generated from the variational distribution $\bm{q}(\cdot;\bm{\Theta}^{(t)})$.
Next, in CM-Step (i), given the current variational parameters $\bm{\Theta}^{(t)}$, the maximization of $Q^{(t+1)}(\bm{\alpha}, \bm{\beta}, \bm{\Psi}, \bm{\Theta}; \bm{X}, \bm{Y})$ is divided into subproblems in $Q_1^{(t+1)}(\bm{\alpha}, \bm{\beta}; \bm{X}, \bm{\Theta}^{(t)})$ and $Q_2^{(t+1)}(\bm{\Psi}; \bm{Y}, \bm{\Theta}^{(t)})$ such that
and
Given the realization of random effects $\bm{w}$, the right-hand-side of (ref) without the expectation operator can be maximized using the iteratively re-weighted least squares (IRLS) method (see e.g. jordan1994hierarchical and fung2019classapplication). To account for the randomness in $\bm{w}$, we adapt the deterministic IRLS procedure to its stochastic version which is described in detailed in (ref). In the meantime, given $z_{ij}^{(t)}$ obtained from the E-Step, the maximization of $Q_2^{(t+1)}(\bm{\Psi}; \bm{Y}, \bm{\Theta}^{(t)})$ over the expert parameters $\bm{\Psi}$ proceeds exactly the same as described in e.g. fung2019classapplication and tseung2021lrmoejl, which is independent of the variational parameters $\bm{\Theta}^{(t)}$. Details are omitted here and we refer interested readers to the cited papers.
Finally, in CM-Step (ii), the complete-data ELBO is maximized over the variational parameters $\bm{\Theta} = \{(\bm{\mu}_l, \bm{\Sigma}_l)\}_{l=1, 2, \dots, L}$ given the updated model parameters $(\bm{\alpha}^{(t+1)}, \bm{\beta}^{(t+1)}, \bm{\Psi}^{(t+1)})$. The variational parameters affect the objective function in (ref) through both the expectation operator $\E_{\bm{w} \sim \bm{q}(\cdot ;\bm{\Theta})} \left[ \cdot \right]$ and the KL divergence term $\textrm{KL}\left[ \bm{q}(\bm{w} ; \bm{\Theta}) || \bm{\phi}(\bm{w}) \right]$. Assuming the mean-field variational family $\bm{q}(\cdot;\bm{\Theta})$, the optimization over $\bm{\Theta}$ can be effectively done by a standard reparameterization technique on the random effects $\bm{w}$ combined with a simple gradient descent. For brevity, details are deferred to (ref).
In addition to the estimated model parameters $(\hat{\bm{\alpha}}, \hat{\bm{\beta}}, \hat{\bm{\Psi}})$, our algorithm also yields the variational parameters $\hat{\bm{\Theta}} = \{(\hat{\bm{\mu}}_l, \hat{\bm{\Sigma}}_l)\}_{l=1, 2, \dots, L}$ which completely specify the approximated posterior distribution of random effects $\bm{w}$. Despite no closed-form formulas for various quantities of interests such as the posterior mean of response $\bm{y}_i$ (see also (ref)), their approximated values can be efficiently calculated by sampling from the variational posterior distribution which is assumed to be multivariate normal.
As with many mixture models, certain restrictions are imposed for the Mixed LRMoE to be identifiable when conducting parameter estimation. In order to avoid label-switching between latent components (see e.g. jiang1999identifiability and fung2019classapplication), we fix $\bm{\alpha}_g = \bm{0}$ and $\bm{\beta}_g = \bm{0}$ as vectors of zeros, so the last latent class serves as a reference class. In addition, we fix $\bm{\beta}_1 = \bm{1}$ as a vector of ones to avoid arbitrary scaling of magnitude and switching of positive and negative signs of the random effects $\bm{w}$. Consequently, we need to estimate the coefficients $\bm{\beta}$ multiplied to the random effects only when there are at least three latent classes (see the examples in (ref)).
Model selection when parameters are estimated using variational inference remains an open problem in general. One may accept the ELBO as a good approximation of the marginal likelihood and use it as the basis of model selection, but this has not been justified in theory (blei2017variational). Other approaches include sequential selection (sato2001online), cross validation (nott2012regression) and Generalized Evidence Bounds (chen2018variational). For the purpose of this paper, we take a more practical approach by using the standard train-test split and examining the approximated loglikelihood and ELBO on the test set to obtain a conservative gauge of goodness-of-fit. Examples are given in the real data analysis in (ref).
\afterpage{ \subfile{tables/simulation-1-summary-table} \subfile{tables/simulation-2-summary-table} }
\afterpage{
}
In this section, we present two simulation studies in order to numerically illustrate the estimation algorithm described in (ref). Our goal is to examine whether the proposed algorithm can correctly estimate the model parameters and make reasonable inference about the posterior distribution of random effects.
In both simulation studies, we consider a sample size of 50,000 observations where the covariate $\bm{x}_i$ consists of an intercept term $x_{i0} = 1$ and an indicator variable $x_{i1} \sim \textrm{Bernoulli}(0.5)$. The one-dimensional response $y_i$ is generated from a mixture of Gamma distributions, where the number of mixture components is two for case I and three for case II. As for the random effects, Simulation I contains one single level $\{w_{1}^{(s)}\}_{s=1, 2, \dots, 200}$ with $S_1=200$ levels randomly assigned to all observations. Simulation II has a more complex nested structure with two random effects such that $\{w_{1}^{(s)}\}_{s=1, 2, \dots, 200}$ has $S_1=200$ levels and $\{w_{2}^{(s)}\}_{s=1, 2, \dots, 2000}$ has $S_2=2000$ levels, both randomly assigned to all observations.
The true and fitted model parameters are summarized in (ref) and (ref), while (ref) and (ref) visualize the simulated versus fitted random effects, latent class probabilities and the marginal distribution of response. Overall, we observe that the estimation algorithm is able to recover the true model parameters $(\bm{\alpha}, \bm{\beta}, \bm{\Psi})$ to a reasonable degree, which results in a close fit to the marginal distribution of the response variable, as indicated by the fitted density and the histogram of simulated data.
In addition to the response, we also investigate the simulated values versus the fitted posterior distribution of the random effects. We examine how well the approximated posterior credible intervals (CI) at different levels (90%, 95%, 97.5% and 99%) can recover the simulated true values of the random effects, which are summarized in the same set of tables and figures. In Simulation I, the random effects are well-recovered, and the plot of 95%-CI shows a high level of alignment with the simulated true values of random effects. The results in Simulation II are noticeably worse but still acceptable, considering the added noises from two different types of random effects and much fewer observations per level. For example, each $w_{2}^{(s)}$ has 25 observations per factor on average, compared with 250 for $w_{1}^{(s)}$, which results in $w_{2}^{(s)}$ having much wider 95%-CI and more cases where the posterior CI does not recover the true simulated values. Still, our algorithm is able to reasonably recover the latent class probabilities in both simulations, which contributes to the nice fit to the marginal distribution of the response.
For comparison, we have experimented the BLUP procedure outlined in yau2003finite and ng2007extension for similar MoE models with random effects. However, the BLUP procedure fails to recover the realizations of random effects, and we arrived at fitted models without any random effects (i.e. all of $\bm{w}$ have degenerated to zero). Compared with the alternative of MCMC methods, our estimation algorithm is highly efficient in terms of computational cost. We have only used 50 iterations of ECM in both simulations to produce the results above, where $M=5$ samples of random effects are used in each ECM iteration for numerical evaluation such as (ref). When implemented as a modification of the LRMoE.jl package, the computation time is 5 minutes for Case I and 15 minutes for Case II on a modern MacBook. We have experimented with standard MCMC algorithms as a benchmark, but our implementation did not converge within an acceptable time frame. A comparable study for MCMC methods in GLMM can be found in hadfield2010mcmc, where the author analyzed a dataset with 828 observations and one single random effects with 106 levels. Their example converges with 60,000 total iterations, 10,000 iterations of burn-in and a thinning interval of 25. Considering our simulation studies are done on a much larger scale, we therefore reasonably expect our VI algorithm to be much more efficient than a comparable implementation of MCMC, in terms of the number of ECM iterations needed to converge. This computational advantage will be more significant in our real data analysis presented in the next section, where the number of levels of the random effect is more than a few thousands and our algorithm typically converges within two days after a few hundreds iterations of ECM.
\afterpage{ \subfile{tables/real-data-overview-table} }
In this section, we apply the Mixed LRMoE model to a real automobile insurance dataset for a posteriori risk classification and ratemaking, and then compare its performance with a number of benchmark models. More specifically, we will investigate whether the Mixed LRMoE model can outperform benchmark models like GLM, GLMM and LRMoE without random effects in terms of goodness-of-fit. We will also investigate whether the Mixed LRMoE produces reasonable results for a posteriori risk classification and ratemaking, that is, policyholders who made claims in the past should generally be considered riskier and should be assigned a higher a posteriori premium.
The dataset contains the Bodily Injury (BI) claim history of 15,492 unique policyholders from policy years 2014 to 2019 (92,952 records in total) of a major North American automobile insurer. For illustration purposes, we have filtered for policyholders with exactly 6 years of history from 2014 to 2019. Practical issues, such as working with policyholders with a shorter history, or people with fractional policy year exposures, are trivial to address in the same modelling framework. Since we are only working with a one-dimensional response, it will be represented by $y_i$ in this section. The description of available covariates $\bm{x}_i$ and the summary statistics of the response $y_i$ are given in (ref). We observe the loss distribution has significant zero inflations and a heavy tail. We divide the entire dataset into training (2014--2018, or 12,394 unique policyholders with 61,968 records), validation (2014--2018, or 3,098 unique policyholders with 15,490 records) and testing (2019, for all 15,492 policyholders in training and validation) sets. Our goal is to fit various model candidates to the 5-year training period and then conduct a posteriori risk classification and rakemaking for the 1-year testing period. The validation set contains a 20% of the unique policyholders randomly selected from the 5-year training period, which is used for selecting the number of latent classes in the (Mixed) LRMoE models.
For illustration purposes, we will model the total amount of loss per year. As a benchmark, we will consider various combinations of GLM and GLMM against which we compare the proposed Mixed LRMoE model. For these benchmark models, we assume independence between claim frequency and severity. We use a probability mass $\delta_{i}(\bm{x}_i)$ at zero for no occurrence of claims and a continuous distribution $g_{i}(y_{i}|\bm{x}_i)$ for the total loss amount given there is at least one claim. Consequently, using $I_{\{y_{i}=0\}}$ and $I_{\{y_{i}>0\}}$ as indicators for the occurrence of claims, the distribution of total loss of policyholder $i$ is given by
where both $\delta_{i}(\bm{x}_i)$ and $g_{i}(y_{i}|\bm{x}_i)$ may be modelled by either GLM or GLMM. In the case of GLMM, we will add policyholder-level random effects with $15,492$ levels which corresponds to the number of unique policyholders in the training dataset.
For the models to investigate, we will consider (mixed) LRMoE with zero-inflated (ZI) lognormal expert functions. With the expert functions fixed, we only need to select the number of latent components for both LRMoE and Mixed LRMoE. We have selected a 4-component LRMoE and a 5-component Mixed LRMoE based on the Akaike Information Criterion (AIC) calculated on the validation dataset. For comparison with LRMoE without random effects, we also include a 4-component Mixed LRMoE in the following discussion.
\afterpage{ \subfile{tables/real-data-model-benchmark} \subfile{tables/real-data-model-LRMoEs}
}
The fitted loglikelihood values of all benchmark models are summarized in (ref). As expected, the GLMM-GLMM model produces the highest loglikelihood since the policyholder-level random effects are used twice. The combinations of GLMM-GLM and GLM-GLMM offer slightly worse fit to data, followed by the GLM-GLM model without any random effects. (ref) summarizes the loglikelihood of the 4-component LRMoE and two Mixed LRMoE models. We see that both mixed LRMoE models offer much better fit to data in terms of loglikelihood on training and testing datasets, and both outperform the LRMoE model without random effects. This demonstrates the flexibility of mixed LRMoE as well as the advantage of incorporating policyholder-level random effects for more accurate modelling of the loss distribution. As for penalization on model complexity, we also include the number of parameters for all model candidates in the tables. It is clear that the Mixed LRMoE models outperforms all benchmark models in terms of AIC on the training set, while they have marginally worse AIC values on the testing set. However, as will be evident in the next subsection, this added model complexity greatly improves a posteriori risk classification and ratemaking, which is the ultimate goal in this context.
Besides loglikelihood values, we also look at how each model candidate fits the probability of claim and the distribution of positive losses. For the probability of claim, all model candidates offer very similar fitting performance. On the training period, all models are able to fit the observed claim probability 0.978996 to the fourth decimal place. However, on the testing period where the observed claim probability is 0.988793, all models candidates have produced a slightly lower prediction, ranging from 0.979844 to 0.980005 (or $-0.91\%$ to $-0.89\%$ of relative error). Meanwhile, the (Mixed) LRMoE models have provided a better fit to the distribution of positive losses, as indicated by (ref) which compares the fitted densities against the empirical distribution. Most notably, the (Mixed) LRMoE models have successfully captured the multimodality in the tail, while GLM and GLMM only fit a unimodal density to the entire distribution of positive losses.
For both the claim probability and the distribution of positive losses, we have observed a potential data drift for testing period. In particular, the claim probability increases in 2019 compared with previous years, while the distribution of positive losses also appears to have changed a little, but the latter is only based on roughly 170 losses observed in the testing period.
\afterpage{ \subfile{tables/real-data-prob-changes} }
\afterpage{
\subfile{tables/real-data-premium-table}
}
For insurance pricing purposes, it is crucial that policyholders' claim history is adequately incorporated in the calculation of premium at policy renewal. In short, higher risks, as reflected by the occurrence of claim and/or higher claim amounts, should lead to a higher a posteriori premium. In this subsection, we compare the model performance in terms of a posteriori risk classification and ratemaking.
For risk classification, the latent classes in (mixed) LRMoE models can be naturally interpreted as different clusters of policyholders based on their risk profile. To compare how risk classification is affected by claim history, we categorize all policyholders into two groups: those with at least one claim and those without any claim during 2014--2018, and summarize their latent class probabilities in (ref). Most notably, with the addition of random effects, the Mixed LRMoE models are able to strongly distinguish risky policyholders who have at least one claim in the past, by assigning almost double the probability to the riskiest latent class. Meanwhile, the LRMoE model without random effects only suggests a slight increase in the risky class probability based solely on covariate information, given the independence assumption for observations across different policy years.
Different decisions in a posteriori risk classification will also lead to differences in ratemaking. For a posteriori ratemaking, we calculate the premium for policy renewals in year 2019 based on the posterior distribution given the claim history in 2014--2018. For illustration purposes, we only consider the pure premium which is equal to the probability of claim multiplied to the expected positive mean loss amount.
On a higher level, we investigate all policyholder based on the same grouping (with and without claims in 2014--2018). The distributions of the predicted posterior premium are shown in (ref) for all model candidates. For models without random effects, i.e. GLM-GLM and LRMoE, the predicted distributions of posterior premium for the two groups appear to be highly overlapping, which indicates that fixed effects alone cannot distinguish policyholders based on claim history. For benchmark models with random effects, namely GLM-GLMM, GLMM-GLM and GLMM-GLMM, there appears to be some difference between the two groups, whereby some policyholders with claim history will have a higher predicted premium. Most notably, the two Mixed LRMoE models show much larger differences between the distributions of predicted premium, which better captures the riskiness of policyholders reflected by their claim history.
On a more detailed level, (ref) summarizes the predicted posterior premium, based on the two groups above in addition to the relative size of incurred total losses. We observe that both Mixed LRMoE models, as well as benchmark GLMM-GLM and GLMM-GLMM, heavily penalizes policyholders who have at least one claim, as shown by the additional premium loadings. For people with claims, only the Mixed LRMoE with five components has further provided a correct ordering of the predicted posterior premium, that is, people with larger incurred claims typically have higher premium at policy renewal, which indicates better performance in a posteriori risk classification and ratemaking. This is because only the 5-component Mixed LRMoE has adequately captured the multimodality in the tail of the positive loss distribution, as indicated by (ref).
For both a posteriori risk classification and ratemaking discussed above, we have primarily focused on differentiating policyholders based on the occurrence of claims and the claim sizes when applicable, whereby the Mixed LRMoE models are shown to have effectively incorporated such information. However, we can still observe the effects of a priori information, i.e. policyholder covariates, when determining the a posteriori premium. Most notably, in (ref), there is a good level of overlap between the histograms of the predicted premium for people with and without claim history, even for all model candidates with random effects. For example, certain policyholders with claim history (lower end of the orange histogram) would still be charged a lower premium than some policyholders without claim history (upper end of the blue histogram), which should be attributed to covariates such as the inherent risk level of certain age groups or the collision rating of a particular group of vehicles.
\afterpage{
\subfile{tables/real-data-Gini-comparison} }
Finally, we examine the model performance using the Gini Index as a measurement of adequacy for insurance risk scoring (see e.g. frees2011summarizing). We first plot the Ordered Lorenz Curve in (ref) for both the training and testing sets, where the $x$-axis represents the cumulative percentage of premium and $y$-axis represents the cumulative percentage of the incurred losses during the training or testing period. The corresponding Gini index values for all model candidates, calculated as twice the area between the Ordered Lorenz Curve and the Line of Equality (45-degree line), as well as their estimated standard error, are summarized in (ref).
On the training set, we see the two mixed LRMoE models have produced Ordered Lorenz Curves farthest from the Line of Equality as well as the largest Gini Index values, which indicates a high degree of differentiation between low- and high-risk policyholders based on their claim history. We also note the second best models in terms of Gini index are GLMM-GLM and GLMM-GLMM, which means the probability of claim may potentially be a more important determinant of policyholders' risk profile compared with claim severity.
However, on the testing set, all model candidates perform quite similarly, and the two mixed LRMoE models do not outperform the classical models. In fact, the estimated standard errors of the Gini Index suggest we cannot conclude whether the performances of all model candidates are significantly different from each other. This may have been caused by the small number of incurred claims, as observed in (ref), but a more important factor might be the potential data drift in year 2019 with a slightly lower claim probability and changed distribution of positive losses. However, such unprecedented data drift is outside the scope of what statistical and predictive models can address based on historical data only.
In this paper, we have proposed to incorporate policyholder-level random effects in a flexible regression framework, called the Mixed LRMoE, which is then applied to the Bonus-Malus problem. Although the addition of random effects has resulted in an intractable marginal likelihood function of the model, we have developed a stochastic variational ECM algorithm for efficient estimation of model parameters and inference of the posterior of random effects, which are crucial for updating policyholders' risk profile based on their claim history. Our numerical simulation and real data analysis have demonstrated the potentials of Mixed LRMoE as a powerful tool for more accurate insurance loss modelling and better a posteriori insurance risk classification and ratemaking. While our current work has already shown promising results, one may consider the following extensions and directions for future work.