Introduction to actuaRE"

In actuaRE: Handling Hierarchically Structured Risk Factors using Random Effects Models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)

wzxhzdk:1

In this document, we give you a brief overview of the basic functionality of the actuaRE package. For a more detailed overview of the functions, you can consult the help-pages. Please feel free to send any suggestions and bug reports to the package author.

Handling hierarchically multi-level factors using random effects models

Multi-level factors (MLFs) are nominal variables with too many levels for ordinary generalized linear model (GLM) estimation [@Ohlsson]. Within the machine learning literature, these type of risk factors are better known as high-cardinality attributes [@Micci2001]. This package focuses on MLFs that exhibit a hierarchical structure and a typical example hereof, within workers' compensation insurance, is the NACE code. In our illustration, we work with a hierarchical MLF that has two hierarchical levels: industry and branch. Figure 1 visualizes this hierarchical structure with a hypothetical example.

knitr::include_graphics("./HierarchicalStructureAdj.png")

Which random effects models can we fit using actuaRE?

With the current version of the actuaRE package, you are able to fit random effects models with the following functional form

\begin{align} g(E[Y_{ijkt} | U_j, U_{jk}]) &= \mu + \boldsymbol{x}{ijkt}^\top \boldsymbol{\beta} + U_j + U{jk} \&= \zeta_{ijkt}.\ \end{align}

Here, $Y_{ijkt}$ denotes the loss cost of risk profile $i$ (based on the company-specific risk factors) operating in branch $k$ within industry $j$ at time $t$. We calculate the loss cost as

\begin{align} Y_{ijkt} = \frac{Z_{ijkt}}{w_{ijkt}} \end{align} where $Z_{ijkt}$ denotes the total claim cost and $w_{ijkt}$ is an appropriate volume measure. $g(\cdot)$ denotes the link function (for example the identity or log link), $\mu$ the intercept, $\boldsymbol{x}{ijkt}$ the company-specific covariate vector and $\boldsymbol{\beta}$ the corresponding parameter vector. With the model parameters $\mu$ and $\boldsymbol{\beta}$ we capture the company-specific effects. To assess the effect of the hierarchical MLF, we introduce the random effects $U_j$ and $U{jk}$ which capture the unobservable effects of the industry and the branch in which the company operates. $U_j$ denotes the industry-specific deviation from $\mu + \boldsymbol{x}{ijkt}^\top \boldsymbol{\beta}$ and $U{jk}$ denotes the branch-specific deviation from $\mu + \boldsymbol{x}{ijkt}^\top \boldsymbol{\beta} + U{j}$. We assume that the random industry effects $U_j$ are independent and identically distributed (i.i.d.) with $E[U_j] = 0$ and $Var(U_j) = \sigma_{I}^2$. Similarly, the random branch effects $U_{jk}$ are assumed to be i.i.d. with $E[U_{jk}] = 0$ and $Var(U_{jk}) = \sigma_{B}^2$.

This package offers three different estimation methods to estimate the model parameters:
- Hierarchical credibility model [@JewellModel]
- Combining the hierarchical credibility model with a GLM [@Ohlsson2008]
- Mixed models [@Molenberghs2005]

Just the code please

Example data set

To illustrate the functions, we make use of two different data sets. We illustrate the hierarchical credibility model of Jewell [@JewellModel] using the Hachemeister [@Hachemeister] data set. The other functions make use of the dataCar data set.

Hierarchical credibility model

To estimate the parameters using the hierarchical credibility model, we use the function hierCredibility. By default, the additive hierarchical credibility model [@Dannenburg] is fit

\begin{align} E[Y_{ijkt} | U_j, U_{jk}] &= \mu + U_j + U_{jk}. \end{align}

capture.output(library(actuaRE), file = tempfile()) # suppress startup message
data("hachemeisterLong")
fitHC = hierCredibility(ratio, weight, cohort, state, hachemeisterLong)
fitHC

To fit the multiplicative hierarchical credibility model [@OhlssonJewell] \begin{align} E[Y_{ijkt} | \widetilde{U}j, \widetilde{U}{jk}] &= \tilde{\mu} \ \widetilde{U}j \ \widetilde{U}{jk} \end{align} you have to specify type = "multiplicative".

fitHCMult = hierCredibility(ratio, weight, cohort, state, hachemeisterLong, type = "multiplicative")
fitHCMult

To get a summary of the model fit, we use the summary function.

summary(fitHC)

To obtain the fitted values, we use the fitted function

fitted(fitHC)

and we use ranef to extract the estimated random effects.

ranef(fitHC)

We can inspect the estimated random effects using the function plotRE.

ggPlots = plotRE(fitHC, plot = FALSE)
ggPlots[[1]]
ggPlots[[2]]

To obtain predictions for a new data frame, we use the predict function.

newDt = hachemeisterLong[sample(1:nrow(hachemeisterLong), 5, F), ]
predict(fitHC, newDt)

Combining the hierarchical credibility model with a GLM

To allow for company-specific risk factors, we extend the multiplicative hierarchical credibility model to \begin{align} E[Y_{ijkt} | \widetilde{U}j, \widetilde{U}{jk}] &= \tilde{\mu} \ \gamma_{ijkt} \ \widetilde{U}j \ \widetilde{U}{jk} = \gamma_{ijkt} V_{jk} \end{align} where $\gamma_{ijkt}$ denotes the effect of the company-specific covariates. To estimate this model using Ohlsson's GLMC algorithm [@Ohlsson2008], we use can either use the function hierCredGLM or hierCredTweedie. hierCredGLM allows the user to specify the power parameter $p$. Conversely, hierCredTweedie estimates the power parameter $p$ along with the other parameters using the cpglm function from the cplm package.

data("dataCar")
fit = hierCredGLM(Y ~ area + (1 | VehicleType / VehicleBody), dataCar, weights = w)
summary(fit)

We use the same syntax as used by the package lme4 to specify the model formula. Here, (1 | VehicleType / VehicleBody) specifies a random effect $U_j$ for VehicleType and a nested random effect $U_{jk}$ for VehicleBody. We extract the estimated parameters using fixef (company-specific effects) and ranef (random effects).

fixef(fit)
ranef(fit)

In addition, the same functions as before can be used.

head(fitted(fit))
predict(fit, newdata = dataCar[1:2, ], type = "response")
ggPlots = plotRE(fit, plot = FALSE)

Mixed models

Alternatively, we can rely on the mixed models framework [@Molenberghs2005] to estimate the model parameters. Here, we can use the \code{\link[cplm]{cpglmm}} function to estimate a Tweedie generalized linear mixed model. Fitting the model, however, takes quite some time. We can speed up the fitting process by providing some initial estimates and this is exactly what the tweedieGLMM function does! Nonetheless, even with the initial estimates the fitting process does take some time (approximately 5 minutes using Windows 10 with an intel i7 and 32 gigabytes of RAM).

fitGLMM = tweedieGLMM(Y ~ area + (1 | VehicleType / VehicleBody), dataCar, weights = w, verbose = TRUE)

Balance property

For insurance applications, it is crucial that the models provide us a reasonable premium volume at portfolio level. Hereto, we examine the balance property [@Buhlmann2005][@Wuthrich] on the training set. That is, \begin{equation} \begin{aligned} \sum_{i, j, k, t} w_{ijkt} \ Y_{ijkt} &= \sum_{i, j, k, t} w_{ijkt} \ \widehat{Y}{ijkt}\ \end{aligned} \end{equation} where $i$ serves as an index for the tariff class. GLMs fulfill the balance property when we use the canonical link (see [@Wuthrich]). For LMMs and hence, the hierarchical credibility model this property also holds. Conversely, most GLMMs do not have this property. To regain the balance property, we introduce a quantity $\alpha$ \begin{equation} \begin{aligned} \alpha &= \frac{\sum{i, j, k, t} w_{ijkt} \ Y_{ijkt}}{\sum_{i, j, k, t} w_{ijkt} \ \widehat{Y}_{ijkt}}\ \end{aligned} \end{equation} which quantifies the deviation of the total predicted damage from the total observed damage. In case of the log link, we can then use $\alpha$ to update the intercept to $\hat{\mu} + \log(\alpha)$ to regain the balance property.

By default, the intercept is updated when fitting models using hierCredGLM, hierCredTweedie and tweedieGLMM. If you do not wish to update the intercept, you can set the argument balanceProperty = FALSE.

fitnoBP  = hierCredGLM(Y ~ area + (1 | VehicleType / VehicleBody), dataCar, weights = w, balanceProperty = F)
yHatnoBP = fitted(fitnoBP)
w        = weights(fitnoBP, "prior")
y        = fitnoBP$y

fitBP  = hierCredGLM(Y ~ area + (1 | VehicleType / VehicleBody), dataCar, weights = w, balanceProperty = T)
yHatBP = fitted(fitBP)

sum(w * y) / sum(w * yHatnoBP)
sum(w * y) / sum(w * yHatBP)

Alternatively, you can use the build-in function BalanceProperty. You can use this function with any object that has the slots fitted, weights and y.

BalanceProperty(fitnoBP)
BalanceProperty(fitBP)

References

actuaRE documentation built on March 31, 2023, 10:16 p.m.