```
import arviz as az
import bambi as bmb
import matplotlib.pyplot as plt
import numpy as np
```

# Wald and Gamma Regression (Australian insurance claims 2004-2005)

```
"arviz-darkgrid")
az.style.use(1234) np.random.seed(
```

## Load and examine Vehicle insurance data

In this notebook we use a data set consisting of 67856 insurance policies and 4624 (6.8%) claims in Australia between 2004 and 2005. The original source of this dataset is the book Generalized Linear Models for Insurance Data by Piet de Jong and Gillian Z. Heller.

```
= bmb.load_data("carclaims")
data data.head()
```

veh_value | exposure | clm | numclaims | claimcst0 | veh_body | veh_age | gender | area | agecat | |
---|---|---|---|---|---|---|---|---|---|---|

0 | 1.06 | 0.303901 | 0 | 0 | 0.0 | HBACK | 3 | F | C | 2 |

1 | 1.03 | 0.648871 | 0 | 0 | 0.0 | HBACK | 2 | F | A | 4 |

2 | 3.26 | 0.569473 | 0 | 0 | 0.0 | UTE | 2 | F | E | 2 |

3 | 4.14 | 0.317591 | 0 | 0 | 0.0 | STNWG | 2 | F | D | 2 |

4 | 0.72 | 0.648871 | 0 | 0 | 0.0 | HBACK | 4 | F | C | 2 |

Let’s see the meaning of the variables before creating any plot or fitting any model.

**veh_value**: Vehicle value, ranges from \$0 to \$350,000.**exposure**: Proportion of the year where the policy was exposed. In practice each policy is not exposed for the full year. Some policies come into force partly into the year while others are canceled before the year’s end.**clm**: Claim occurrence. 0 (no), 1 (yes).**numclaims**: Number of claims.**claimcst0**: Claim amount. 0 if no claim. Ranges from \$200 to \$55922.**veh_body**: Vehicle body type. Can be one of bus, convertible, coupe, hatchback, hardtop, motorized caravan/combi, minibus, panel van, roadster, sedan, station wagon, truck, and utility.**veh_age**: Vehicle age. 1 (new), 2, 3, and 4.**gender**: Gender of the driver. M (Male) and F (Female).**area**: Driver’s area of residence. Can be one of A, B, C, D, E, and F.**agecat**: Driver’s age category. 1 (youngest), 2, 3, 4, 5, and 6.

The variable of interest is the claim amount, given by `"claimcst0"`

. We keep the records where there is a claim, so claim amount is greater than 0.

`= data[data["claimcst0"] > 0] data `

For clarity, we only show those claims amounts below \$15,000, since there are only 65 records above that threshold.

`"claimcst0"] > 15000].shape[0] data[data[`

`65`

```
"claimcst0"] <= 15000]["claimcst0"], bins=30)
plt.hist(data[data["Distribution of claim amount")
plt.title("Claim amount ($)"); plt.xlabel(
```

And this is when you say: “Oh, there really are ugly right-skewed distributions out there!”. Well, yes, we’ve all been there :)

In this case we are going to fit GLMs with a right-skewed distribution for the random component. This time we will be using Wald and Gamma distributions. One of their differences is that the variance is proportional to the cubic mean in the case of the Wald distribution, and proportional to the squared mean in the case of the Gamma distribution.

## Wald family

The Wald family (a.k.a inverse Gaussian model) states that

\[ \begin{array}{cc} y_i \sim \text{Wald}(\mu_i, \lambda) & g(\mu_i) = \mathbf{x}_i^T\beta \end{array} \]

where the pdf of a Wald distribution is given by

\[ f(x|\mu, \lambda) = \left(\frac{\lambda}{2\pi}\right)^{1/2}x^{-3/2}\exp\left\{ -\frac{\lambda}{2x} \left(\frac{x - \mu}{\mu} \right)^2 \right\} \]

for \(x > 0\), mean \(\mu > 0\) and \(\lambda > 0\) is the shape parameter. The variance is given by \(\sigma^2 = \mu^3/\lambda\). The canonical link is \(g(\mu_i) = \mu_i^{-2}\), but \(g(\mu_i) = \log(\mu_i)\) is usually preferred, and it is what we use here.

## Gamma family

The default parametrization of the Gamma density function is

\[ \displaystyle f(x | \alpha, \beta) = \frac{\beta^\alpha x^{\alpha -1} e^{-\beta x}}{\Gamma(\alpha)} \]

where \(x > 0\), and \(\alpha > 0\) and \(\beta > 0\) are the shape and rate parameters, respectively.

But GLMs model the mean of the function, so we need to use an alternative parametrization where

\[ \begin{array}{ccc} \displaystyle \mu = \frac{\alpha}{\beta} & \text{and} & \displaystyle \sigma^2 = \frac{\alpha}{\beta^2} \end{array} \]

and thus we have

\[ \begin{array}{cccc} y_i \sim \text{Gamma}(\mu_i, \sigma_i), & g(\mu_i) = \mathbf{x}_i^T\beta, & \text{and} & \sigma_i = \mu_i^2/\alpha \end{array} \]

where \(\alpha\) is the shape parameter in the original parametrization of the gamma pdf. The canonical link is \(g(\mu_i) = \mu_i^{-1}\), but here we use \(g(\mu_i) = \log(\mu_i)\) again.

## Model fit

In this example we are going to use the binned age, the gender, and the area of residence to predict the amount of the claim, conditional on the existence of the claim because we are only working with observations where there is a claim.

`"agecat"`

is interpreted as a numeric variable in our data frame, but we know it is categorical, and we wouldn’t be happy if our model takes it as if it was numeric, would we?

We have two alternatives to tell Bambi that this numeric variable must be treated as categorical. The first one is to wrap the name of the variable with `C()`

, and the other is to pass the same name to the categorical argument when we create the model. We are going to use the first approach with the Wald family and the second with the Gamma.

The `C()`

notation is taken from Patsy and is encouraged when you want to explicitly pass the order of the levels of the variables. If you are happy with the default order, better pass the name to `categorical`

so tables and plots have prettier labels :)

### Wald

```
= bmb.Model("claimcst0 ~ C(agecat) + gender + area", data, family = "wald", link = "log")
model_wald = model_wald.fit(tune=2000, target_accept=0.9, idata_kwargs={"log_likelihood": True}) fitted_wald
```

```
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [lam, Intercept, C(agecat), gender, area]
```

```
Sampling 2 chains for 2_000 tune and 1_000 draw iterations (4_000 + 2_000 draws total) took 23 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
```

`; az.plot_trace(fitted_wald)`

` az.summary(fitted_wald)`

mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|

C(agecat)[2] | -0.164 | 0.102 | -0.345 | 0.030 | 0.004 | 0.003 | 765.0 | 822.0 | 1.01 |

C(agecat)[3] | -0.258 | 0.102 | -0.446 | -0.061 | 0.004 | 0.003 | 775.0 | 838.0 | 1.01 |

C(agecat)[4] | -0.263 | 0.101 | -0.442 | -0.063 | 0.004 | 0.003 | 803.0 | 899.0 | 1.01 |

C(agecat)[5] | -0.377 | 0.109 | -0.575 | -0.166 | 0.004 | 0.003 | 797.0 | 776.0 | 1.01 |

C(agecat)[6] | -0.320 | 0.122 | -0.550 | -0.086 | 0.004 | 0.003 | 888.0 | 1005.0 | 1.01 |

Intercept | 7.724 | 0.099 | 7.529 | 7.895 | 0.004 | 0.003 | 539.0 | 614.0 | 1.01 |

area[B] | -0.030 | 0.073 | -0.170 | 0.103 | 0.002 | 0.002 | 1056.0 | 1032.0 | 1.00 |

area[C] | 0.070 | 0.067 | -0.052 | 0.194 | 0.002 | 0.001 | 1244.0 | 1324.0 | 1.00 |

area[D] | -0.022 | 0.089 | -0.189 | 0.148 | 0.002 | 0.002 | 1305.0 | 1630.0 | 1.00 |

area[E] | 0.147 | 0.103 | -0.056 | 0.330 | 0.003 | 0.002 | 1347.0 | 1370.0 | 1.00 |

area[F] | 0.368 | 0.137 | 0.124 | 0.641 | 0.003 | 0.003 | 1621.0 | 1347.0 | 1.00 |

gender[M] | 0.154 | 0.050 | 0.053 | 0.244 | 0.001 | 0.001 | 2273.0 | 1414.0 | 1.00 |

lam | 722.857 | 14.623 | 694.186 | 750.000 | 0.337 | 0.238 | 1890.0 | 1335.0 | 1.00 |

If we look at the `agecat`

variable, we can see the log mean of the claim amount tends to decrease when the age of the person increases, with the exception of the last category where we can see a slight increase in the mean of the coefficient (-0.307 vs -0.365 of the previous category). However, these differences only represent a slight tendency because of the large overlap between the marginal posteriors for these coefficients (see overlaid density plots for `C(agecat)`

.

The posterior for `gender`

tells us that the claim amount tends to be larger for males than for females, with the mean being 0.153 and the credible interval ranging from 0.054 to 0.246.

Finally, from the marginal posteriors for the areas, we can see that F is the only area that clearly stands out, with a higher mean claim amount than in the rest. Area E may also have a higher claim amount, but this difference with the other areas is not as evident as it happens with F.

### Gamma

```
= bmb.Model(
model_gamma "claimcst0 ~ agecat + gender + area",
data,="gamma",
family="log",
link="agecat",
categorical
)= model_gamma.fit(tune=2000, target_accept=0.9, idata_kwargs={"log_likelihood": True}) fitted_gamma
```

```
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [alpha, Intercept, agecat, gender, area]
```

```
Sampling 2 chains for 2_000 tune and 1_000 draw iterations (4_000 + 2_000 draws total) took 29 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
```

`; az.plot_trace(fitted_gamma)`

` az.summary(fitted_gamma)`

mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|

Intercept | 7.718 | 0.062 | 7.588 | 7.823 | 0.002 | 0.001 | 959.0 | 1220.0 | 1.0 |

agecat[2] | -0.182 | 0.063 | -0.299 | -0.060 | 0.002 | 0.002 | 829.0 | 1106.0 | 1.0 |

agecat[3] | -0.276 | 0.062 | -0.389 | -0.156 | 0.002 | 0.002 | 803.0 | 1207.0 | 1.0 |

agecat[4] | -0.268 | 0.062 | -0.385 | -0.152 | 0.002 | 0.001 | 874.0 | 1177.0 | 1.0 |

agecat[5] | -0.388 | 0.070 | -0.510 | -0.245 | 0.002 | 0.002 | 999.0 | 1434.0 | 1.0 |

agecat[6] | -0.314 | 0.078 | -0.461 | -0.171 | 0.002 | 0.002 | 1125.0 | 1223.0 | 1.0 |

alpha | 0.762 | 0.014 | 0.737 | 0.788 | 0.000 | 0.000 | 2389.0 | 1185.0 | 1.0 |

area[B] | -0.025 | 0.051 | -0.120 | 0.069 | 0.001 | 0.001 | 1706.0 | 1278.0 | 1.0 |

area[C] | 0.070 | 0.046 | -0.022 | 0.152 | 0.001 | 0.001 | 1498.0 | 1300.0 | 1.0 |

area[D] | -0.017 | 0.062 | -0.125 | 0.100 | 0.002 | 0.001 | 1358.0 | 1220.0 | 1.0 |

area[E] | 0.150 | 0.068 | 0.031 | 0.286 | 0.002 | 0.001 | 2016.0 | 1432.0 | 1.0 |

area[F] | 0.369 | 0.075 | 0.229 | 0.508 | 0.002 | 0.001 | 1984.0 | 1538.0 | 1.0 |

gender[M] | 0.166 | 0.034 | 0.105 | 0.229 | 0.001 | 0.000 | 2309.0 | 1276.0 | 1.0 |

The interpretation of the parameter posteriors is very similar to what we’ve done for the Wald family. The only difference is that some differences, such as the ones for the area posteriors, are a little more exacerbated here.

## Model comparison

We can perform a Bayesian model comparison very easily with `az.compare()`

. Here we pass a dictionary with the InferenceData objects that `Model.fit()`

returned and `az.compare()`

returns a data frame that is ordered from best to worst according to the criteria used.

```
= {"wald": fitted_wald, "gamma": fitted_gamma}
models = az.compare(models)
df_compare df_compare
```

rank | elpd_loo | p_loo | elpd_diff | weight | se | dse | warning | scale | |
---|---|---|---|---|---|---|---|---|---|

wald | 0 | -38581.390685 | 12.850506 | 0.000000 | 1.0 | 106.056885 | 0.000000 | False | log |

gamma | 1 | -39628.917221 | 26.536036 | 1047.526536 | 0.0 | 104.988496 | 35.778307 | False | log |

`=False); az.plot_compare(df_compare, insample_dev`

By default, ArviZ uses *loo*, which is an estimation of leave one out cross-validation. Another option is the widely applicable information criterion (WAIC). Since the results are in the log scale, the better out-of-sample predictive fit is given by the model with the highest value, which is the Wald model.

```
%load_ext watermark
%watermark -n -u -v -iv -w
```

```
Last updated: Sun May 26 2024
Python implementation: CPython
Python version : 3.11.9
IPython version : 8.24.0
bambi : 0.13.1.dev39+gb7d6a6cb
matplotlib: 3.8.4
numpy : 1.26.4
arviz : 0.18.0
Watermark: 2.4.3
```