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.
108,553 characters · 37 sections · 62 citation commands
Machine Learning with High-Cardinality Categorical Features in Actuarial Applications
\normalem
The advances in machine learning (ML) over the past two decades have transformed many disciplines. Within the actuarial literature, we see an explosion of works that apply and develop ML methods for various actuarial tasks richman2021b, richman2021a. Meanwhile, there are many distinctive challenges with insurance data that are not addressed by general-purpose ML algorithms.
One prominent example is the usual presence of high-cardinality categorical features (i.e. categorical features with many levels or categories). These features often represent important risk factors in actuarial data. Examples include the occupation in commercial property insurance, or the cause of injury in workers' compensation insurance. Unfortunately, ML algorithms cannot “understand” (process) categorical features on their own.
The classical textbook approach to this problem, and also the most popular one in practice, is one-hot encoding. It turns a categorical feature of $q$ unique categories into numeric representations by constructing $q$ binary attributes, one for each category; for example, a categorical feature with three unique categories will be represented as $[1, 0, 0]$, $[0, 1, 0]$ and $[0, 0, 1]$.
One-hot encoding works well with a small number of independent categories, which is the case with most applications in the ML literature hastie2009. However, issues arise as cardinality expands: (i) The orthogonality (i.e. independence) assumption of one-hot encoding no longer seems appropriate, since the growing number of categories will inevitably start to interact as the feature space gets more and more crowded, (ii) The resultant high-dimensional feature matrix entails computational challenges, especially when used with already computationally expensive models such as neural networks; (iii) The often uneven distribution of data across categories makes it difficult to learn the behaviour of the rare categories.
For a motivating example, consider the workers' compensation scheme of the stateofnewyork2022. In the four million claims observed from 2000 to 2022, the cause of injury variable has 78 unique categories, with more than 200,000 observations (about 5%) for the most common cause (lifting) and less than 1000 observations (0.025%) for the 10 least common causes (e.g. crash of a rail vehicle, or gunshot). The uneven coverage of claims presents a modelling challenge and has been documented by actuarial practitioners pettifer2012. Furthermore, the cause of injury variable alone may not serve well to differentiate claim severities. Injuries from the same cause can result in vastly different claims experiences. It seems natural to consider its interaction with the nature of injury and part of body variables, each with 57 reported categories. Exploring $4446~(78 \times 57)$ interaction categories is clearly infeasible with one-hot encoding.
The example above highlights that one-hot encoding is an inadequate tool for handling high-cardinality categorical features. A few alternatives exist but also have their own drawbacks, as listed below. It is worth pointing out that while the following approaches appear very distinct, they share one common feature. They all encourage the categories to “communicate” with each other, in the sense that the information learned from one category should also help with learning of the others.
Some recent work has appeared in the ML literature aiming to extend GLMMs to take advantage of the more capable ML models, by embedding a ML model into the linear predictor of the GLMMs. The most notable examples include GPBoost sigrist2022 that combines GLMMs with a gradient boosting algorithm, and LMMNN simchoni2022 that combines linear mixed models (LMM) with neural networks. Unfortunately, these models present their own limitations in the face of insurance data: GPBoost assumes overly optimistic random effects variance, and LMMNN is limited to a Gaussian-distributed response. While the Gaussian assumption offers computational advantages (in terms of analytical tractability), it is often ill suited to the more general distributions that are required for the modelling of financial and insurance data.
None of the existing techniques appears satisfactory for the modelling of high-cardinality categorical features. This paper seeks an alternative approach that more effectively leverages information in the categories, offers greater transparency than entity embeddings, and demonstrates improved predictive performance. In the next section, we introduce our approach.
Our work presents two main contributions.
First, we take inspiration from the latest ML developments and propose a generalised mixed effects neural network called GLMMNet. The GLMMNet is an extension to the LMMNN proposed by simchoni2022. In a similar vein, the GLMMNet fuses a deep neural network to the GLMM structure: The Network component of the model provides it with the flexibility required to capture complex non-linear relationships in the data, and the GLMM-like component of the model allows a probabilistic interpretation and offers full transparency on the categorical variables (through the random effects estimates). In contrast with the Gaussian-based LMMNN, the key enhancement of the GLMMNet is its ability to model the entire class of exponential dispersion (ED) family distributions. The proposed extension requires a significant effort as it involves moving away from analytical results. The flexibility it brings, however, makes the GLMMNet approach widely applicable to many insurance contexts, such as in the prediction of claim frequency or the estimation of pure risk premiums.
Secondly, we provide a systematic empirical comparison of some of the most popular approaches for modelling high-cardinality categorical features. Although the difficulty with such variables has been a long-standing issue in actuarial modelling, to the best of our knowledge, this paper is the first piece of work to extensively compare a range of the existing approaches (including our own GLMMNet). We believe that such a benchmarking study will be valuable, especially to practitioners facing many options available to them.
Importantly, the methodological extensions proposed in this paper, while motivated by actuarial applications, are not limited to this context and can have much wider applicability. The extension of the LMMNN to the GLMMNet would suit any applications where the response cannot be sufficiently modelled by a Gaussian distribution. This unparalleled flexibility with the form of the response distribution, in tandem with the computational efficiency of the algorithm (due to the fast and highly scalable variational inference used to implement GLMMNet), make the GLMMNet a promising tool for many practical ML applications.
In Section (ref), we introduce the GLMMNet, which fuses neural networks with GLMMs to enjoy both the predictive power of deep learning models and the statistical strength---specifically, the interpretability and likelihood-based estimation---of GLMMs. In Section (ref), we illustrate and compare the proposed GLMMNet against a range of alternative approaches across a spectrum of simulation environments, each to mimic some specific elements of insurance data commonly observed in practice. Section (ref) presents a real-life insurance case study and Section (ref) concludes.
The code used in the numerical experiments is available on \url{https://github.com/agi-lab/glmmnet}.
This section describes our GLMMNet in detail. We start by defining some notation in Section (ref). We then provide an overview of GLMMs in Section (ref) to set the stage for the introduction of the GLMMNet. Sections (ref)--(ref) showcase the architectural design of the GLMMNet and provide the implementation details.
We first introduce some notation to formalise the problem discussed above.
We write \(Y\) to denote the response variable, which is typically one of claim frequency, claim severity or risk premium for most insurance problems. We assume that the distribution of \(Y \in \mathcal{Y}\) depends on some covariates or features. In this work, we draw a distinction between standard features \(\mathbf{x} \in \mathcal{X}\) (i.e. numeric features and low-dimensional categorical features) and the high-cardinality categorical features \(\mathbf{z} \in \mathcal{Z}\) consisting of \(K\) features each with \(q_i~(i = 1, \cdots, K)\) categories. The latter is the subject of our interest. All categorical features have been one-hot encoded in the representations $\mathbf{x}$ and $\mathbf{z}$.
We assume \(\mathbf{x} \in \mathbb{R}^p\) and \(\mathbf{z} \in \mathbb{R}^q\), where \(p \in \mathbb{Z}\) and \(q = \sum_{i=1}^K q_i \in \mathbb{Z}\) represent the respective dimensions of the vectors. The empirical observations (i.e. data) are denoted \(\mathcal{D}=(y_i, \mathbf{x}_i, \mathbf{z}_i)_{i=1}^n\), where \(n\) is the number of observations in the dataset. For convenience, let us also write \(\mathbf{y}=[y_1, y_2, \cdots, y_n]^\top \in \mathbb{R}^n, \mathbf{X} = [\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_n]^\top \in \mathbb{R}^{n \times p}\) and \(\mathbf{Z} = [\mathbf{z}_1, \mathbf{z}_2, \cdots, \mathbf{z}_n]^\top \in \mathbb{R}^{n \times q}\). In general, we use bold font to denote vectors and matrices.
The purpose of any predictive model is to learn the conditional distribution \(p(y|\mathbf{x}, \mathbf{z})\). It is generally assumed that the dependence of the response \(Y\) on the covariates is through a true but unknown function \(\mu:(\mathcal{X}, \mathcal{Z}) \to \mathcal{Y}\), such that
where \(\boldsymbol{\phi}\) represents a vector of any additional or auxiliary parameters of the likelihood. Most commonly $\mu(\mathbf{x}, \mathbf{z})$ will be the conditional mean \(\mathbb{E}(Y|\mathbf{x}, \mathbf{z})\), although it can be interpreted more generally. The predictive model then approximates \(\mu(\cdot)\) by some parametric function \(\hat{\mu}(\cdot)\). Specifically, a linear regression model will assume a linear structure for \(\mu(\cdot)\) with \(\boldsymbol{\phi} = \sigma^2_\epsilon\) to account for the error variance. Indeed, for most ML models, the function $\mu(\cdot)$ is the sole subject of interest (as opposed to the entire distribution $p(y|\mathbf{x}, \mathbf{z})$). There are many different options for choosing the function $\mu(\cdot)$. We discuss and compare the main ones in this paper; see also Section (ref).
In what follows, for notational simplicity we will now assume that \(\mathbf{z}\) consists of only one (high-cardinality) categorical feature, i.e. \(K = 1\) and \(q = q_1\). The extension to multiple categorical features, however, is very straightforward; see also Section (ref).
Mixed models were originally introduced to model multiple correlated measurements on the same subject (e.g. longitudinal data), but they also offer an elegant solution to the modelling of high-cardinality categorical features frees2014. We now give a brief review of mixed models in the latter context. For further details, the books by mcculloch2001, gelman2007, denuit2019glm are all excellent references for this topic.
Let us start with the simplest linear model case. A one-hot encoded linear model assumes
where \(\mathbf{x}\) are standard features (e.g. sum insured, age), \(\mathbf{z}\) is a \(q\)-dimensional binary vector representation of a high-cardinality categorical variable (e.g. occupation), $\beta_0 \in \mathbb{R}$, $\boldsymbol{\beta} = [\beta_1, \beta_2, \cdots, \beta_p]^\top \in \mathbb{R}^p$, and \(\boldsymbol{\alpha} = [\alpha_1, \alpha_2, \cdots, \alpha_q]^\top \in \mathbb{R}^q\) with \(\sum_{i=1}^q \alpha_i = 0\) are the regression coefficients to be estimated. Note that the additional constraint \(\sum_{i=1}^q \alpha_i = 0\) is required to restore identifiability of parameters.
In high-cardinality settings, the vector \(\mathbf{z}\) becomes high-dimensional, and hence so must be \(\boldsymbol{\alpha}\). High cardinality therefore introduces a large number of parameters to the estimation problem. Moreover, in most insurance problems, the distribution of categories will be highly imbalanced, meaning that some categories will have much more exposure (i.e. data) than others. The rarely observed categories are likely to produce highly variable parameter estimates. As established in earlier sections, all this makes it extremely difficult to estimate the \(\boldsymbol{\alpha}\) accurately gelman2007.
Equation (ref) is also referred to as the no pooling model by antonio2014, as each $\alpha_i$, $i = 1, 2, \cdots, q$ has to be learned independently. At the other end of the spectrum, there is a complete pooling model, which simply ignores the \(\mathbf{z}\) and assumes
In the middle ground between the overly noisy estimates produced by (ref) and the over-simplistic estimates from (ref), we can find the (linear) mixed models, which assume
where \(\mathbf{u} = [u_1, \cdots, u_q]^\top\) with \(u_i \stackrel{iid}{\sim} \mathcal{N}(0, \sigma_u^2)\) are known as the random effects characterising the deviation of individual categories from the population mean, in contrast with the fixed effects \(\beta_0\) and \(\boldsymbol{\beta}\). Model (ref) is also known as a random intercept model.
The central idea in (ref) is that instead of assuming some fixed coefficient for each category, we assume that the effects of individual categories come from a distribution, so we only need to estimate (far fewer) parameters that govern the distribution of random effects rather than learning an independent parameter per category. This produces the equivalent effect of partial pooling gelman2007: categories with a smaller number of observations will have weaker influences on parameter estimates, and extreme values get regularised towards the collective mean (modelled by the fixed effects, i.e. \(\beta_0 + \mathbf{x}^\top \boldsymbol{\beta}\)).
In the same way that linear mixed models extend linear models, GLMMs extend GLMs by adding random effects capturing between-category variation to complement the fixed effects in the linear predictor.
Suppose that we have \(q\) categories in the data, each with \(n_i~(i = 1,2, \cdots, q)\) observations (so the total number of observations is \(n = \sum_{i=1}^q n_i\)). A GLMM is then defined by four components:
As explained earlier, the addition of random effects allows GLMMs to account for correlation within the same category, without overfitting to an individual category as is the case of a one-hot encoded GLM. The shared parameters that govern the distribution of random effects essentially enable information to transfer between categories, which is helpful for learning.
In Section (ref), we briefly reviewed previous attempts at embedding ML into mixed models. In particular, the LMMNN structure proposed by simchoni2022 is the closest to what we envision for a mixed effects neural network. However, it suffers from two major shortcomings that restrict its applicability to insurance contexts: First, the LMMNN assumes a Gaussian response with an identity link function, for which analytical results can be obtained, but which is ill suited to modelling skewed, heavier-tailed distributions that dominate insurance and financial data. Second, LMMNNs model the random effects in a non-linear way (through a sub-network in the structure). This complicates the interpretation of random effects, which carries practical significance in many applications.
The GLMMNet aims to address these concerns. The architectural skeleton of the model is depicted in Figure (ref). We adopt a very similar structure to that of simchoni2022, except that we remove the sub-network that they used to learn a non-linear function of \(\mathbf{Z}\). As noted earlier, the main purpose of our modification is to retain the interpretability of random effect predictions from the GLMM's linear predictor. In addition, we also want to avoid an over-complicated structure for the random effects, whose role is to act as a regulariser gelman2007.
In mathematical terms, we assume that the conditional \(\mathbf{y}|\mathbf{u}\) follows an ED family distribution, with
where \(g(\cdot)\) is the link function and \(f(\cdot)\) is learned through a neural network.
In the following, we give a detailed description of each component in the architecture plotted above, followed by a discussion on how the collective network can be trained in Section (ref). We remark that, in this work we have allowed \(f(\cdot)\) to be fully flexible to take full advantage of the predictive power of neural networks. We recognise that in some applications, interpretability of the fixed effects may also be desirable. In such cases, it is possible to replace the fixed effects module (described in Section (ref)) by a Combined Actuarial Neural Network wuethrich2019 structure.
The biggest limitation of GLMMs lies in their linear structure in (ref): similar to GLMs, features need to be hand-crafted to allow for non-linearities and interactions richman2021b, richman2021a. The proposed GLMMNet addresses this issue by utilising a multi-layer network component for the fixed effects, which is represented as the shaded structure in Figure (ref). For simplicity here we consider a fully connected feedforward neural network (FFNN), although many other network structures, e.g. convolutional neural networks (CNN) can also be easily accommodated.
A FFNN consists of multiple layers of neurons (represented by circles in the diagram) with non-linear activation functions to capture potentially complex relationships between input and output vectors; we refer to goodfellow2016 for more details. Formally, for a network with \(L\) hidden layers and \(q_l\) neurons in the \(l\)-th layer for \(l=1,\cdots,L\), the \(l\)-th layer is defined in Equation (ref) below.
with
where \(\varphi^{(l)}: \mathbb{R} \to \mathbb{R}\) is called the activation function for the \(l\)-th layer, \(b_j^{(l)} \in \mathbb{R}\) and \(\mathbf w_j^{(l)} \in \mathbb{R}^{q_{l-1}}\) respectively represent the bias term (or intercept in statistical terminology) and the network weights (or regression coefficients) for the \(j\)-th neuron in the \(l\)-th layer, and \(q_0\) is defined as the number of (preprocessed) covariates entering the network.
The network therefore provides a mapping from the (preprocessed) covariates \(\mathbf x\) to the desired output layer through a composition of hidden layers:
Training a network involves making many specific choices; Appendix (ref) gives more details.
Below we give a description of the random effects component of GLMMNet, which is the top (unshaded) structure in Figure (ref). For illustration, we stay with the case of a single high-cardinality categorical feature with \(q\) unique levels. The column of the categorical feature is first one-hot transformed into a binary matrix \(\mathbf{Z}\) of size \(n \times q\), which forms the input to the random effects component of the GLMMNet.
Here comes the key difference from the fixed effects network: The weights in the random effects layer, denoted by \(\mathbf{u}=[u_1, u_2, \cdots, u_q]^\top\), are not fixed values but are assumed to come from a distribution. This way, instead of having to estimate the effects of every individual category---many of which may lack sufficient data to support a reliable estimate, we only need to estimate the far fewer parameters that govern the distribution of $\mathbf{u}$ (in our case, just a single parameter \(\sigma_u^2\) in (ref)).
A less visible yet important distinction of our GLMMNet from the LMMNN simchoni2022 is that our model takes a Bayesian approach to the estimation of random effects (as opposed to an exact likelihood approach in the LMMNN). The Bayesian approach helps sidestep some difficult numerical approximations of the (marginal) likelihood function (as detailed in Section (ref)). We should point out that, although the Bayesian approach is popular for estimating GLMMs brms2017, there are further computational challenges that prevent them from being applied to ML mixed models; we elaborate on this in Section (ref).
Following the literature mcculloch2001, antonio2007, we assume that the random effects follow a normal distribution with zero mean, and that the random effects are independent across categories:
or equivalently, \(u_i \stackrel{iid}{\sim} \mathcal{N}(0, \sigma_u^2)\) for \(i = 1, 2, \cdots, q\). This is taken as the prior distribution of the random effects \(\mathbf{u}\).
In the context of our research problem, we are interested in the predictions for \(\mathbf{u}\), which indicate the deviation of individual categories from the population mean, or in other words, the excess risk they carry. Under the Bayesian framework, given the posterior distribution \(p(\mathbf{u} | \mathcal{D})\), we can simply take the posterior mode \(\hat{\mathbf{u}} = \mathop{\mathrm{argmax}}_\mathbf{u} p(\mathbf{u} | \mathcal{D})\) as a point estimate, which is also referred to as the maximum a posteriori (MAP) estimate. Alternatively, we can also derive interval predictions from the posterior distribution to determine whether the deviations as captured by \(\mathbf{u}\) carry statistical significance.
We should note that the estimation of the posterior distribution, which we treated as given in the above, represents a major challenge in Bayesian inference. We come back to this topic and discuss how we tackle it in Section (ref).
It should also be pointed out that the extension to multiple categorical features is very straightforward: all we need is to add another random effects layer (and one more variance parameter to estimate for each additional feature).
As explained earlier, one novel contribution of GLMMNet is its ability to accommodate many non-Gaussian distributions, notably gamma (commonly used for claim severity), (over-dispersed) Poisson (for claim counts), and Bernoulli (for binary classification tasks). Our extension to non-Gaussian distributions is an important pre-requisite for insurance applications.
Precisely, we assume that the responses \(Y_i, i = 1, 2, \cdots, n\), given the random effect \(u_{j[i]}\) for the category it belongs to (where $j[i]$ means that the $i$-th observation belongs to the $j$-th category), are conditionally independent with a distribution in the ED family. Equation (ref) implies that the conditional distribution of \(Y_{i} | u_{j[i]}\) is completely specified by the conditional mean \(\mu_{i}\) and the dispersion parameter \(\phi\), so we can write\footnote{The assumption of constant dispersion parameter $\phi$ is made to comply with the standard assumptions of GLM/GLMMs, but can be relaxed if desired by a slight modification. For the purpose of illustration in this work, we do not pursue this path down further but note that it can be an interesting direction for future research.} \[y_i | u_{j[i]} \sim \mathrm{ED}({\mu}_i, \phi).\]
In the GLMMNet, we assume that a link function \(g(\cdot)\) connects the conditional mean \(\mu_{i}\) with the non-linear predictor \(\eta_{i}\) formed by adding together the output from the fixed effects module and the random effects layer:
In implementation, we use the inverse link \(g^{-1}(\cdot)\) as the activation function for the final layer, to map the output \([f(\mathbf{x}_i), u_{j[i]}]\) from the fixed effects module and the random effects layer to a prediction for the conditional mean response \(\mu_i = g^{-1}(f(\mathbf{x}_i) + u_{j[i]}).\)
Finally, as GLMMNet operates under probabilistic assumptions, the dispersion parameter \(\phi\) (independent of covariates in our setting) can be estimated under maximum likelihood as part of the network. This is implemented by adding an input-independent but trainable weight to the final layer of GLMMNet keras, whose value will be updated by (stochastic) gradient descent on the loss function.
In a LMMNN simchoni2022, the network is trained by minimising the negative log-likelihood of \(\mathbf{y}\). As the LMMNN assumes a Gaussian density for \(\mathbf{y} | \mathbf{u}\) and \(\mathbf{u} \sim \mathcal{N}(\mathbf{0}, \boldsymbol{\Sigma})\), the marginal likelihood of \(\mathbf{y}\) can be derived analytically, i.e. \(\mathbf{y} \sim \mathcal{N}(f(\mathbf X), \mathbf Z \boldsymbol\Sigma \mathbf Z^\top + \sigma_\epsilon^2 \mathbf I)\). The network weights and biases, as well as the variance and covariance parameters \(\sigma_\epsilon^2\) and \(\sigma_u^2\) are learned by minimising the negative log-likelihood
where \(\mathbf{V} = \mathbf Z \boldsymbol\Sigma \mathbf Z^\top + \sigma_\epsilon^2 \mathbf I\) and \(\boldsymbol \Sigma := \operatorname{Cov}(\mathbf u) = \sigma_u^2 \mathbf{I}\).
In making the extension to GLMMNet, however, we no longer obtain a closed-form expression for the marginal likelihood:
where \(\boldsymbol{\beta}\) denote the weights and biases in the fixed effects component of the GLMMNet (Section (ref)), \(j[i]\) indicates the category to which the \(i\)-th observation belongs, \(\sigma_u^2\) is the variance parameter of the random effects, and \(\phi\) is the usual dispersion parameter for the ED density for the response.
We take a Bayesian approach to circumvent the difficult numerical approximations required to compute the integral in Equation (ref). Under the Bayesian framework, \(\pi(u_j) = f(u_j | \sigma_u) = \mathcal{N}(0, \sigma_u^2), ~j = 1, \cdots, q\), is taken as the prior for the random effects. Our goal is to make inference on the posterior of the random effects, given by
where \(\pi(\mathbf{u}) = \prod_{j=1}^q \pi(u_j)\) and $\mathcal{D}$ represents the data. Note that we again encounter an intractable integral in (ref). The traditional solution is to use Markov chain Monte Carlo (MCMC) to sample from the posterior without computing its exact form; however, the computational burden of MCMC techniques restricts its applicability to large and complex models like a deep neural network. This partly explains why Bayesian methods, although frequently used for estimating GLMMs brms2017, tend not to be as widely adopted among ML mixed models sigrist2022,hajjem2014.
With our GLMMNet, we apply variational inference to solve this problem. Variational inference is a popular approach among ML researchers working with Bayesian neural networks, due to its computational efficiency blei2017,zhang2019. Variational inference does not try to directly estimate the posterior, but proposes a surrogate parametric distribution \(q \in \mathcal{Q}\) to approximate it, therefore turning the difficult estimation into a highly efficient optimisation problem. In practice, the choice of the surrogate distribution is often made based on simplicity or computational feasibility. One particularly popular option is to use a mean-field distribution family blei2017, which assumes independence among all latent variables. In this work, we consider a diagonal Gaussian distribution for the surrogate posterior, \[q_{\boldsymbol{\vartheta}}(\mathbf{u})=\mathcal{N}({\boldsymbol{\mu}}_u, {\boldsymbol{\Sigma}}_u),\] where \(\boldsymbol{\mu}_u\) is a \(q\)-dimensional vector of the surrogate posterior mean of the \(q\) random effects, and \(\boldsymbol{\Sigma}_u = \mathrm{diag}(\sigma_1^2, \sigma_2^2, \cdots, \sigma_q^2)\) is a diagonal matrix of the posterior variance. More complex distributions, such as those that incorporate dependencies among latent variables, may provide a more accurate approximation. In this work, however, as we are mainly concerned with estimating the mean and dispersion of the posterior random effects, the chosen diagonal Gaussian surrogate distribution should suffice for this purpose.
The vector \(\boldsymbol{\mu}_u\) and the diagonal elements of \(\boldsymbol{\Sigma}_u\) together form the variational parameters of the diagonal Gaussian, denoted by \(\boldsymbol{\vartheta}\). The task here is to find the variational parameters \(\boldsymbol{\vartheta}^*\) that make the approximating density \(q_{\boldsymbol{\vartheta}} \in \mathcal{Q}\) closest to the true posterior, where closeness is defined in the sense of minimising the Kullback-Leibler (KL) divergence kullback1951. In mathematical terms, this means
We take the term inside the outermost bracket, which is called the evidence lower bound (ELBO) loss or the negative of variational free energy in neal1998, to be the loss function for the GLMMNet:
We see that the loss function used to optimise the GLMMNet is a sum of two parts: a first part that captures the divergence between the (surrogate) posterior and the prior, and a second part that captures the likelihood of observing the given data under the posterior. This decomposition echoes the trade-off between the prior beliefs and the data in any Bayesian analysis blei2017.
The loss function in (ref) can be calculated via Monte Carlo, that is, by generating posterior samples under the current estimates of the variational parameters $\boldsymbol{\vartheta}$, evaluating the respective expressions at each of the sampled values, and taking the average to approximate the expectation. The derivatives of the loss function can be computed via automatic differentiation keras and used by standard gradient descent algorithms for optimisation kingma2014.
Recall that \(\mathbf{x}\) denotes the standard input features and \(\mathbf{z} \in \mathbb{R}^q\) is a binary representation of the high-cardinality feature. Let \((\mathbf{x}^*, \mathbf{z}^*)\) represent a new data observation for which a prediction should be made. Note that when making predictions, it is possible to encounter a category that was never seen during training, in which case \(\mathbf{z}^*\) will be a zero vector (as it does not belong to any of the already-seen categories) and the predictor in Equation (ref) reduces to \(\eta^* = f(\mathbf{x}^*)\).
To make predictions from GLMMNet, we take expectations under the posterior distribution on the random effects blundell2015,jospin2022:
where \(q_{\boldsymbol{\vartheta}^{*}}(\cdot)\) denotes the optimised surrogate posterior. As a result of Equation (ref), any functional of \(y^* | \mathbf{x}^*, \mathbf{z}^*\), e.g. the expectation, can be computed in an analogous way.
The expectation in (ref) can be approximated by drawing Monte Carlo samples from the (surrogate) posterior, as outlined in Algorithm (ref). When there are multiple high-cardinality features, each will have its own binary vector representation and optimised surrogate posterior, and the expectation in Equation (ref) will be taken with respect to the joint posterior (which is simply a product of the surrogate posteriors when we assume independence between these features).
It is also worth pointing out that the model averaging in line (ref) has the equivalent effect of training an ensemble of networks, which is usually otherwise computationally impractical blundell2015.
In this section, we compare the performance of the GLMMNet against the most popular existing alternatives for handling high-cardinality categorical features. We list the candidate models in Section (ref) and the performance metrics we use to evaluate the models in Section (ref). Control over the true data generating process via simulation allows us to highlight the respective strengths of the models. We challenge the models under environments of varying levels of complexity, e.g. low signal-to-noise ratio, highly imbalanced distribution of categories, and/or skewed response distributions. The simulation environments detailed in Section (ref) are designed to mimic these typical characteristics of insurance data. Section (ref) summarises the results. Ensuing insights are used to inform the real data case study in Section (ref).
Listed below are the leading models widely used in practice that allow for high-cardinality categorical features. They will be implemented and used for comparison with the GLMMNet in this experiment.
Each model has its own features and properties, which renders the collection of models suitable for different contexts; Table (ref) gives a high-level summary of their different characteristics. In the table, “interpretability" refers to the ease with which a human can understand a model, as well as the ability to derive insights from the model. For instance, GLMs have high interpretability, because GLMs have a small number of parameters and each parameter carries a physical meaning, which contrasts with the case of a standard deep neural network characterised by a large number of parameters that are meaningless on their own. “Operational complexity" refers to the resources and time required to implement, train and fine-tune the model.
Based on this qualitative comparison only, the GLMMNet has good potential to improve on neural networks with entity embeddings, because it can similarly capture complex relationships while also offering a wider range of response distributions and better interpretability. However, further investigations are needed to understand and illustrate in which circumstances the GLMMNet approach can outperform other approaches, which is the focus of the rest of this section.
For all experiments presented in this paper (including the real data case study in Section (ref)), we randomly split the data into a training and a test set. The training set is used to select hyperparameters (by further splitting into an inner-training and a validation set) and fit the pool of candidate models. The different models are then evaluated and compared based on their performance on the test dataset. Specifically, to quantify the predictive performance, we consider the metrics listed below:
We generate \(n\) observations from a (generalised) non-linear mixed model with \(q\) categories (\(q \le n\)). We assume that the responses \(y_i, i = 1, \cdots, n\), each belonging to a category \(j[i] \in \{1, \cdots, q\}\), given the random effects \(\mathbf{u}=(u_j)_{j=1}^q\), are conditionally independent with a distribution in the ED family (e.g. gamma); denoted as \(y_i | u_{j[i]} \sim \mathrm{ED}({\mu}_i, \phi),\) where \({\mu}_i=\mathbb{E}(y_i | u_{j[i]})\) represents the mean and \(\phi\) represents a dispersion parameter shared across all observations. A link function \(g(\cdot)\) connects the conditional mean \({\mu}_i\) with a non-linear predictor in the form of
where \(\mathbf{x}_i\) is a vector of fixed effects covariates for the \(i\)-th observation, \(\mathbf{z}_i\) is a vector of random effects variables (in our case, a binary vector representation of a high-cardinality categorical feature), the random effects \({u}_{j} \stackrel{iid}{\sim} \mathcal{N}(0, \sigma_u^2)\) measure the deviation of individual categories from the population mean, and \(f(\cdot)\) is some complex non-linear function of the fixed effects.
Here, we consider
which was studied in friedman1991 and has become a standard simulation function for regression models (used in e.g. kuss2005; also available in the scikit-learn Python library by sklearn). In our application, we use the Friedman function to test the model's ability to filter out irrelevant predictors (note that \(x_6, \cdots, x_{10}\) is not used in \(f\)) as well as detect interactions and non-linearities of input features, all in the presence of a high-cardinality grouping variable that is subject to random effects.
We set up the desired simulation environments by adjusting the following parameters; refer to Appendix (ref) for details.
Table (ref) lists the simulation environments considered in our experiment. Each of the environments in experiments 2--5 has been configured to mimic one specific characteristic of practical insurance data, such as a low signal-to-noise ratio (the specification of $[8,1,4]$ was selected based on estimated parameters from the real insurance case study in Section (ref)), an imbalanced distribution across categories, or a highly skewed response distribution for claim severity. We give a brief description of each environment below.
For each scenario, we generate 5,000 training observations and 2,500 testing observations.
In all scenarios studied, the GLMMNet outperforms or at least performs on par with the target benchmark model of an entity embedded neural network. In low to medium noise level environments, the GLMMNet usually outperforms the second best model (i.e. the entity embedded neural network) by a significant margin. The two outperform all the other mixed models, including GLMM and GPBoost.
Figure (ref) compares the out-of-sample performance of the candidate models (listed in Section (ref)) over 50 simulation runs (of 5,000 training and 2,500 testing observations each) under a base scenario (experiment 1 in Table (ref)); the proposed GLMMNet is highlighted in green. For all metrics considered, smaller values mean a better fit.
From Figure (ref), we see that GLMMNet is leading in all metrics. In particular, on average it performs better than the current state-of-the-art neural network with entity embeddings (NN_ee), which is also the next best model. A paired sample Wilcoxon signed rank test on the out-of-sample MAE rejects the null hypothesis---that GLMMNet and NN_ee perform the same---with strong evidence (\(p < 0.001\)), in favour of the alternative that GLMMNet produces a smaller error. The same conclusion holds for all metrics tested. Indeed, upon closer examination of the paired sample data, we find that for every single simulation run, the performance of GLMMNet surpasses that of \texttt{NN_ee}. This suggests that optimising the networks via the (more suitable) mixed model likelihood provides a significant benefit to their predictive performance, in terms of both point predictions and (more notably) probabilistic predictions.
Ignoring GLMMNet, the best performing models in panels (a)--(c) are NN_ee and GBM_GLMM_enc (a boosting model with GLMM encoding). Both models significantly outperform the linear family of models (GLM or GLMM), showing that a flexible model structure is required to capture the non-linearities in the data.
More importantly, without a suitable way of accounting for the high-cardinality categorical feature, a flexible model structure alone is not enough to achieve good performance. For example, the struggle of tree-based models in dealing with categorical features has been well documented prokhorenkova2018. In this experiment, with the usual one-hot encoding, the GBM model (GBM_one_hot) performs even worse than its GLM counterpart (GLM_one_hot); although the reverse is true when both are fitted without the categorical variable altogether. This confirms the motivation for this research: more flexible ML models, when dealing with financial or insurance data, can often be constrained by their capability to model categorical variables with many levels.
The ranking of models in Figure (ref)(d)---which is designed to compare the models' ability to explain between-category variations---is slightly different. GLM_one_hot and GLMM perform better, moving up in the ranking to join the (consistently) top-performing GLMMNet and NN_ee. This is a result of the balance property (for GLMs with canonical links) which ensures that the sum of fitted values must equal exactly the sum of observed values for any level of the categorical features in the training set wuethrich2021. In this case, the training RMSE_avg will be zero by definition. In practice, it is remarkable that this quality often carries over to an out-of-sample set as well, as demonstrated by the performance of GLM_one_hot and GLMM here.
Figure (ref) gives a more insightful picture of per-category predictive accuracy by plotting the densities of predictions from selected models against the true underlying density that generated the observations (leftmost). The plot shows a selection of representative categories (top three with the highest response values, two in the middle, and bottom three with the lowest), but the conclusions drawn here hold generally for other categories too. The colours roughly correspond to the range of values under which the predictions fall; darker colours represent smaller values (relative to an average observation), and lighter represent larger values. Ideally, the predicted density shapes should match as closely as possible to the true densities, and the colours should also agree as much as possible.
In line with the previous observations, GLMMNet produces very accurate density estimates and models the variation in the data extremely well. NN_ee, on the other hand, clearly overestimates the response for Category 9, likely due to an overfit to the training data points. Unsurprisingly, ignoring the categorical variable (GLM_ignore_cat) almost eliminates the between-category differences and leads to very poor estimates. When equipped with one-hot encoding, the GLM (GLM_one_hot) is able to sense the differences in the mean (as indicated by the colours), as the balanced distribution of categories provides sufficient data points for learning the between-category differences. However, GLM_one_hot tends to produce narrower densities than the truth. This indicates a failure to capture the remaining, or within-category, variation in the data, confirming that its linear structure is too restrictive for this dataset.
The results from experiments 2--3 are mainly consistent with what we have seen in experiment 1 (see Appendix (ref)). GLMMNet takes a strong lead in both cases and showcases a statistically significant advantage over its closest competitor NN_ee in all cases (\(p < 0.01\) from the corresponding Wilcoxon signed rank tests).
The result in experiment 2 (gamma-distributed response) is much within expectation; the GLMMNet by design takes good care of the response distribution via its loss function. The result in experiment 3 (skewed distribution of categories) further suggests that the predictive strength of the GLMMNet is, at least to some degree, immune to the shape of the categorical distribution.
As the noise level increases (experiments 4--6), the differences in model performance become smaller, and the predictive advantage of the GLMMNet starts to fade away (see Appendix (ref)). In experiments 5--6, the GBM family of models surpasses both GLMMNet and NN_ee in CRPS measures. One possible explanation is that the variance of the random effects here is too small with respect to the error variance to warrant the need for accurately modelling the categorical variable. Although GLMMNet may not always provide a boost to the overall predictive performance, it remains one of the top performers in terms of its ability to capture the individual category means (as measured by RMSE_avg). The latter may be of practical significance in many applications.
The outperformance of GBM over deep learning models as observed in experiments 5 and 6 was already discussed in the prior literature, see, e.g. borisov2021, though it remains an open question what causes this gap in performance. It is possible that neural networks are more likely to overfit to the high noise in the background, or that they may get stuck in a suboptimal local minimum and cannot gain enough momentum from the data to escape. In general, the high level of noise in the environment prevents the network-based models (including the proposed GLMMNet) from fully realising their predictive power. Further investigations suggest that adding regularisation helps improve the overall performance of the GLMMNet to a great extent, as shown in panel (a) of Figure (ref), though it seems to bias the point predictions (panel b).
Contrasting experiments 5 and 6 shows that the GLMMNet ranks higher when the response moves away from Gaussian (see Figure (ref) in Appendix (ref)). The predictive power that the GLMMNet regains in experiment 6 is likely a result of the change in the true distribution of the response. When the shape of the response looks less normal, the outperformance of GBM models over the GLMMNet becomes less distinctive. This comparison highlights the advantage of the GLMMNet in using a loss function that is better aligned with the distribution of the response. We expect this effect to be more pronounced when we work with heavier-tailed distributions that deviate even more from the normal.
Across all experiments, we find that the GLMMNet, NN_ee, and GBM_GLMM_enc consistently outperform the other models in terms of predictive performance. In particular, we find that our GLMMNet outperforms or at least performs on par with the target benchmark model of an entity embedded neural network in all scenarios studied. That said, each model has its own strengths and limitations, and the choice of which model to use will depend on the specific needs of the modelling problem. Table (ref) lists some of these considerations. We should also note that GBM_ee has similar properties to GBM_GLMM_enc and is therefore not included in the table for brevity of presentation (except for the fact that GBM_ee is more complicated to implement as it requires training an entity-embedded neural network beforehand).
In this section, we apply GLMMNet to a real (proprietary) insurance dataset (described in Section (ref)), and discuss its performance relative to the other models we have considered (Section (ref)). We also demonstrate, in Section (ref), how such an analysis can potentially inform practitioners in the context of technical pricing.
The data for this analysis were provided by a major Australian insurer. The original data cover 27,351 commercial building and contents reported claims by small and medium-sized enterprises (SME) over the period between 2010 and 2015. The analysis seeks to construct a regression model to predict the ultimate costs of the reported claims based on other individual claim characteristics that can be observed early in the claims process (e.g. policy-level details and claim causes). A description of the (post-engineered) input features is given in Appendix (ref).
The response variable is the claim amount (claim severity), which has a very skewed distribution, as shown in Figure (ref): even on a log scale, we can still observe some degree of positive skewness. In search for a suitable distribution to model the severity, we fit both lognormal and loggamma distributions to the (unlogged) marginal. Figure (ref) indicates that both models may be suitable, with loggamma being slightly more advantageous. We will test both lognormal and loggamma distributions in our experiments below (Section (ref)); specifically, we fit Gaussian and gamma models to the log-transformed response.
The high-cardinality categorical variable of interest is the occupation of the business, which is coded through the Australian and New Zealand Standard Industrial Classification (ANZSIC) system. The ANZSIC system hierarchically classifies occupations into Divisions, Subdivisions, Groups and Classes (from broadest to finest). See Appendix (ref) for an example.
We will look at occupations at the Class level.\footnote{We chose not to incorporate the hierarchical information. Instead, we directly modelled the lowest level (i.e. Class) as a flat categorical feature. The decision was made in light of (1) the increased complexity of estimation in hierarchical models, and (2) the lack of alignment between the hierarchical system built for general-purpose occupation classification and the factors that differentiate claim severities.} There are over 300 unique levels of such occupation classes in this data set. Panel (a) of Figure (ref) highlights the skewed distribution of occupations: the most common occupation has more than double the number of observations than the second most common, and the number of observations decays rapidly for the rarer classes. Panel (b) shows the heterogeneity in claims experiences between different occupation classes; calculations reveal a coefficient of variation of 160% for the occupation means. One challenge in modelling the occupation variable lies in determining how much confidence we can have in the observed claims experiences; intuitively, we should trust the estimates more when there are more data points, and less when there are few data points. As explained in Section (ref), the mixed effects models are one solution to this problem, which we will explore further in this section.
The dataset under consideration displays many typical characteristics of insurance data that we studied in Section (ref) (see Section (ref) in particular). For example, it is characterised by an extremely skewed response distribution (Figure (ref)), a low signal-to-noise ratio (there is little information that covariates can add to the marginal fits, as implied by Figure (ref)) and a highly disproportionate distribution of categories (Figure (ref)). These features of the dataset best align with those of Experiment 6 in Section (ref). If the conclusions generalise, we expect that a regularised version of our GLMMNet will show strong performance.
We consider the same candidate models as before; see Section (ref). Having noted the underperformance of GLMMNet in the presence of high noise and the benefit of adding regularisation to it, we also consider an \(\ell_2\)-regularised GLMMNet (GLMMNet_l2) in the hope that the regularisation would help improve performance.
We follow the learning pipeline described at the start of Section (ref) and split the data into a 90% training set and a 10% test set. Table (ref) presents the out-of-sample metrics for the top performing lognormal and loggamma models respectively, compared against a one-hot encoded GLM as a baseline (full results in Appendix (ref)).\footnote{In Table (ref) and the relevant tables in Appendix (ref), the median absolute error (MedAE) and NLL are calculated on the original scale of the data (i.e. after the predictions are back-transformed). The CRPS is calculated on the log-transformed scale, because there is no closed-form formula for the loggamma CRPS. As a result of this, the CRPS as presented here will be more tolerant of wrong predictions for extreme observations than the NLL.} Note that we have selected slightly different metrics from the simulation experiments to better accommodate the heavier-tailed response. The RMSE, for example, is no longer a suitable measure as it penalises very heavily for predictions that are far from the observed values, which can happen very frequently when there are a few extreme observations in the data.
The three metrics tell consistent stories most of the time. The results for CRPS and NLL are almost perfectly aligned, which is not surprising given that both are measures of how far the observed values depart from the probabilistic predictions. In general, we consider the CRPS and NLL measures to be more reliable estimates of the goodness-of-fit of the models than the MedAE, which only takes into account the point predictions but not the uncertainty associated with them.
The results mostly match our expectation. Among the lognormal models, the regularised GLMMNet_l2 takes the lead, followed by GBM_GLMM_enc, NN_ee, GBM_ee and GLMM, whose performances are so close together that it is impossible to distinguish between the four. The GLM family of models performs poorly on this data, much worse than their mixed model counterpart, i.e. \texttt{GLMM}, confirming that the presence of the high-cardinality categorical variable interferes with their modelling capabilities. Among the loggamma models, although \texttt{NN_ee} outperforms \texttt{GLMMNet_l2} in CRPS, \texttt{GLMMNet_l2} remains the most competitive model in terms of NLL (on par with \texttt{GLMM}). One practical consideration to note is that when the environment is highly noisy and the underlying model is uncertain, such as here, the interpretability of the model becomes even more important. In such cases, models that are more interpretable are preferred, which in turn, highlights the advantages of \texttt{GLMMNet} over \texttt{NN_ee}.
An important observation is that adding the \(\ell_2\) regularisation term clearly helps GLMMNet navigate the data better. In the lognormal models, the regularised network represents a 1.5% improvement in the CRPS and NLL from the original GLMMNet, which suggests that the original GLMMNet overfits to the training data (as can be confirmed from its deterioration in the test performance). Indeed, the amelioration of the score proves statistically significant ($p < 0.001$) under the Diebold-Mariano test gneiting2014. The results for the loggamma models lead to similar conclusions.
Comparing the lognormal and loggamma models, we find that assuming a loggamma distribution for the response improves the model fit in almost all instances (except the GLM models) and across all three performance metrics. While the CRPS and NLL values are reasonably similar between the two tables, the median absolute error (MedAE) halves when moving from lognormal to loggamma. An intuitive explanation is that the lognormal models struggle more with fitting to the extreme claims and thus tend to over-predict the middle-ranged values in an attempt to boost the predictions. The loggamma models, on the other hand, are more comfortable with the extremes, as those can be captured reasonably well by the long tail. These results confirm the need to consider alternative distributions beyond the Gaussian---one major motivation for our proposed extension of the LMMNN to the GLMMNet. Even within the Gaussian framework, as demonstrated by the lognormal models presented here, the GLMMNet can still be useful.
Overall, as we have seen in the simulation experiments, the differences in model performance are relatively small due to the low signal-to-noise ratio. Only a limited number of risk factors were observed, and they are not necessarily the true drivers of the claims cost. In the experiments above, a simpler structure like GLMM sometimes performs comparably with the more complex and theoretically more capable GLMMNet. That said, the results clearly demonstrate that the GLMMNet is a promising model in terms of predictive performance. We expect to see better results with less noisy data or a greater volume of data, e.g. when more information becomes available as a policy progresses and claims develop, reducing the amount of noise in the system. This has been demonstrated in the simulation experiments (Section (ref)).
One main advantage of mixed effects models in dealing with high-cardinality categorical features is the transparency they offer on the effects of individual categories. This is particularly important as it is often the case that the high-cardinality categorical feature is itself a key risk factor. For example, in this data set, there is a lot of heterogeneity in claims experience across different occupations, but the lack of data for some makes it hard to determine how much trust can be given to the experience.
The GLMMNet can provide useful insights in this regard. Figure (ref) shows the posterior predictions for the effects of individual occupations (ordered by decreasing $z$-scores), plotted with 95% confidence intervals. Unsurprisingly, occupations with a larger number of observations generally have tighter intervals (i.e. lower standard deviations); see Figure (ref) of Appendix (ref). This prevents us from overtrusting extreme estimates, which may happen simply due to the small sample size. In the context of technical pricing, such an analysis can be helpful for identifying occupation classes that are statistically more or less risky. This is not possible with the often equally or less high-performing neural networks with entity embeddings. With the latter, we can analyse the embedding vectors to identify occupations that are “similar” to each other in the sense of producing similar output, but we cannot go any further.
High-cardinality categorical features are often encountered in real-world insurance applications. The high dimensionality creates a significant challenge for modelling. As discussed in Section (ref), the existing techniques prove inadequate in addressing this issue. Recent advancements in ML modeling hold promise for the development of more effective approaches.
In this paper, we took inspiration from the latest developments in the ML literature and proposed a novel model architecture called GLMMNet that targets high-cardinality categorical features in insurance data. The proposed GLMMNet fuses neural networks with GLMMs to enjoy both the predictive power of deep learning models and the transparency and statistical strength of GLMMs.
GLMMNet is an extension to the LMMNN proposed by simchoni2022. The LMMNN is limited to the case of a Gaussian response with an identity link function. GLMMNet, on the other hand, allows for the entire class of ED family distributions, which are better suited to the modelling of skewed distributions that we see in insurance and financial data.
Importantly, in making this extension, we have made several modifications to the original LMMNN architecture, including:
In our experiments, we found that the GLMMNet often performed better than or at least comparably with the benchmark model of an entity embedded neural network, which is the most popular approach among actuarial researchers working with neural networks. Although they can both be outperformed by boosting models in an environment of low signal-to-noise ratio, adding regularisation to the GLMMNet was sufficient to bring the model back to the top. The GLMMNet's outperformance over the entity embedded neural network is an achievement worth celebrating, as it comes with the additional benefits of transparency on the random effects as well as the ability to produce probabilistic predictions (e.g. uncertainty estimates).
To sum up, the proposed GLMMNet has at least the following advantages over the existing approaches:
All these make GLMMNet a worthy addition to the actuaries' modelling toolbox for handling high-cardinality categorical features.
This work was presented at the 2022 Australasian Actuarial Education and Research Symposium (AAERS) in November 2022 (Canberra, Australia). The authors are grateful for the constructive comments received from colleagues present at the event.
This research was supported under Australian Research Council's Discovery Project DP200101859 funding scheme. Melantha Wang acknowledges financial support from UNSW Australia Business School. The views expressed herein are those of the authors and are not necessarily those of the supporting organisations.
The authors declare that they have no conflicts of interest.
The code used in the numerical experiments, as well as the simulation datasets, is available on \url{https://github.com/agi-lab/glmmnet}.
\setcounter{section}{0} \setcounter{subsection}{0} \setcounter{table}{0} \setcounter{figure}{0} \gdef\thetable{\Alph{table}} \gdef\thefigure{\Alph{figure}} \gdef\thesection{\Alph{section}} \setcounter{section}{0} \counterwithin{figure}{section} \counterwithin{table}{section}
Below we outline a few practical considerations that arise from the implementation of GLMMNet.
The random effects layer of GLMMNet (Section (ref)) requires the specification of a prior distribution. We mentioned that we would use a Gaussian prior as per common practice with GLMMs, but we did not specify the exact parameters for it. With GLMMNet, it is possible to have a trainable prior whose parameters can be found by gradient descent with respect to the loss function in (ref); this approach of estimating the prior from the data is known as empirical Bayes casella1985. We experimented with both fixed and trainable priors in our implementation. We found that having a trainable prior leads to worse performance in general. blundell2015 also reached the same conclusion from their experiments with Bayesian neural networks and they speculated that the algorithm would be more tempted to update the prior parameters than the posterior when the prior were trainable. We choose to work with fixed priors.
As to what values to use for the fixed prior, we reference the Stan documentation by gelman2020, which recommends the use of weakly informative priors, i.e. priors that will give way to the likelihood in the presence of sufficient data but will dominate in the absence of data. The exact prior parameters to use will depend on the data and the task at hand.
We use the softplus bijector function, defined as \(\sigma(\lambda) = \log(1+\exp(\lambda)), ~\lambda \in \mathbb{R}\), to reparametrise the scale parameters (\(\sigma_j, ~j=1, \cdots, q\)) in the diagonal Gaussian distribution (i.e. the surrogate posterior). This ensures that the surrogate scale parameters always stay within their support (i.e. in the positive region).
Furthermore, we found that it is important to shrink the value of the \(\sigma_j\)'s at the start of the learning process to help guide the algorithm in the right direction. In our implementation, this is achieved by adding a constant multiplier (e.g. 0.01) to the parametrisation of \(\sigma_j\)'s. In the absence of such a constraint, the algorithm tends to return unreasonably large \(\sigma_j\) values, which significantly deteriorates the performance of GLMMNet. It appears that the unguided GLMMNet converges to strange local minima where the model attributes all the variation to the noise.
We note that, in theory, adding a constant multiplier as we did does not modify the solution space and thus should return the same results regardless. We were able to empirically verify that if we gave the model enough training time, it was able to figure out the right range of values for the parameters at the end. The tweak we made, however, helps the algorithm immediately converge to a better local minimum. One possible explanation is that by reducing the size of the gradients on the scale parameters, it helps the model focus on learning the more important location (mean parameters) of the posterior random effects.
The network weights and biases in GLMMNet are initialised by Tensorflow's default Glorot uniform initialiser glorot2010. The dispersion parameter for the ED family, which is also learned as part of GLMMNet, is initialised with a fixed estimate roughly calibrated to the data at hand. It does not seem to have a major impact on GLMMNet's performance, but we still consider it good practice to run some small scale experiments with this choice of initial value.
GLMM encoding can be regarded as an easier and more flexible alternative to the machine learning mixed models---e.g. GPBoost sigrist2021,sigrist2022 or GLMMNet, which often present a convoluted estimation procedure. While simple enough, pargent2022 found that this encoding scheme performed very well across a range of prediction tasks and a variety of datasets.
In the numerical experiments, we implemented a cross-validated version of GLMM encoding, as presented in Algorithm (ref). The encoded values \(\mathbf{z}'\) then take the place of the original categorical feature to enter any subsequent ML model, e.g. a gradient boosting machine.
Training a network involves making many specific choices. Below we briefly describe the specifications we use for the two network models, NN_ee and GLMMNet. To allow a fair comparison, where possible, we keep those choices consistent across both networks.
As per common practice, the data are split into training and testing sets. The training set is used for learning (i.e. fitting of the models) and the testing set is used for assessing and comparing the performance of the optimised candidate models.
In this work, we only perform a minimal search for hyperparameters, as hyperparameter tuning is not the main purpose of this study. We use the validation set approach. For models that involve hyperparameters, we split the training set into an inner-training set and a validation set. Models are fitted on the inner-training set, and hyperparameters are selected based on their performance on the validation set. This approach is less systematic than cross validation, but from experience, it usually yields reasonably similar results at a much more sustainable computational cost.
We set up the desired simulation environments in Section (ref) by adjusting the following parameters: