Generalized linear models (GLMs) and generalized linear mixed models (GLMMs) are fundamental tools for predictive analysis. These statistical tools have motivated many modern insurance-pricing techniques (Nelder and Verrall, 1997; Antonio and Beirlant, 2007; Wuthrich, 2019). Meanwhile, the credibility theory is a field within actuarial science that focuses on calculating risk premiums. It originated from Mowbray (1914) and Whitney (1918), who aimed to determine a premium that balanced the experience of an individual policyholder with that of a group of similar policyholders (Nelder and Verrall, 1997). While GLMMs and credibility theory are both widely used in risk classification in non-life insurance sectors, such as auto insurance, practitioners in these sectors tend to be more familiar with credibility theory.

GLMMs have several advantages over credibility theory. They offer a wider range of models and can be easily fitted and tested using standard statistical software. As a result, GLMMs provide additional modeling options for various actuarial problems, such as premium ratings and claims reserves. Furthermore, it is well known that credibility theory is closely related to the framework of GLMMs (Nelder and Verrall, 1997). While the insurance literature has some initial discussion on credibility theory based on GLMMs (Ohlsson, 2008; Lai and Sun, 2012), most assume above-average knowledge of the statistical aspects of GLMM, which makes its application inaccessible to traditional insurance practitioners.

This study aims to explain GLMs and GLMMs in plain actuarial language and demonstrate how to obtain a simple and fair premium. The remainder of the article is organized as follows. Section 2 introduces the notations used. Section 3 provides the general settings of risk classification procedures for insurance using the language of GLMs. Sections 4 and 5 explain the modeling and prediction methods for GLMs and GLMMs, respectively. Section 6 provides the theoretical link between GLMMs and credibility theory. Finally, the results are discussed in Section 7.

We denote ℕ, ℕ₀, ℝ, and ℝ⁺ as the set of natural numbers, the set of non-negative integers, the set of real numbers, and the set of positive real numbers, respectively. We define N_it as the number of claims in the t^th policy year for the i^th policyholder. We denote the realization of N_it as n_it. Actuarial science literature often refers to N_it as the frequency of insurance claims.

We differentiate a fixed effect from a random effect. The former is based on the policyholders’ observed risk characteristics, whereas the latter is based on the policyholders’ unobserved risk characteristics after controlling for fixed effects. Let the data matrix X_i of the i^th policyholder have a dimension of τ × (m + 1), where τ is the number of observed years and m is the total number of explanatory variables. Denote

where X_it := (1, X_it1, . . ., X_itm) is the (m+1)-dimensional row vector of the a priori risk characteristics of the i^th policyholder in the t^th year. Let γ_i denote the policyholder-specific random effect of the i^th policyholder. For simplicity, we assume that the random effect γ_i does not have the subscript t, because we assume that the a priori risk characteristics remain constant over time. Under GLM and GLMM, the distribution within the exponential dispersion family (EDF) has been addressed by McCullagh and Nelder (1989) and Breslow and Clayton (1993). EDF can be viewed as an abstraction of many well-known distributions, including normal, gamma, and Poisson distributions. This provides a unified approach to statistical modeling. However, owing to its abstract form, it may be difficult to understand intuitively. Therefore, we do not explore the benefits of a generalized form of EDF in this study, nor do we attempt to identify the best parametric assumption for the distribution that fits the data. Instead, we focus primarily on Poisson distribution because the general idea can easily be extended to other distributions within the EDF.

We use Pois(λ) to denote a Poisson distribution with mean λ. Using gamma(μ, ψ), we denote a gamma distribution whose mean and variance parameters are μ and μ²ψ. Gamma(μ, ψ) can also be interpreted as a gamma distribution with the shape parameter 1/ψ and scale parameter μψ. We use N(μ,σ²) to denote a normal distribution, whose mean and variance are given by μ and σ², respectively. For Y ~ N(μ,σ²), we define the distribution of exp(Y) as lognormal(μ,σ²).

In terms of data, for each policyholder i = 1, . . ., z in τ policy years, with m explanatory variables, the dataset has the following format:

where z is the total number of policyholders and x_it = (x_it1, x_it2, . . ., x_itm). An example of a typical dataset is listed in Table 1.

First, we are interested in a fair premium, defined as , for the i^th policyholder in calendar year τ^* = τ + 1. Throughout this study, although we do not specify so in the mathematical notation, we assume that the explanatory variables are always known and conditioned. Hence, for example, the fair premium should be understood as

The fair premium for each policy year is determined solely by the fixed effect. In insurance, a fair premium based on the fixed effect is called the a priori rate and the risk classification process based on the fixed effect is called the a priori risk classification procedure. The premium determined by the a priori risk classification procedure can be further updated based on the claim history as follows:

which is called the a posteriori rate. This premium adjustment mechanism based on claim history is called the a posteriori risk classification procedure. Typically, in insurance, the risk characteristics for a year are known at the beginning of each year. Hence, in predicting the premium in the τ^*th year, the explanatory variables x_i1, . . ., x_iτ^* are assumed to be known. In the following sections, we examine the a posteriori risk classification procedure in (3.2) under various model assumptions.

GLMs are a flexible method that extend the fixed-effect linear regression model to response variables, which are not normally distributed. The expected value of N_it depends on the observed risk characteristics with the following relation:

where X_itβ is the linear predictor and η(·) is the link function, which relates to a linear predictor.

4.1. Modelling and estimation of the GLM

Model 1. (Poisson GLM)Consider z different policyholders with the observation given in the form of (3.1) and assume mutual independence among policyholders. Furthermore, we assume the following joint distribution of (N_i1, . . ., N_iτ^*):

• The distribution of the i^thpolicyholder’s frequency N_itin the t^thyear is specified as

  Nit~Pois (Λit),

  where

  (4.1)Λit:=exp (Xitβ)=exp (β0+β1Xit1+⋯+βmXitm)

  with the coefficient parameters (fixed effects) β := (β₀, β₁, . . ., β_m)′. Here, we assume the link function η to be a log function that gives the exponential expression in (4.1).

• N_i1, . . ., N_iτ^*are independent.

In Model 1, the only unknown parameters are the coefficients β shared by each policyholder. Hence, once the coefficients β are known, the conditional distribution of frequency for each policyholder in the future calendar year τ^* = t + 1 can be determined exactly as

with the a priori rate for the i^th policyholder in calendar year τ^* given by .

We assume that we have a dataset in the form of (3.1) with z policyholders and that each policyholder is observed for τ years. Then, the log-likelihood function for Model 1 is written as

where f_β (n_it) is the probability mass function of the Poisson distribution given by

where Λ_it = exp(X_itβ). For the estimation of unknown coefficients β, typically, maximum likelihood estimation (MLE) is considered and the estimator of β is defined as

4.2. Prediction in the GLM

In this subsection, we predict the a posteriori rate in (3.2), which is the premium in the τ^*th year. given the claim history up to the τ^th year. Under the settings in Model 1, the a posteriori rate defined in (3.2) can be obtained as

based on the assumption of independence in Model 1. The a posteriori rate in (4.2) coincides with the a priori rate. As (4.2) shows, the a posteriori rate is not a function of the past claim history under the setting in Model 1 and there is no premium adjustment by the claim history n_i1, . . ., n_iτ.

4.3. Data example

Throughout this study, we use the dataset for collision coverage on vehicles collected by the Wisconsin Local Government Property Insurance Fund (LGPIF) from 2006 to 2011 (Frees et al., 2016). Collision coverage covers the impact of a vehicle on an object, the impact of a vehicle on an attached vehicle, or the overturning of a vehicle. The datasets contain 497 government entities. Each policyholder has six levels of entity types and average coverage over time as explanatory variables as shown in Table 2. For further explanation of the dataset, refer to Oh et al. (2020).

To fit the GLM in Model 1, we can use the glm function in R, whose form is given by

 > glm(formula, family=familytype(link=linkfunction), data=dataset,. . .). 

The specific form of the function under the setting in Model 1 is given by

 > glm(Frequency~factor(entityType)+avgCoverage, family=poisson(link=log), data=mydata). 

If the categorical variable (entity type) is presented as a dummy variable (e.g., TypeCity, TypeCounty) in a dataset, the formula specifying the fixed effects part of the function is modified as follows:

 > glm(Frequency~TypeCity+TypeCounty+· · · +avgCoverage, family=poisson(link=log), data=mydata). 

The estimation results of this model are summarized in Table 3. The asterisk (*) signs indicate that the parameters are significant at the significance level of 0.05. Based on the table, the a posteriori rate, which is equivalent to the a priori rate, can be obtained as

Hence, the estimated parameter β̂₁, for example, can be interpreted as the effect of the entity type of the county compared to the baseline, miscellaneous entity type, in the logarithm of the estimated a posteriori rate for the same average coverage. The table indicates that the county entity type is expected to incur more accidents than any other entity type having the same amount of coverage. Within the same entity type, a higher average coverage corresponds to a higher expected accident count.

Moreover, the third policyholder in the dataset, for example, is in a county and has an average coverage of 2.495, so that under Model 1, the estimated expectation is

GLMMs are a popular method to model correlated data; they extend GLMs by including random effects in the predictor. The conditional expected value of N_it|γ_i depends on the observed risk characteristics using the following relation:

where γ_i is the random effect and η(·) is a link function that converts into a linear predictor X_itβ. The random effect γ_i measures the difference between the average score of policyholder i and the average score in the portfolio.

5.1. Modelling and estimation of the GLMM

In Model 1, we assume that the selected explanatory variable X_i1, . . ., X_iτ^* can explain the distribution of the frequency. However, in reality, it is almost impossible to describe human behavior and contributors are often not observed or overlooked. These unobserved components can be used in a regression setting, as shown in the following example:

Example 1. Consider Model 1 with three explanatory variables: X_it1, X_it2, and X_it3. If we assume that they are observed, then the a priori premium is

Now, we assume that only the first explanatory variable X_it1 is observed and the other explanatory variables, X_it2 and X_it3, are not. Then, a reasonable description of the a priori premium for frequency under this assumption is

If the unobserved explanatory variables are denoted as

where the explanatory variables are assumed invariant over time, we have the following conditional model:

By further defining R_i := exp(γ_i), the conditional distribution in (5.1) can be represented as

Regardless of whether observed or unobserved, the explanatory variables can also be assumed to be random variables. By placing an additional distributional assumption on the unobserved explanatory variables, Example 5.1 can be formally represented as follows:

Model 2. (Poisson GLMM)Consider z policyholders with the observation given in the form of (3.1), assuming mutual independence among policyholders. Furthermore, we assume the following conditional joint distribution of (N_i1, . . ., N_iτ^* ):

• Conditional on the random effect R_i, the distribution of the i^thpolicyholder’s frequency N_itin the t^thyear is specified as

  Nit∣Ri~Pois (RiΛit),

  where

  (5.2)Λit:=exp (xitβ)=exp (β0+β1Xit1+⋯+βmXitm)

  with the coefficient parameters (fixed effects) β := (β₀, β₁, . . ., β_m)′. Here, we assume the link function η to be a log function that gives the exponential expression in (5.2).

• The random effect R_ifollows distribution function G with a parameter ψ.

• Furthermore, we assume that N_i1, . . ., N_iτ^*are independent conditional on R_i.

In the case of a normal random effect, that is, γ_i ~ N(0, σ₂), we have R_i ~ lognormal(0, σ₂). Otherwise, if the random effect distribution is non-normal, the mixed model is called a hierarchical generalized linear mixed model (HGLMM) (Lee and Nelder, 1996). However, we do not differentiate between GLMMs and HGLMMs in this study. Although it is possible to assume that some of the parameters in (5.2) are random effects, known as random slopes in the classical GLMM literature, we do not pursue this direction for the sake of simplicity. We refer the reader to Zimmer (2018) for information on the impact of random slopes on insurance pricing.

In Model 2, the mean of the random effect R_i is often set to be one for the identifiability issue. For example, we may assume

with the mean of 1 and the dispersion parameter ψ. We refer to Section 2 for the details of parameterization. Note that using the random effect R_i introduces longitudinal dependence among N_its. Specifically, in Model 2, the law of total covariance is derived as

However, in Model 1, we have

An important implication in the comparison of the covariance in (5.4) under Model 2 and the covariance in (5.5) under Model 1 is that the random effect induces dependence among the observations. Thus, the random effect not only explains the unobserved explanatory components but also introduces dependence among observations. Hence, because the observations of the i^th policyholder,

are not independent under Model 2, the corresponding likelihood function under Model 2 cannot be written as the product of the likelihood function. Instead, we have the following log-likelihood function for the observations:

which requires one-dimensional integration. Here, f_β(· | r_i) is the probability mass function of the Poisson distribution defined by

and g_ψ is the probability density function of the random effect in (5.3). From a purely mathematical perspective, the introduction of a random effect in Model 2 is for modeling dependence within the observations. In this sense, the random-effect model is linked to the copula model (Frees and Valdez, 1998; Embrechts, 2009; Krupskii and Joe, 2013).

Next, for data in the form of (3.1), with z policyholders and each policyholder observed for τ years, the log-likelihood function for Model 2 is written as

For the estimation of unknown coefficients β and parameter ψ, MLE is usually considered, and the estimators of β and ψ are typically defined as

5.2. Prediction in the GLMM

The goal of the a posteriori risk classification is to predict the net premium for each policyholder. In Model 2, by definition, we have the following conditional expectation:

Because the random effect for the i^th policyholder, R_i, is an unknown random variable, the conditional expectation in (5.8) cannot be used directly to predict the a posteriori rate in (3.2). However, in combination with the prediction of a random effect,

the conditional expectation in (5.8) can calculate the a posteriori rate in (3.2) as follows:

Hence, the prediction of the a posteriori rate in Model 2 is reduced to the estimation of predicting the random effect in (5.9). Note that the prediction of the random effect in (5.9) allows for a closed-form expression depending on the choice of distributional assumption for the random effect. Otherwise, the random effect in (5.9) can be predicted numerically, either based on numerical integration or a Monte Carlo simulation. The next subsection explains a specific method for predicting the random effect in (5.9).

5.3. Parametric assumptions for the random effect in the generalized linear mixed model

Model 3. Considering the settings in Model 2, we further assume that the parametric assumption for R_iin (5.3) is

In Model 3, the likelihood function in (5.6) can be analytically expressed as

where λ_it denotes the realization of Λ_it in (4.1).

Hence, the statistical estimation procedure is based on the simple optimization of (5.11), which does not require complicated numerical integration. Furthermore, since the gamma distribution is a conjugate prior of the Poisson distribution, the posterior mean of the random effect in (5.9) can be analytical obtained as follows:

from the fact that the conditional distribution R_i|n_i1, . . ., n_iτ follows a gamma distribution with shape and scale parameters given by

respectively. In the conjugate families, the posterior mean of the individual’s random effect in (5.12) can be represented by a weighted average between the prior mean of R_i (weighted by its precision) and the ratio of the total number of insurance claims across τ years for policyholder i to the sum of the respective rate parameters, (∑t=1τnit)/(∑t=1τλit) (weighted by the sum of the rate parameters):

Because later acts as a multiplier of the expectation from the fixed effects alone (see the example at the end of Section 5, under Model 3), suggests no modification in the premium. Additionally, the ratio (∑t=1τnit)/(∑t=1τλit) indicates how the observed number of claims relates to the expected number of claims over time for policyholder i. For example, if ∑t=1τnit>∑t=1τλit, then this would lead to a ratio greater than one, resulting in .

So far, we have analyzed the properties of GLMMs in Model 2, in which the random effect R_i is assumed to follow a gamma distribution. This distributional assumption is convenient because it provides an analytical expression for the likelihood function in (5.11) and the analytical conditional distribution of R_i. However, inherited from the linear mixed model, in which both the response and random variables are assumed to follow a normal distribution, the more popular choice of distributional assumption for R_i is lognormal distribution. The following model provides a GLMM model with lognormal assumption for R_i.

Model 4. We now consider the settings in Model 2 and we further assume that the parametric assumption for R_iis

We have

under the parameterization in Model 4. Unlike the gamma distributional assumption for the random effect as in Model 3, the lognormal distributional assumption for the random effect does not lead to a closed-form expression for the likelihood function in (5.6). Hence, for the statistical estimation procedure in (5.7), numerical integration is required in (5.6). Standard numerical integration methods, such as the Laplace approximation (Laplace, 1986) and Gaussian quadrature (Golub and Welsch, 1969), can efficiently approximate the likelihood function for GLMMs. We explain that the a posteriori rate in (5.10) can be approximated as an affine function of observations in Section 6.1. In the next subsection the use of statistical programming for GLMMs is explained.

5.4. Data example

To fit the GLMMs in Model 2, we consider the hglm function in the hglm package (Rönnegård et al., 2010) in R. The function is for fitting HGLMs via extended quasi-likelihood (Lee et al., 2006), which extends GLMs by relaxing the assumption that error components are independent. Its grammatical form is as follows:

> hglm(formula, family=familytype(link=linkfunction),
            random=var, rand.family=familytype(link=linkfunction),
            data=dataset,. . .). 

This is similar to the glm function in R but with additional terms for the random effect: random and rand.family. First, we specify the policyholder-specific random effect of interest as random=~1|Polic yNum. To easily specify the distribution of the random effect, we consider the conditional expectation of frequency as

and having a log link function. To use the gamma distribution for the random effect in the function, the random effect γ_i is given by γ_i = log(R_i) and R_i ~ gamma(1, ψ). Thus, we need a log link function to specify the distribution for the random effect as rand.family=gamma(link=log). Similarly, for the log-normal distribution, we can use the code rand.family=gaussian(link=identity), which implies that γ_i ~ N(0, ψ); that is, R_i ~ lognormal(0, ψ). The specific form of the function under the setting of Model 3 is

> hglm(Frequency~TypeCity+· · · +avgCoverage,
                      family=poisson(link=log),
                      random=~1|PolicyNum,
                      rand.family=gamma(link=log),
                      data=mydata).

Under the settings in Model 4,

> hglm(Frequency~TypeCity+· · · +avgCoverage,
                      family=poisson(link=log),
                      random=~1|PolicyNum,
                      rand.family=gaussian(link=identity),
                      data=mydata).

Alternative functions to fit GLMM are the glmer function of the lme4 package and the glmmPQL function of the MASS package. The
glmer function is commonly used, and it fits the model via maximum likelihood. Meanwhile, glmmPQL uses the penalized quasi-likelihood approach to evaluate the likelihood of GLMM. The functions in Model 4 are

> glmer(Frequency~TypeCity+ · · · +avgCoverage+(1|PolicyNum),
                      family=poisson(link=log),
                      data=mydata),

and

> glmmPQL(Frequency~TypeCity+· · · +avgCoverage,
                           family=poisson(link=log),
                           random=~1|PolicyNum,
                           data=mydata).

The estimation results of these models for the parameters are summarized in Table 4. The estimates for the fixed effects differ slightly depending on the distribution of the random effect; however, we have the same interpretation for the tendency in risk characteristics. In this example, the function provides an estimate of the dispersion parameters for the random effect. Because the dispersion parameter of the gamma distribution is the reciprocal of the shape parameter and has a mean of 1, ψ is implied. Furthermore, using the code print(summary(res), print.ranef = TRUE), we obtain estimates of random effects for all policyholders. Table 5 lists the estimates for selected policyholders in Models 3 and 4. This function provides a different form of estimates depending on the link function of the random effect. From the estimation results in Tables 4 and 5, we obtain the estimated conditional expectation as follows:

Here, γ̂_i changes over time, based on historical observations, whereas the corresponding expectation of GLMs in (4.3) remains invariant with respect to time and historical observations. For example, the third policyholder in the dataset is in a county and has an average coverage of 2.495; thus, under Model 3, the estimated conditional expectation is

and under Model 4, it is

Although GLMMs offer a unified approach for calculating the a posteriori insurance rate, they often lack closed-form expressions, except for a few distributional assumptions; for instance, see Model 3. In this section, we demonstrate the incorporation of the credibility theory into the framework of GLMMs by providing a closed-form approximation of the a posteriori rate.

For simplicity, we consider the settings in Model 2 with the further assumption that the risk characteristics are the same across years; that is, X_i1 = · · · = X_iτ^* . While one of the main interests of insurance is the prediction of the net premium of each policyholder defined in (5.8) using the a posteriori rate defined in (3.2), the closed form of the a posteriori rate defined in (3.2) is generally not possible. Alternatively, we may consider approximating the a posteriori rate for the i^th policyholder defined in (3.2) as the affine function of the historical observations, by the following:

where

Following the classical procedure in the Bühlmann premium (Bühlmann and Gisler, 2006), the optimal premium for the i^th policyholder in (6.1) is obtained as

where Z_τ, the Bühlmann factor, is given by

Furthermore, , and can be obtained as

from equations

where

In this study, we present an intuitive explanation of the structure of GLMs and GLMMs and their relationship with insurance pricing methods; and we explain their connection with the Bühlmann method in insurance. With recent advances in computing tools, such as deep neural networks, insurance pricing can become more accurate through the use of complex predictive expressions. However, the traditional credibility method, which seeks an intuitive representation of insurance prices, seems at odds with modern statistical techniques. An interesting future topic is to explore how to adapt modern statistical computing tools to the credibility method in order to improve accuracy while preserving its intuitive and simple form.

Here, we mainly focus on the Poisson distribution in GLM and GLMMs to model the frequency of insurance claims, which are commonly used for count data. However, based on the assumption that the variance and mean are the same in the distribution, there are limitations, such as incorrect inference when overdispersion is present in the dataset. To address overdispersion in the Poisson model, overdispersion adjustments can be made using the quasi-likelihood method (quasi-Poisson model) or an alternative distribution method, such as negative binomial distribution. We note that, with some limitation (Lee et al., 2020), the problem of overdispersion can be reduced for Poisson distribution by using the GLMM framework.