I Fought the (distribution) Law (and the Law did not win)

A few days ago, I was asked if we should spend a lot of time to choose the distribution we use, in GLMs, for (actuarial) ratemaking. On that topic, I usually claim that the family is not the most important parameter in the regression model. Consider the following dataset

> db <- data.frame(x=c(1,2,3,4,5),y=c(1,2,4,2,6))
> plot(db,xlim=c(0,6),ylim=c(-1,8),pch=19)

To visualize a regression model, use the following code

> nd=data.frame(x=seq(0,6,by=.1))
> add_predict = function(reg){
+ prd1=predict(reg,newdata=nd,se.fit = TRUE,type="response")
+ y1=prd1$fit
+ y1_upp=prd1$fit+prd1$residual.scale*1.96*
prd1$se.fit   
+ y1_low=prd1$fit-prd1$residual.scale*1.96*
prd1$se.fit 
+ polygon(c(nd$x,rev(nd$x)),c(y1_upp,
rev(y1_low)),col="light green",angle=90,
density=40,border=NA)
+ lines(nd$x,y1,col="red",lwd=2)
+ }

For instance, with a Poisson regression (with a log link function) we get

> plot(db)
> reg1=glm(y~x,family=poisson(link="log"),
+ data=db)
> add_predict(reg1)

while, with a Gaussian regresion (but still with a log link function), we get

> plot(db)
> reg2=glm(y~x,family=gaussian(link="log"),
+ data=db)
> add_predict(reg2)

If we just care about the expected value of our prediction, the output is more or less the same

> plot(db)
> lines(nd$x,predict(reg1,newdata=nd,
+ type="response"),col="red",lwd=1.5)
> lines(nd$x,predict(reg2,newdata=nd,
+ type="response"),col="blue",lwd=1.5)

So, indeed, forget about the (distribution) law when running a GLM. Not convinced? Consider – on the same dataset – a Poisson regression (with an identity link function this time)

> plot(db)
> reg1=glm(y~x,family=poisson(link="identity"),
+ data=db)
> add_predict(reg1)

while, with a Gaussian regresion (but still with an identity link function), we get

> plot(db)
> reg2=glm(y~x,family=gaussian(link="identity"),
+ data=db)
> add_predict(reg2)

Again, if we just plot the expected value of our prediction, the output is more or less the same

> plot(db)
> lines(nd$x,predict(reg1,newdata=nd,
+ type="response"),col="red",lwd=1.5)
> lines(nd$x,predict(reg2,newdata=nd,
+ type="response"),col="blue",lwd=1.5)

So clearly, the simplistic message you should not care too much about the (distribution) law seems to be valid…

But if we compare those graphs (and the confidence intervals) it looks like the regression can be sensitive to outliers. Or not.

To visualize the impact of a change in the value of some observations, use

> change_pred_poisson=function(x){
+   prd=function(y){
+   db_rev=db
+   db_rev$y[3]=y
+   reg1=glm(y~x,family=poisson(link=log),
+ data=db_rev)
+   prd1=predict(reg1,newdata=data.frame(x=x),
+ se.fit = TRUE,type="response")
+   y1=prd1$fit
+   y1_upp=prd1$fit+prd1$residual.scale*
+ 1.96*prd1$se.fit   
+   y1_low=prd1$fit-prd1$residual.scale*1.96*
+ prd1$se.fit 
+   c(y1,y1_low,y1_upp)
+   }
+   vx=seq(1,8,by=.1)
+   s=sapply(vx,prd)
+   plot(vx,vx,col="white",xlab="",ylab="")
+   polygon(c(vx,rev(vx)),c(s[2,],rev(s[3,])),
+ col="light green",border=NA)
+   lines(vx,s[1,],col="red",lwd=2)
+   invisible(s)
+ }

Here we change the third observation. Its https://latex.codecogs.com/gif.latex?y-value was 4. Now we change it, from 1 up to 8.  And we look at the predicted value, with a Poisson regression, for the same value, the third one

> s_poisson=change_pred_poisson(3)

Now, we can write a similar code,  with a Gaussian regression

> s_gaussian=change_pred_gaussian(3)

Even with a Gamma regression

> s_gamma=change_pred_gamma(3)

If we compare now the three models, we get

> plot(seq(1,8,by=.1),s_gauss[1,],
+ type="l",lwd=2,xlab="",ylab="")
> lines(seq(1,8,by=.1),s_poisson[1,],
+ col="red",lwd=2)
> lines(seq(1,8,by=.1),s_gamma[1,],col=
+ "blue",lwd=2)
> legend("topleft",c("Gaussian","Poisson",
+ "Gamma"),col=c("black","red","blue"),lty=1)

So, the (distribution) law has an impact on our prediction. When https://latex.codecogs.com/gif.latex?y-value is close to 4 (the original value), the three models are equivalent, but not if we have a very small value (e.g. close to 1), where the Gaussian model is far away, compared with the other two.

Nevertheless, I still have the feeling that the choice of the distribution is not the most important problem. I used to say that the choice of the link function is clearly a more important problem. Especally in actuarial ratemaking

But I wonder if actually it is really such a big deal. If we use splines, and nonparametric transformations. Consider a Poisson regression, with some identity link function, and some spline smoothing function

> library(splines)
> plot(db)
> reg1=glm(y~bs(x),family=poisson(link=
+ "identity"),data=db)
> add_predict(reg1)
Warning message:
In bs(x, degree = 3L, knots = numeric(0), Boundary.knots = c(1,  :
  some 'x' values beyond boundary knots may cause ill-conditioned bases

(since we have only 6 points, it is a bit stupid to consider such a complex model, but here, it is just a toy dataset, to illustrate)

With a log link function (and still some spline smoothing)

> plot(db)
> reg2=glm(y~bs(x),family=poisson(link=
+ "log"),data=db)
> add_predict(reg2)
Warning message:
In bs(x, degree = 3L, knots = numeric(0), Boundary.knots = c(1,  :
  some 'x' values beyond boundary knots may cause ill-conditioned bases

The impact on the prediction is clearly smaller than what we had, without any smoothing parameter

Again, if you’re not convinced, try to generate a larger dataset. For instance, generate a Poisson regression, with a log link.

> set.seed(1)
> x=runif(1000,0,6)
> L=exp(-0.07357+0.35076*x)
> y=rpois(1000,L)
> db=data.frame(x,y)
> plot(db)

And fit a Gaussian regression, with an identity link, but with some spline smoothing function,

> reg1=glm(y~x,family=poisson(link="log"),
+ data=db)
> y1=predict(reg1,newdata=nd,type="response")
> lines(nd$x,y1,col="red",lwd=2)
> reg2=glm(y~bs(x),family=gaussian(link=
+ "identity"),data=db)
> y2=predict(reg2,newdata=nd,type="response")
> lines(nd$x,y2,col="blue",lwd=2)

So it looks like neither the (distribution) law, nor the link function, have an impact on a prediction. Sad, isn’t it?


7 thoughts on “I Fought the (distribution) Law (and the Law did not win)”

  1. Well, we know that any continuous function can be well approximated by polynomials. Indeed, this is the reason behind the nonparametric sieves, right?

    I wonder if the same would happen if you have discrete regressors, or the regressor dimension increases.

    Interesting post!

  2. That’s an entertainingly subversive and thought-provoking story – thanks. Can you recommend any published studies that explore this ?

  3. Not sad — It is a cool demonstration . The distribution and the link function do not matter much when you are within the range of your data, but it seems that if you try to extrapolate beyond the range of your data, they do matter.

  4. Super article

    J’aurais néanmoins 2 petites questions.

    Le choix de la fonction de variance n’impacte pas l’estimation mais semble relativement impactante sur sa volatilité Quid des tests de significativité pour sélectionner les variables d’un modèle? La fonction variance n’a t-elle aucun impact sur les p-value des tests ?

    Sur la fonction de lien, l’exemple montre nettement son faible impact, mais on ne considère qu’un seul régresseur. Si on envisage plusieurs variables explicatives qui sont liées à la variable réponse de manière additive ou multiplicative, j’aurais tendance à penser que la fonction de lien devient relativement déterminante?

    Bonne journée et encore merci pour le super travail sur votre blog.

    1. ce sont des très bonnes questions (qui vont nécessiter d’autres billets !!)
      oui, clairement, il y a un impact sur la significativité, c’est pour ça que j’ai toujours voulu rajouter les intervalles de confiance, pour ne pas trop mentir sur les conclusions (et éviter des simplifier trop outrageusement !)

Leave a Reply

Your email address will not be published. Required fields are marked *