Adversarial Variational Bayes Methods for Tweedie Compound Poisson Mixed Models

TL;DR

This paper introduces an adversarial variational inference method for Bayesian Tweedie mixed models, effectively handling intractable likelihoods and hierarchical structures, with improved bias and predictive performance demonstrated on real and simulated data.

Contribution

It develops a novel adversarial variational inference approach for Tweedie mixed models, addressing intractability and hierarchical complexities with enhanced accuracy.

Findings

Smaller estimation bias on random effects compared to MCMC.

Achieves state-of-the-art predictive performance.

Provides richer variance function estimation.

Abstract

The Tweedie Compound Poisson-Gamma model is routinely used for modeling non-negative continuous data with a discrete probability mass at zero. Mixed models with random effects account for the covariance structure related to the grouping hierarchy in the data. An important application of Tweedie mixed models is pricing the insurance policies, e.g. car insurance. However, the intractable likelihood function, the unknown variance function, and the hierarchical structure of mixed effects have presented considerable challenges for drawing inferences on Tweedie. In this study, we tackle the Bayesian Tweedie mixed-effects models via variational inference approaches. In particular, we empower the posterior approximation by implicit models trained in an adversarial setting. To reduce the variance of gradients, we reparameterize random effects, and integrate out one local latent variable of Tweedie. We also employ a flexible hyper prior to ensure the richness of the approximation. Our method is evaluated on both simulated and real-world data. Results show that the proposed method has smaller estimation bias on the random effects compared to traditional inference methods including MCMC; it also achieves a state-of-the-art predictive performance, meanwhile offering a richer estimation of the variance function.

Figures

Tables

Table 1. Table 1 : The pairwise Gini index comparison with standard error based on 20 20 20 random splits

	GLM[29]	PQL [29]	Laplace [15]	AGQ[16]	MCMC[17]	TDBoost[30]	AVB
GLM (autoclaim)	$/$	$-{2.97}_{6.28}$	${1.75}_{5.68}$	${1.75}_{5.68}$	$-{15.02}_{7.06}$	${1.61}_{6.32}$	${9.84}_{5.80}$
PQL	${7.37}_{5.67}$	$/$	${7.50}_{6.26}$	${7.50}_{5.72}$	${6.73}_{5.95}$	${0.81}_{6.22}$	${9.02}_{6.07}$
Laplace	${2.10}_{4.52}$	$-{1.00}_{5.94}$	$/$	${8.84}_{5.36}$	${4.00}_{4.61}$	${21.45}_{4.84}$	${20.61}_{4.54}$
AGQ	${2.10}_{4.52}$	$-{1.00}_{5.94}$	${8.84}_{5.36}$	$/$	${4.00}_{4.61}$	${21.45}_{4.84}$	${20.61}_{4.54}$
MCMC	${14.75}_{6.80}$	$-{1.06}_{6.41}$	${3.12}_{5.99}$	${3.12}_{5.99}$	$/$	${7.82}_{5.83}$	${11.88}_{5.50}$
TDBoost	${17.52}_{4.80}$	${17.08}_{5.36}$	${19.30}_{5.19}$	${19.30}_{5.19}$	${11.61}_{4.58}$	$/$	${20.30}_{4.97}$
AVB	$-{0.17}_{4.70}$	${0.049}_{5.62}$	${3.41}_{4.94}$	${3.41}_{4.94}$	${0.86}_{4.62}$	${11.49}_{4.93}$	$/$
GLM (fineroot)	$/$	${23.18}_{9.30}$	${35.87}_{6.52}$	${35.87}_{6.52}$	$-{15.73}_{10.63}$	${35.71}_{6.52}$	${35.64}_{6.53}$
PQL	$-{21.61}_{8.34}$	$/$	${30.90}_{8.86}$	${30.90}_{8.86}$	$-{21.45}_{8.38}$	${24.19}_{8.36}$	${28.97}_{8.92}$
Laplace	$-{12.55}_{7.37}$	$-{14.72}_{8.81}$	$/$	${15.07}_{6.41}$	$-{12.55}_{7.37}$	${15.33}_{7.36}$	${10.61}_{7.20}$
AGQ	$-{12.55}_{7.37}$	$-{14.72}_{8.81}$	${15.07}_{6.41}$	$/$	$-{12.55}_{7.37}$	${15.33}_{7.36}$	${10.61}_{7.20}$
MCMC	${17.27}_{10.25}$	${22.53}_{9.31}$	${35.10}_{6.54}$	${35.10}_{6.54}$	$/$	${35.10}_{6.54}$	${34.87}_{6.55}$
TDBoost	${22.47}_{6.80}$	${8.50}_{9.09}$	${22.63}_{6.80}$	${22.63}_{6.80}$	${11.61}_{6.80}$	$/$	${22.39}_{6.80}$
AVB	$-{8.261}_{7.66}$	$-{10.886}_{8.98}$	${2.13}_{7.28}$	${2.13}_{7.28}$	$-{8.26}_{7.66}$	${11.00}_{7.74}$	$/$

Table 2. Table 2 : Estimation on 𝒫 𝒫 \mathcal{P} and σ b 2 superscript subscript 𝜎 𝑏 2 \sigma_{b}^{2} for both AutoClaim and FineRoot data with random effects.

	${𝒫}_{\text{auto}}$	${\sigma }_{\text{b,auto}}^2$	${𝒫}_{\text{root}}$	${\sigma }_{\text{b,root}}^2$
PQL	$/$	$8.94\cdot {10}^{-5}$	$/$	$1.0346\cdot {10}^{-3}$
LAPLACE	$1.3394$	$4.06\cdot {10}^{-4}$	$1.4202$	$6.401\cdot {10}^{-3}$
AGQ	$1.3394$	$4.027\cdot {10}^{-4}$	$1.4202$	$6.401\cdot {10}^{-3}$
MCMC	$1.3458$	$4.575\cdot {10}^{-3}$	$1.4272$	$3.490\cdot {10}^{-2}$
AVB	$1.3403$	$3.344\cdot {10}^{-3}$	$1.4237$	$2.120\cdot {10}^{-2}$

Equations

$$P(Y,N=n|\lambda ,\alpha ,\beta )=$$

$$+$$

$$\{\begin{array}{cc}\lambda =& \frac{{\mu }^{2-\mathcal{P}}}{\varphi (2-\mathcal{P})}\\ \alpha =& \frac{2-\mathcal{P}}{\mathcal{P}-1}\\ \beta =& \varphi (\mathcal{P}-1){\mu }^{\mathcal{P}-1}\end{array},\{\begin{array}{cc}\mu =& \lambda \alpha \beta \\ \mathcal{P}=& \frac{\alpha +2}{\alpha +1}\\ \varphi =& \frac{{\lambda }^{1-\mathcal{P}}(\alpha \beta {)}^{2-\mathcal{P}}}{2-\mathcal{P}}.\end{array}$$

$$\theta^{*}=\arg\max_{\theta}\operatorname{E}_{Q_{\theta}(\bm{z})}\bigg{[}-\log\frac{Q_{\theta}(\bm{z})}{P(\bm{z})}+\log P(D|\bm{z})\bigg{]}$$

$$\displaystyle\phi^{*}=\arg\min_{\phi}\bigg{[}-$$

$$-$$

$$P(y_i|\mathbf{w},\mathbf{b},\mathbf{\Sigma })$$

$$=P(\mathbf{b},\mathbf{\Sigma })\cdot \sum _{j=1}^{T}P(y_i|n_j,\mathbf{w},\mathbf{b},\mathbf{\Sigma })\cdot P(n_j|\mathbf{w},\mathbf{b},\mathbf{\Sigma })$$

Full text

**Index Terms— ** Tweedie model, variational inference, insurance policy pricing

1 Introduction

Tweedie models [1, 2] are special members in the exponential dispersion family; they specify a power-law relationship between the variance and the mean: $\text{Var}(Y)=\text{E}(Y)^{\mathcal{P}}$. For arbitrary positive $\mathcal{P}$, the index parameter of the variance function, outside the interval of $(0,1)$, Tweedie corresponds to a particular stable distribution, e.g., Gaussian ($\mathcal{P}=0$), Poisson ($\mathcal{P}=1$), Gamma ($\mathcal{P}=2$), Inverse Gaussian ($\mathcal{P}=3$). When $\mathcal{P}$ lies in the range of $(1,2)$, the Tweedie model is equivalent to Compound Poisson–Gamma Distribution [3], hereafter Tweedie for simplicity. Tweedie serves as a special Gamma mixture model, with the number of mixtures determined by a Poisson-distributed random variable, parameterized by$\{\lambda,\alpha,\beta\}$ and denoted as: $Y=\sum_{i=1}^{T}G_{i},T\sim\text{Poisson}(\lambda),G_{i}\overset{i.i.d}{\sim}\text{Gamma}(\alpha,\beta).$ Tweedie is heavily used for modeling non-negative heavy-tailed continuous data with a discrete probability mass at zero (see Fig. 1a). As a result, Tweedie gains its importance from multiple domains [4, 5], including actuarial science (aggregated loss/premium modeling), ecology (species biomass modeling), meteorology (precipitation modeling). On the other hand, in many field studies that require manual data collection, for example in insurance underwriting, the sampling heterogeneity from a hierarchy of groups/populations has to be considered. Mixed-effects models can represent the covariance structure related to the grouping hierarchy in the data by assigning common random effects to the observations that have the same level of a grouping variable; therefore, estimating the random effects is also an important component in Tweedie modeling.

Despite the importance of Tweedie mixed-effects models, the intractable likelihood function, the unknown variance index $\mathcal{P}$, and the hierarchical structure of mixed-effects (see Fig. 1b) hinder the inferences on Tweedie models. Unsurprisingly, there has been little work devoted to full-likelihood based inference on the Tweedie mixed model, let alone Bayesian treatments. In this work, we employ variational inference to solve Bayesian Tweedie mixed models. The goal is to introduce an accurate and efficient inference algorithm to estimate the posterior distribution of the fixed-effect parameters, the variance index, and the random effects.

2 Related Work

To date, most practice of Tweedie modeling are conducted within the quasi-likelihood (QL) framework [6] where the density function is approximated by the first and second order moments. QL does not come with the schema of estimating the variance function, even though the estimation of the index parameter plays an important role in hypothesis testing and uncertainty evaluation [7]. Extended quasi-likelihood (EQL) [8], together with the profile likelihood method [9] is proposed to fill the gap of index estimation for QL. Nonetheless, EQL is not capable of handling exact zero values.

It was not until recently that numerical approximations to the Tweedie density function appeared. Well-known approaches include series expansion method [10] and Fourier inversion method [11]. Yang et al. [12] and Chen et al. [13] applied the gradient boosted trees that try to incorporate non-linear interactions in estimating the mean value of the Tweedie model. Qian et al. [14] investigated the grouped elastic net and Lasso methods for the Tweedie model. Nonetheless, the dependence on density approximations also makes the above algorithms vulnerable to numerical issues. Both series expansion and Fourier inversion methods are noted to suffer from unbounded increment in the number of required terms.

When the Tweedie model is expended to incorporate the random effects, the profiled quasi-likelihood (PQL) method still suffers from the numerical restrictions. Likelihood-based inference methods require approximation of the integral over the distribution of the random effect. Advances in numerical approximation like Laplace approximation [15] and adaptive Gauss-Hermite quadrature (AGQ) [16] are then applied to solve the mixed model. However the complex density of the Tweedie model stills bothers the profiled likelihood methods.

The expectation of a latent variable for an intractable likelihood function can typically be estimated using the EM algorithm in an empirical Bayes approach or MCMC methods in full Bayesian approach. Simsekli et al. [3, 17] proposed solving the Tweedie mixed-effect models by explicitly exploiting the latent variable in the Bayesian formulation. Both MCMC and Monte Carlo EM methods have been applied. Zhang et al. [17] compared the latent variable approach with the density function approximations aforementioned, and found that MCMC are computationally demanding on the high dimensional dataset, and dedicated supervisions are needed in order to adjust the Monte Carlo error and the sampling variance [18].

3 Tweedie Mixed-Effects Models

Tweedie EDM $f_{Y}(y|\mu,\phi,\mathcal{P})$ equivalently describes the compound Poisson–Gamma distribution when $\mathcal{P}\in(1,2)$. Tweedie assumes Poisson arrivals of the events, and Gamma–distributed ``cost" for each individual event. Judging on whether the aggregated number of arrivals $N$ is zero, the joint density for each observation can be written as:

[TABLE]

where $d_{0}(\cdot)$ is the Dirac Delta function at zero. The connection of the parameters $\{\lambda,\alpha,\beta\}$ of Tweedie Compound Poisson model with the parameters of the general Tweedie EDM model, parameterized by the mean, the dispersion and the index parameter $\{\mu,\mathcal{P},\phi\}$ is denoted as:

[TABLE]

A mixed-effects model contains both fixed effects and random effects; graphical model is shown in Fig. 1b. We denote the mixed model as: $\kappa(\mu_{i})=f_{\bm{w}}(\bm{X}_{i})+\bm{U}_{i}^{\top}\bm{b}_{i},\bm{b}_{i}{\sim}\mathcal{N}(\bm{0,\Sigma_{b}})$ where $\bm{X}_{i}$ is the $i$-th row of the design matrix of fixed-effects covariates, $\bm{w}$ are parameters of the fixed-effects model which could represent linear function or deep neural network. $\bm{U}_{i}$ is the $i$-th row of the design matrix associated with random effects, $\bm{b}_{i}$ is the coefficients of the random effect which is usually assumed to follow a multivariate normal distribution with zero mean and covariance $\bm{\Sigma_{b}}$, $\kappa(\cdot)$ is the link function, and $\mu_{i}$ is the mean of the $i$-th response variable $Y$. In this work, we have considered the random effects on the intercept.

In the context of conducting Bayesian inference on Tweedie mixed models, we define 1) the observed data $D=(\bm{x}_{i},\bm{u}_{i},y_{i})_{i=1,\ldots,M}$. 2) global latent variables $\{\bm{w},\bm{\sigma}_{\bm{b}}\}$, we assume $\bm{\Sigma_{b}}$ is a diagonal matrix with its diagonal elements $\bm{\sigma}_{\bm{b}}$; 3) local latent variables $\{n_{i},\bm{b}_{i}\}$, indicating the number of arrivals and the random effect. The parameters of Tweedie is thus denoted by $(\lambda_{i},\alpha_{i},\beta_{i})=f_{\lambda,\alpha,\beta}(\bm{x}_{i},\bm{u}_{i}|\bm{w},\bm{\sigma}_{\bm{b}})$. The latent variables thus contain both local and global ones $\bm{z}=(\bm{w},n_{i},\bm{b}_{i},\bm{\sigma}_{\bm{b}})_{i=1,\ldots,M}$, and they are assigned with prior distribution $P(\bm{z})$. The joint log-likelihood is computed by summing over the number of observations $M$ by $\sum_{i=1}^{M}\log\left[P(y_{i}|n_{i},\lambda_{i},\alpha_{i},\beta_{i})\cdot P(n_{i}|\lambda_{i})\cdot P(\bm{b}_{i}|\bm{\sigma}_{\bm{b}})\right]$. The goal here is to find the posterior distribution of $P(\bm{z}|D)$ and make future predictions via $P(y_{\text{pred}}|D,\bm{x}_{\text{pred}},\bm{u}_{\text{pred}})=\int P(y_{\text{pred}}|\bm{z},\bm{x}_{\text{pred}},\bm{u}_{\text{pred}})P(\bm{z}|D)\operatorname{d}\bm{z}$.

4 Methods

Adversarial Variational Bayes. Variational Inference (VI) [19] approximates the posterior distribution, often complicated and intractable, by proposing a class of probability distributions $Q_{\theta}({\bm{z}})$ (so-called inference models), and then finding the best set of parameters $\theta$ by minimizing the KL divergence between the proposal and the true distribution, i.e., $\text{KL}(Q_{\theta}({\bm{z}})||P(\bm{z}|D))$. Minimizing the KL divergence is equivalent to maximizing the evidence of lower bound (ELBO) [20], expressed as Eq. 3. Optimizing the ELBO requires the gradient information of $\nabla_{\theta}\text{ELBO}$. In our experiments, we find that the model-free gradient estimation method – REINFORCE [21] fails to yield reasonable results due to the unacceptably-high variance issue [22], even equipped with baseline trick [23] or local expectation [24]. We also attribute the unsatisfactory results to the over-simplified proposal distribution in traditional VI. Here we try to solve these two issues by employing the implicit inference models with variance reduction tricks.

[TABLE]

AVB [25, 26] empowers the VI methods by using neural networks as the inference models; this allows more efficient and accurate approximations to the posterior distribution. Since a neural network is black-box, the inference model is implicit thus have no closed form expression, even though this does not bother to draw samples from it. To circumvent the issue of computing the gradient from implicit models, the mechanism of adversarial learning is introduced; an additional discriminative network $T_{\phi}(\bm{z})$ is used to model the first term in Eq. 3. By building a model to distinguish the latent variables that are sampled from the prior distribution $p(\bm{z})$ from those that are sampled from the inference network $Q_{\theta}(\bm{z})$, namely, optimizing the blow equation (where $\sigma(x)$ is the sigmoid function):

[TABLE]

the ratio is estimated as $\operatorname{E}_{Q_{\theta}(\bm{z})}[\log\frac{Q_{\theta}(\bm{z})}{P(\bm{z})}]=\operatorname{E}_{Q_{\theta}(\bm{z})}[T_{\phi^{*}}(\bm{z})]$. AVB considers optimizing Eq. 3 and Eq. 4 as a two-layer minimax game. We apply stochastic gradient descent alternately to find a Nash-equilibrium. Such Nash-equilibrium, if reached, is a global optimum of the ELBO [26].

Reparameterizable Random Effect. In mixed-effects models, the random effect is conventionally assumed to be $\bm{b}_{i}{\sim}\mathcal{N}(\bm{0,\Sigma_{b}})$. In fact, they are reparameterisable. As such, we model the random effects by the reparameterization trick [27]; $\bm{b}_{i}$ is now written as $\bm{b}(\bm{\epsilon}_{i})=\bm{0}+\bm{\sigma}_{\bm{b}}\odot\bm{\epsilon}_{i}$ where $\bm{\epsilon}_{i}\sim\mathcal{N}(\bm{0,\bm{I}})$. The $\bm{\sigma}_{\bm{b}}$ is a set of latent variables generated by the inference network and the random effects become deterministic given the auxiliary noise $\bm{\epsilon}_{i}$.

Theorem 1** (Reparameterizable Random Effect).**

Random effect is naturally reparameterizable. With the reparameterization trick, the random effect is no longer restricted to be normally distributed. For any ``location-scale" family distributions (Students't, Laplace, Elliptical, etc.)

Proof.

See the Section 2.4 in [27] ∎

Note that when computing the gradients, we no longer need to sample the random effect directly, instead, we can back-propagate the path-wise gradients which could dramatically reduce the variance [28].

Hyper Priors. The priors of the latent variables are fixed in traditional VI. Here we however parameterize the prior $P_{\psi}(\bm{w})$ by $\psi$ and make $\psi$ also trainable when optimizing Eq. 3. We refer $\psi$ as a kind of hyper prior to $\bm{w}$. The intuition is to not constraint the posterior approximation by one over-simplified prior. We would like the prior to be adjusted so as to make the posterior $Q_{\theta}(\bm{z})$ close to a set of prior distributions, or a self-adjustable prior; this could further ensure the expressiveness of $Q_{\theta}(\bm{z})$. We can also apply the same trick if the class of prior is reparameterizable.

Variance Reduction. We find that integrating out the latent variable $n_{i}$ by Monte Carlo in Eq. 1 gives significantly lower variance in computing the gradients. As $n_{i}$ is a Poisson generated sample, the variance will explode in the cases where $Y_{i}$ is zero but the sampled $n_{i}$ is positive, and $Y_{i}$ is positive but the sampled $n_{i}$ is zero. This also accounts for why the direct application of REINFORCE algorithm fails to work. In practice, we find limiting the number of Monte Carlo samples between $2-10$, dependent on the dataset, has the similar performance as summing over to larger number.

[TABLE]

5 Experiments and Results

We compare our method with six traditional inference methods on Tweedie. We evaluate our method on two public datasets of modeling the aggregate loss for auto insurance polices [31], and modeling the length density of fine roots in ecology [32]. We separate the dataset in $50\%/25\%/25\%$ for train/valid/test respectively. Considering the sparsity and right-skewness of the data, we use the ordered Lorenz curve and its corresponding Gini index [33, 34] as the metric. Assuming for the $i_{th}/N$ observations, $y_{i}$ is the ground truth, $p_{i}$ to be the results from the baseline predictive model, $\hat{y}_{i}$ to be the predictive results from the model. We sort all the samples by the relativity $r_{i}=\hat{y_{i}}/{p_{i}}$ in an increasing order, and then compute the empirical cumulative distribution as $(\hat{F}_{p}(r)=\frac{\sum_{i=1}^{n}p_{i}\cdot\mathbbm{1}(r_{i}\leqslant r)}{\sum_{i=1}^{n}p_{i}},\hat{F}_{y}(r)=\frac{\sum_{i=1}^{n}y_{i}\cdot\mathbbm{1}(r_{i}\leqslant r)}{\sum_{i=1}^{n}y_{i}})$. The plot of $(\hat{F}_{p}(r),\hat{F}_{y}(r))$ is an ordered Lorenz curve, and Gini index is the twice the area between the Lorenz curve and the $45^{\circ}$ line. The Gini index is not symmetrical if different baselines are used. In short, a model with a larger Gini index indicates its superior capability of separating the observations. In the insurance context, Gini index profiles the model's capability of distinguishing policy holders with different risk levels.

Gini comparison is presented in Table. 1. Here only fixed effects are considered. AVB method shows the state-of-the-art performance on both AutoClaim and FineRoot datasets, even comparatively better than the boosting tree method that is considered to have a strong performance. By inspecting the posterior estimation of $\mathcal{P}$ in Fig. 2, we find AVB shows a richer representation of the index parameter, with two modals at $1.20$ and $1.31$. As $\mathcal{P}$ uniquely characterizes a Tweedie distribution, compared with one single value that traditional profiled likelihood method offers, flexible choices of $\mathcal{P}$ enable customized descriptions for insurance policies that may have different underlying risk profiles. Also note that AVB uses a neural sampler that does not involve any rejection procedures; unlike MCMC, it holds the potentials for large–scale predictions on high dimensional dataset.

To compare the estimations on random effects, we add CARTYPE" in the AutoClaim data as the random effect with 6 groups and PLANT" in the FineRoot data with 8 groups. From table 2, we can see that all the algorithms estimate $\mathcal{P}$ in a reasonable similar range while the results from MCMC and AVB are closer. The estimation of $\sigma_{b}^{2}$ has a large variance. The estimation from AGQ and Laplace methods are around ten times smaller than MCMC and AVB results while PQL produces even smaller estimations. This is consistent with our finding in the simulation study that AGQ and Laplace methods tends to underestimate the random effect variance. The AVB and MCMC estimations are considered more accurate; as the estimated variance is comparatively smaller than MCMC, the AVB estimation can be believed to be more accurate.

6 Conclusions

We present Adversarial Variational Bayes in solving Bayesian Tweedie mixed models. To empower the posterior approximation, we employ implicit models and training in an adversarial setting. To reduce variance of the gradients, we sum over the local latent variable of Tweedie model that is not parameterizable. We also make the prior distribution trainable to guarantee the expressiveness of the posterior. Our method outperforms all traditional methods on both simulated data and real-world data; meanwhile, the estimations on the variance parameter and the random effects are competitively accurate and more stable. Tweedie model is one fundamental model that is widely used in insurance policy pricing. As far as we are concerned, this is the first work that leverages recent progresses of deep learning methods to tackle insurance problems.