Tweedie regression, or Poisson-Gamma regressions ?

Yesterday, I was chating with a young and enthousiastic actuary, who asked a nice (and classical) question: is it the same, or not to use a Tweedie regression, or two regressions (Poisson, and Gamma). For distributions, the two are equivalent, but when we have heterogeneity and explanatory variable, I really think that using all information, and running two regressions is much more interesting.

Homogeneous case

In the homogenous case, without any explanatory variable, the Tweedie distribution and compound Poisson-gamma distribution are equivalent representation (i.e., it is simply a reparametrization)

Consider a Tweedie distribution, with variance function power $p\in(1,2)$, mean $\mu$ and scale parameter $\phi$, then it is a compound Poisson model,

• $N\sim\mathcal{P}(\lambda)$ with $$\lambda=\displaystyle{\frac{\phi \mu^{2-p}}{2-p}}$$
• $Y_i\sim\mathcal{G}(\alpha,\beta)$ with $$\alpha=\displaystyle{-\frac{p-2}{p-1}}\text{~and~}\beta=\displaystyle{\frac{\phi \mu^{1-p}}{p-1}}$$

Conversely, consider a compound Poisson model $N\sim\mathcal{P}(\lambda)$ and $Y_i\sim\mathcal{G}(\alpha,\beta)$, then

• variance function power is $$p=\displaystyle{\frac{\alpha+2}{\alpha+1}}$$
• mean is $$\mu=\displaystyle{\frac{\lambda \alpha}{\beta}}$$
• scale (nuisance) parameter is
$$\phi=\displaystyle{\frac{[\lambda\alpha]^{\frac{\alpha+2}{\alpha+1}-1}\beta^{2-\frac{\alpha+2}{\alpha+1}}}{\alpha+1}}$$

So the two are equivalent…

Heterogeneous case

Now, in the context of regression$$N_i\sim\mathcal{P}(\lambda_i)\text{ with }\lambda_i=\exp[\boldsymbol{x}_i^\top\boldsymbol{\beta}_{\lambda}]$$
and$$Y_{j,i}\sim\mathcal{G}(\mu_i,\phi)\text{ with }\mu_i=\exp[\boldsymbol{x}_i^\top\boldsymbol{\beta}_{\mu}]$$
Then $S_i=Y_{1,i}+\cdots+Y_{N,i}$ has a Tweedie distribution

• variance function power is $$p=\displaystyle{\frac{\phi+2}{\phi+1}}$$
• mean is $\lambda_i \mu_i$
• scale parameter is$$\displaystyle{\frac{\lambda_i^{\frac{1}{\phi+1}-1}}{\mu_i^{\frac{\phi}{\phi+1}}}\left(\frac{\phi}{1+\phi}\right)}$$

There are $1+2\text{dim}(\boldsymbol{X})$ degrees of freedom here. And a Tweedie regression is

• variance function power is $p\in(1,2)$
• mean is $$\mu_i=\exp[\boldsymbol{x}_i^{\top}\boldsymbol{\beta}_{\text{Tweedie}}]$$
• scale parameter is $\phi$

There are now $2+\text{dim}(\boldsymbol{X})$ degrees of freedom.

In the actuarial terminology

• $N$ is the annual claim frequency
• $Y$ is the cost of single claims
• $S$ is the annual cost for a single insurance policy

As explained in our book, frequency and costs can be explained by different features, so that itself is a motifivation to consider two models. But consider the following simulated data

n = 1e4 a=2 set.seed(123) x = runif(n) etan = exp(-2+a*x) N = rpois(n,etan) dfn = data.frame(y=N,x=x) I=rep(1:n,N) etaz = exp(2-a*x[I]) Z = rgamma(sum(N),etaz,20) dfz = data.frame(y=Z,x=x[I]) S=tapply(Z,as.factor(I),sum) V=as.numeric(S[as.character(1:n)]) V[is.na(V)]=0 dfy = data.frame(y=V,x=x)

We can run two regressions, for the frequency, and for the costs

regn = glm(y~x, family=poisson(link="log"),data=dfn) regz = glm(y~x, family=Gamma(link="log"),data=dfz)

For the tweedie regression, let us find the optimal power parameter

library(statmod) library(tweedie) glmtw = function(t){ m = glm(y~x, family=tweedie(var.power = t, link.power = 0),data=dfy) d = NULL if(t == 1) d = 1 AICtweedie(m, dispersion = d) } vt = seq(1.01,1.99,length=251) vg = Vectorize(glmtw)(vt) plot(vt,vg,log="y",type="l") i=which.min(vg)

and consider the associated Tweedie regression.

regy = glm(y~x, family=tweedie(var.power = vt[i], link.power = 0),data=dfy)

For frequency, there is a clear increase of the average frequency with $x$ (and significant)

 summary(regy)

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -3.00822 0.04101 -73.356 <2e-16 *** x           -0.02226 0.07154  -0.311  0.756 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Tweedie family taken to be 0.6516459)

For the individual costs, there is a clear decline of the average cost with $x$ (and highly significant)

 summary(regn)

Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.01508 0.04135 -48.73 <2e-16 *** x            1.99036 0.05887  33.81 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Now, if we consider the average cost for the policy, we have

summary(regy)

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -3.00822 0.04101 -73.356 <2e-16 *** x           -0.02226 0.07154  -0.311  0.756 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Tweedie family taken to be 0.6516459)

I.e., the average annual cost for a single policy does not depend on $x$ (it is clearly not significant). As the product of the frequency and the average costs tells more or less the same story…

If the outcome, the price, is the same, one could agree that having here the two regressions is much more informative for risk management (if one wants to introduce deductibles for instance).

Model selection, AIC and Tweedie regression

Just some simple codes to illustrate some points we will discuss this week, for the last course on GLMs, before the final exam.  We have mentioned that the Gamma distribution belongs to the exponential, so we can run a regression, and compute the associated AIC,

> set.seed(123)
> test.data = rgamma(n=2000, scale=1, shape=1)
> m1 = glm( test.data~1, family=Gamma(link=log))
> AIC(m1)
[1] 3997.332


The Gamma distribution is also a special case of the Tweedie distribution, with power 2

> library(statmod)
> library(tweedie)
> m2 = glm( test.data~1, family=tweedie(link.power=0, var.power=2) )
> AIC(m2)
[1] NA


Unfortunately, we cannot compute the AIC, and we need a trick (with the appropriate R function).

> AICtweedie(m2)
[1] 3997.332


Of course, we can do the same with the Poisson distribution, which also belongs to the exponential family

> test.data = rpois(n=2000, lambda=1)
> m3 = glm( test.data~1, family=poisson(link=log))
> m4 = glm( test.data~1, family=tweedie(link.power=0, var.power=1) )
> AIC(m3)
[1] 5124.61


Here, we have a problem with the AICtweedie function

> AICtweedie(m4)
[1] Inf


because we need to specify the dispersion parameter

> AICtweedie(m4, dispersion=1)
[1] 5124.61


We can now check: we generate some Gamma sample, and fit various Tweedie distribution, changing simply the variance function (which is a power function)

> set.seed(123)
> test.data = rgamma(n=2000, scale=1, shape=1)
> glmtw = function(t){
+ m1 = glm( test.data~1, family=tweedie(link.power=0, var.power=t) )
+ d = NULL
+ if(t == 1) d = 1
+ AICtweedie(m1, dispersion = d)
+
+ }
> vt = seq(1,2.7,length=100)
> vg = Vectorize(glmtw)(vt)
> plot(vt,vg,log="y",type="l")


The minimum of the AIC is close to 2, corresponding to the Gamma distribution

We can also try with a Poisson

> set.seed(123)
> test.data = rpois(n=2000, lambda=1)
> glmtw = function(t){
+ m1 = glm( test.data~1, family=tweedie(link.power=0, var.power=t) )
+ d = NULL
+ if(t == 1) d = 1
+ AICtweedie(m1, dispersion = d)
+
+ }
> vt = seq(1,2,length=100)
> vg = Vectorize(glmtw)(vt)
> plot(vt,vg,log="y",type="l")


The minimum is now close to 1, corresponding to the Poisson distriubtion (the variance is equal to the average)

Let us now try some compound Poisson distribution,

> rcpd=function(n,lambda,shape,scale){
+ N=rpois(n,lambda)
+ X=rgamma(sum(N),shape=shape, scale=scale)
+ I=as.factor(rep(1:n,N))
+ S=tapply(X,I,sum)
+ V=as.numeric(S[as.character(1:n)])
+ V[is.na(V)]=0
+ return(V)}


Let us generate some compound Poisson random variables, with Poisson distribution with average 1, and with Gamma jumps, with average and variance 1,

> set.seed(123)
> test.data = rcpd(n=2000, 1,1,1)
> glmtw = function(t){
+ m1 = glm( test.data~1, family=tweedie(link.power=0, var.power=t) )
+ d = NULL
+ if(t == 1) d = 1
+ AICtweedie(m1, dispersion = d)
+ }
> vt = seq(1.1,1.9,length=100)
> vg = Vectorize(glmtw)(vt)
> plot(vt,vg,log="y",type="l")


The optimal value for the power function is here 1.5, based on the AIC (relationships between Tweedie parameters and the compound Poisson ones are given in the slides)

We can now play a little bit with the variance of the jumps: they still have aveage 1, but they now have a smaller variance

> set.seed(123)
> test.data = rcpd(n=2000, 1,3,1/3)
> vt = seq(1.05,1.95,length=100)
> vg = Vectorize(glmtw)(vt)
> plot(vt,vg,log="y",type="l")


The optimal power for the Tweedie is closer to one, closer to the Poison case

while if we increase the variance of the jumps

> set.seed(123)
> test.data = rcpd(n=2000, 1,1/3,3)
> vt = seq(1.05,1.95,length=100)
> vg = Vectorize(glmtw)(vt)
> plot(vt,vg,log="y",type="l")


the optimal power is higher, closer to the Gamma distribution.

Usually, when I give a course on GLMs, I try to insist on the fact that the link function is probably more important than the distribution. In order to illustrate, consider the following dataset, with 5 observations

x = c(1,2,3,4,5) y = c(1,2,4,2,6) base = data.frame(x,y)

Then consider several model, with various distributions, and either an identity link (and in that case $\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=\mathbf{x}^T\mathbf{\beta}$) or a log link function (so that $\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=e^{\mathbf{x}^T\mathbf{\beta}}$)

regNId = glm(y~x,family=gaussian(link="identity"),data=base) regNlog = glm(y~x,family=gaussian(link="log"),data=base) regPId = glm(y~x,family=poisson(link="identity"),data=base) regPlog = glm(y~x,family=poisson(link="log"),data=base) regGId = glm(y~x,family=Gamma(link="identity"),data=base) regGlog = glm(y~x,family=Gamma(link="log"),data=base) regIGId = glm(y~x,family=inverse.gaussian(link="identity"),data=base) regIGlog = glm(y~x,family=inverse.gaussian(link="log"),data=base

One can also consider some Tweedie distribution, to be even more general

library(statmod) regTwId = glm(y~x,family=tweedie(var.power=1.5,link.power=1),data=base) regTwlog = glm(y~x,family=tweedie(var.power=1.5,link.power=0),data=base)

Consider the prediction obtained in the first case, with the linear link function

library(RColorBrewer) darkcols = brewer.pal(8, "Dark2") plot(x,y,pch=19) abline(regNId,col=darkcols[1]) abline(regPId,col=darkcols[2]) abline(regGId,col=darkcols[3]) abline(regIGId,col=darkcols[4]) abline(regTwId,lty=2)

The predictions are very very close, aren’t they ? In the case of the exponential prediction, we obtain

plot(x,y,pch=19) u=seq(.8,5.2,by=.01) lines(u,predict(regNlog,newdata=data.frame(x=u),type="response"),col=darkcols[1]) lines(u,predict(regPlog,newdata=data.frame(x=u),type="response"),col=darkcols[2]) lines(u,predict(regGlog,newdata=data.frame(x=u),type="response"),col=darkcols[3]) lines(u,predict(regIGlog,newdata=data.frame(x=u),type="response"),col=darkcols[4]) lines(u,predict(regTwlog,newdata=data.frame(x=u),type="response"),lty=2)

We can actually look closer. For instance, in the linear case, consider the slope obtained with a Tweedie model (that will include all the parametric familes mentioned here, actually)

pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base))$coefficients[2,1:2] Vgamma = seq(-.5,3.5,by=.05) Vpente = Vectorize(pente)(Vgamma) plot(Vgamma,Vpente[1,],type="l",lwd=3,ylim=c(.965,1.03),xlab="power",ylab="slope") The slope here is always very very close to one ! Even more if we add a confidence interval plot(Vgamma,Vpente[1,]) lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2) lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2) Heuristically, for the Gamma regression, or the Inverse Gaussian one, because the variance is a power of the prediction, if the prediction is small (here on the left), the variance should be small. So, on the left of the graph, the error should be small with a higher power for the variance function. And that’s indeed what we observe here erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base),newdata=data.frame(x=1),type="response")-y[x==1] Verreur = Vectorize(erreur)(Vgamma) plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(-.1,.04),xlab="power",ylab="error") abline(h=0,lty=2) Of course, we can do the same with the exponential models pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base))$coefficients[2,1:2] Vpente = Vectorize(pente)(Vgamma) plot(Vgamma,Vpente[1,],type="l",lwd=3)

or, if we add the confidence bands, we obtain

plot(Vgamma,Vpente[1,],ylim=c(0,.8),type="l",lwd=3,xlab="power",ylab="slope") lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2) lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2)

So here also, the “slope” is rather similar… And if we look at the error we make on the left part of the graph, we obtain

erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base),newdata=data.frame(x=1),type="response")-y[x==1] Verreur = Vectorize(erreur)(Vgamma) plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(.001,.32),xlab="power",ylab="error")

So my point is that the distribution is usually not the most important point on GLMs, even if chapters of books on GLMs are distribution based… But as mentioned in an another post, if you consider a nonlinear transformation, like we have with GAMs, the story is more complicated…

Actuariat de l’Assurance Non-Vie #9

Cette semaine, nous avons fini les modèles de tarification, avec une extention sur les modèles globaux (sans passer par une séparation entre fréquence et coût moyen), et sur le choix de variables, et le choix de modèles. Les slides introductifs sont en ligne.

Actuariat de l’Assurance Non-Vie #9

Pour le neuvième chapitre du cours d’actuariat de l’assurance non-vie à l’ENSAE, un petit fourre-tout avant d’attaquer la modélisation du passif, en parlant un peu de modèles Tweedie (modèle collectif vs. modèles individuels), de choix de variables, et de choix de modèles. Les slides sont en ligne (la version pdf téléchargeable est comme souvent plus complète que celle sur slideshare)

“Improving Segmentation” (using Lorenz curves, or sort of)

This afternoon, André did send me an interesting graph about the use of Lorenz curve in the context of insurance pricing (and modeling)

It is some sort of Lorenz curve, upside-down, with on the x-axis the proportion of the population, and on the y-axis the proportion of the losses. The important point is that the population is sorted according the their risk, i.e. their premium. The code to generate such a curve is actually quite simple,

L <- function(u,varx="premium",vary="losses"){
base=base[order(base[,varx],decreasing=TRUE),]
base$cum=(1:nrow(base))/nrow(base) return(sum(base[base$cum<=u,vary])/
sum(base[,vary]))}

vu=seq(0,1,by=.01)
vv=Vectorize(function(u) L(u))(vu)

My concern was more on two labels on the figure, with on the top-left “perfect pricing” and on the first diagonal “average pricing“. What could that possibly mean? Is there even such a thing as a “perfect pricing“? In order to understand what we plot here, let us generate some dataset, and fit some model. Including things that might be seen as the “perfect model“: the price base on the parameters used to generate the data, and the model used to generate the data, fitted on the data.