# Actuariat de l’Assurance non-Vie #7

Après la modélisation de la fréquence de sinistres, nous allons commencer les modèles de coûts individuels des sinistres. Les slides sont en ligne.

# Actuariat de l’Assurance non-Vie #4

Lundi prochain, suite du cours d’actuariat de l’assurance non-vie. Nous avons terminé la partie sur la classification (modèle logistique, arbres, forêts, bagging, etc), et nous allons aborder la section sur la modélisation de la fréquence, et la régression de Poisson. Je rajouter quelques slides sur la présentation des GLM, qui seront utiles pour parler un peu de sur-dispersion,

# Regression on factors

Most of our intuitions about regression models come from the Gaussian standard linear model. One interesting feature is that, when we have a factor explanatory variable, the sum of predictions per class is the sum of observations of the endogeneous variable, per class. To be more specific, consider some factor variable $x_1\in\{0,1\}$, and a regression model

$y_i=\beta_0+\beta_1 \boldsymbol{1}(x_1=1)+\beta_2 x_2+\varepsilon_i$

Use ordinary least squares to fit that model

$\widehat{y}_i=\widehat{\beta}_0+\widehat{\beta}_1 \boldsymbol{1}(x_1=1)+\widehat{\beta}_2 x_2$

Then for all $x\in\{0,1\}$

$\sum_{i:x_i=x} y_i = \sum_{i:x_i=x} \widehat{y}_i$

> n=200 > X1=rep(0:1,each=n/2) > set.seed(1) > X2=runif(2*n) > L=X1-X2 > B=data.frame(Y=rnorm(n,L),X1=as.factor(X1),X2=X2) > pd=aggregate(x=B$Y,by=list(B$X1),mean)$x > pd [1] -0.4881735 0.5341301 > fit=lm(Y~X1+X2,data=B) > B2=data.frame(x=B$X1,y=predict(fit)) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] -0.4881735 0.5341301 # Data Science for Actuaries, Regression Models with R After an introduction to Advanced R, we will discuss for the last part of our crash course visualization and graphs (from the previous set of slides), and I just uploaded additional slides on regression models (including some pdf version) # Simple Distributions for Mixtures? The idea of GLMs is that given some covariates$X$$Y|X$ has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that $Y$ has a similar distribution… so there is no reason to test for a Gamma model for $Y$ before running a Gamma regression, for instance. But are there cases where it might work? That the non-conditional distribution is the same (same family at least) than the conditional ones? For instance, if $(X,Y)$ has a joint Gaussien distribution, then both marginals are Gaussian, but also $Y|X$. So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for $Y$. The econometric interpretation is that with a standard Gaussian linear model, if$X$ is normally distributed, not only the conditional distribution $Y|X$ is Gaussian but also the non-conditional distribution of $Y$. > set.seed(1) > n=1e3 > X=rnorm(n,10,2) > Y=1+3*X+rnorm(n) > plot(X,Y,xlim=c(4,20)) Indeed, here the distribution of $Y$ is also Gaussian > library(nortest) > ad.test(Y) Anderson-Darling normality test data: Y A = 0.23155, p-value = 0.802 > shapiro.test(Y) Shapiro-Wilk normality test data: Y W = 0.99892, p-value = 0.8293 (not only from a statistical point of view, the thoery of Gaussian random vectors confirms that the non-conditional distribution is Gaussian actually) Here $X$ is continuous. What if we consider a finite mixture here, i.e.$X$ takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for $Y$. But in practice, would we really reject the Gaussian assumption, for $Y$? If the number of classes is to small, yes. But with a large number of classes (a sufficiently large number of mixture components), it is possible, > pv=function(k=2){ + n=1e4 + X=rnorm(n,10,2) + Q=quantile(X,(0:k)/k) + Q[1]=0 + Xc=cut(X,Q,labels=1:k) + XcN=tapply(X,Xc,mean) + Xn=XcN[as.numeric(Xc)] + Y=1+3*Xn+rnorm(n) + ad.test(Y)$p.value}

> plot(2:100,Vectorize(pv)(2:100),type="l")
> abline(h=.05,col="red")

So here, it could be possible to have also a Gaussian distribution, for $Y$. As least to accept that assumption, statistically.

In the context of a Poisson regression, it is well know that it’s not possible to have at the same time $Y|X$ that is Poisson distributed (that’s a Poisson regression) and also $Y$ that is Poisson distributed. That simply comes from the fact that

$\mathbb{E}[Y]=\mathbb{E}[\mathbb{E}[Y\vert X]]$

while

$\text{Var}[Y]=\mathbb{E}[\text{Var}[Y\vert X]]+\text{Var}[\mathbb{E}[Y\vert X]]$

and because of the conditional Poisson distribution, then

$\text{Var}[Y\vert X]=\mathbb{E}[Y\vert X]$

Thus,

$\text{Var}[Y]=\mathbb{E}[Y]+\underbrace{\text{Var}[\mathbb{E}[Y\vert X]]}_{>0}$

So $Y$ cannot be Poisson distribution. But again, it could be possible, if heterogeneity is not too large, to accept the null assumption of a Poisson distribution for $Y$.

More generally, it is very difficult to have a distribution family for  $Y|X$ that is also the distribution of the non-conditional variable $Y$. In the context of a finite mixture ($X$ takes a finite number of values),Teicher (1963) proved that it was not not possible, neither for the Gaussian distribution nor the Gamma distribution. An to go further, check Monfrini (2002) (thanks Romuald for point out the reference).

Hence, as a keep saying, before running a regression model on$Y|X$ with some given family, it is never a good idea to check if the non-conditional distribution $Y$ has the same distribution. Because there is no reason, usually, to remain in the same family.

# Actuariat de l’Assurance Non-Vie #7

Pour le septième chapitre du cours d’actuariat de l’assurance non-vie à l’ENSAE, on abordera la modélisation des coûts individuels, aussi bien dommage que responsabilité civile.. Les slides sont en ligne (la version pdf téléchargeable est comme souvent plus complète que celle sur slideshare)

# Actuariat de l’Assurance Non-Vie #5

Pour le quatrième cours d’actuariat de l’assurance non-vie à l’ENSAE, la semaine prochaine, on abordera la notion de sur-dispersion dans les modèles de comptage. Les slides sont en ligne (la version pdf téléchargeable est plus complète que celle sur slideshare)

# Actuariat de l’Assurance Non-Vie #4

Pour le troisième cours d’actuariat de l’assurance non-vie à l’ENSAE, on abordera la théorie des modèles linéaires généralisés Les slides sont en ligne (la version pdf téléchargeable est plus complète que celle sur slideshare)

# Actuariat de l’Assurance Non-Vie #2

Pour le second cours d’actuariat de l’assurance non-vie à l’ENSAE, qui aura lieu lundi après midi, les slides présentant les modèles classiques pour prédire des variables factorielles (classification) sont en ligne,

# Computational Time of Predictive Models

Tuesday, at the end of my 5-hour crash course on machine learning for actuaries, Pierre asked me an interesting question about computational time of different techniques. I’ve been presenting the philosophy of various algorithm, but I forgot to mention computational time. I wanted to try several classification algorithms on the dataset used to illustrate the techniques

> rm(list=ls())
"http://freakonometrics.free.fr/myocarde.csv",
> levels(myocarde$PRONO)=c("Death","Survival") But the dataset is rather small, with 71 observations and 7 explanatory variables. So I decided to replicate the observations, and to add some covariates, > levels(myocarde$PRONO)=c("Death","Survival")
> idx=rep(1:nrow(myocarde),each=100)
> TPS=matrix(NA,30,10)
> myocarde_large=myocarde[idx,]
> k=23
> M=data.frame(matrix(rnorm(k*
+ nrow(myocarde_large)),nrow(myocarde_large),k))
> names(M)=paste("X",1:k,sep="")
> myocarde_large=cbind(myocarde_large,M)
> dim(myocarde_large)
[1] 7100   31
> object.size(myocarde_large)
2049.064 kbytes

The dataset is not big… but at least, it does not take 0.0001 sec. to run a regression.  Actually, to run a logistic regression, it takes 0.1 second

> system.time(fit< glm(PRONO~.,
+ data=myocarde_large, family="binomial"))
user      system     elapsed
0.114       0.016       0.134
> object.size(fit)
9,313.600 kbytes

And I was surprised that the regression object was 9Mo, which is more than four times the size of the dataset. With a large dataset, 100 times larger,

> dim(myocarde_large_2)
[1] 710000     31

it takes 20 sec.

> system.time(fit<-glm(PRONO~.,
+ data=myocarde_large_2, family="binomial"))
utilisateur     système      écoulé
16.394       2.576      19.819
> object.size(fit)
90,9025.600 kbytes

and the object is ‘only’ ten times bigger.

# Modelling Occurence of Events, with some Exposure

This afternoon, an interesting point was raised, and I wanted to get back on it (since I did publish a post on that same topic a long time ago). How can we adapt a logistic regression when all the observations do not have the same exposure. Here the model is the following: ,

• the occurence of an event $Y_i^\star$ on the period $[0,1]$ is unobserved
• the occurence of an event $Y_i$ on $[0,E_i]$ is observed (as well as $E_i$)

If we assume that the ‘occurence of an event’ is the first occurence of a Poisson processus, we can prove that

$\mathbb{P}(Y^\star=0)=\mathbb{P}(Y^\star_{[0,1]}=0)=\mathbb{P}(Y^\star_{[0,E]}=0)\cdot \mathbb{P}(Y^\star_{[E,1]}=0)$

i.e. no event occur on $[0,1]$ if no event occur on $[0,E]$ and no event occur on $[E,1]$. Assuming independence between the two, we can prove that we have

$\mathbb{P}(Y^\star=0) = \mathbb{P}(Y=0)^E$

With words, it means that the probability of not having a claim in the first six months of the year is the square root of not have a claim over a year. Which makes sense.

# Choosing a Classifier

In order to illustrate the problem of chosing a classification model consider some simulated data,

> n = 500
> set.seed(1)
> X = rnorm(n)
> ma = 10-(X+1.5)^2*2
> mb = -10+(X-1.5)^2*2
> M = cbind(ma,mb)
> set.seed(1)
> Z = sample(1:2,size=n,replace=TRUE)
> Y = ma*(Z==1)+mb*(Z==2)+rnorm(n)*5
> df = data.frame(Z=as.factor(Z),X,Y)

A first strategy is to split the dataset in two parts, a training dataset, and a testing dataset.

> df1 = training = df[1:300,]
> df2 = testing  = df[301:500,]
• The Holdout Method: Training and Testing Datasets

The two datasets can be visualised below, with the training dataset on top, and the testing dataset below

> plot(df1$X,df1$Y,pch=19,col=c(rgb(1,0,0,.4),
+ rgb(0,0,1,.4))[df1$Z]) # 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… # Visualising a Classification in High Dimension, part 2 A few weeks ago, I published a post on Visualising a Classification in High Dimension, based on the use of a principal component analysis, to get a projection on the first two components. Following that post, I was wondering what could be done in the context of a classification on categorical covariates. A natural idea would be to consider a correspondance analysis, and to run a similar code. Consider here the dataset used in a recent post, > source("http://freakonometrics.free.fr/import_data_credit.R") If we consider a correspondance analysis, we get > library(FactoMineR) > acm=MCA(train.db,quali.sup = + which(names(train.db,)=="class"),ncp=10) For the covariates (including also the variable we want to model, considered here as some supplementary variable), the visualisation – on the first two components – is and for the individuals # Classification with Categorical Variables (the fuzzy side) The Gaussian and the (log) Poisson regressions share a very interesting property, $\frac{1}{n}\sum_{i=1}^n \widehat{Y}_i=\frac{1}{n}\sum_{i=1}^n Y_i$ i.e. the average predicted value is the empirical mean of our sample. > mean(predict(lm(dist~speed,data=cars))) [1] 42.98 > mean(cars$dist)
[1] 42.98

One can prove that it is also the prediction for the average individual in our sample

> predict(lm(dist~speed,data=cars),
+ newdata=data.frame(speed=mean(cars$speed))) 42.98 The geometric interpretation is that the regression line passes through the centroid, > plot(cars) > abline(lm(dist~speed,data=cars),col="red") > abline(h=mean(cars$dist),col="blue")
> abline(v=mean(cars$speed),col="blue") > points(mean(cars$speed),mean(cars\$dist))

But in all other cases, it is no longer the case. Consider for instance the case of a logistic regression. And to ask for something even more complicated, consider the case where we have only categorical explanatory variables. In that context, it is more difficult to get a prediction for the “average individual”. Unless we consider some fuzzy interpretation of the regression.