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)

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)

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,

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()) > myocarde=read.table( "http://freakonometrics.free.fr/myocarde.csv", head=TRUE,sep=";") > 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.

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 on the period is unobserved
- the occurence of an event on is observed (as well as )

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

i.e. no event occur on if no event occur on and no event occur on . Assuming independence between the two, we can prove that we have

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.

Continue reading Modelling Occurence of Events, with some Exposure

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])

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

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

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

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

Continue reading Visualising a Classification in High Dimension, part 2

The Gaussian and the (log) Poisson regressions share a very interesting property,

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.

Continue reading Classification with Categorical Variables (the fuzzy side)

Yesterday, I did upload a post where I tried to show that “standard” regression models where not performing bad. At least if you include splines (multivariate splines) to take into accound joint effects, and nonlinearities. So far, I do not discuss the possible high number of features (but with boostrap procedures, it is possible to assess something related to variable importance, that people from machine learning like).

But my post was not complete: I was simply plotting the prediction obtained by some model. And it “looked like” the regression was nice, but so were the random forrest, the -nearest neighbour and boosting algorithm. What if we compare those models on new data?

Continue reading Regression Models, It’s Not Only About Interpretation

When you start discussing with people in machine learning, you quickly hear something like *“forget your econometric models, your GLMs, I can easily find a machine learning ‘model’ that can beat yours”. *I am usually very sceptical, especially when I hear “*easily*” or “*always*“. I have no problem about the fact that I use old econometric models, but I had the feeling that things aren’t that easy. I can understand that we might have problems when we do have a lot of features (I am still working on that, I’ll get back to this point soon), but **I have the feeling that I can still capture interactions, and non-linearities with standard econometric models as well as any machine learning algorithm**.

Just to illustrate, consider the following ‘*model*‘

where is (just to illustrate)

> n <- 5000 > rtf <- function(x1, x2) { sin(x1+x2)/(x1+x2) } > xgrid <- seq(1,6,length=31) > ygrid <- seq(1,6,length=31) > zgrid <- outer(xgrid,ygrid,rtf) > persp(xgrid,ygrid,zgrid,theta=30, phi=30, + col="green", ticktype="detailed",shade=TRUE)

So far, when discussing classification, we’ve been playing on my toy-dataset (actually, I should no claim it’s mine, it is inspired by the one used in the introduction of Boosting, by Robert Schapire and Yoav Freund). But in ral life, there are more observations, and more explanatory variables.With more than two explanatory variables, it starts to be more complicated to visualise. For instance, consider

MYOCARDE=read.table( "http://freakonometrics.free.fr/saporta.csv", head=TRUE,sep=";")

where we have observations from people in E.R., for infarctus, and we want to understand who did survive, to get a predictive model. But before running some classifier, let us visualise our data. Since we have seven explanatory variables and our class (survival or death), we can go for a PCA.

library(FactoMineR) # ACP (sur les var continues) X=MYOCARDE[,1:7] acp=PCA(X)

To add the death/survival variable, treat it as numerical 0/1 variable (at least to get a direction)

MYOCARDE2=MYOCARDE MYOCARDE2$PRONO=(MYOCARDE2$PRONO=="SURVIE")*1 acp=PCA(MYOCARDE2,quanti.sup=8,graph=TRUE)

The nice thing is that we see here where variables are colinear with that one. It is also possible to visualise individuals, and classes, too

acp=PCA(MYOCARDE,quali.sup=8,graph=TRUE) plot(acp, habillage = 8,col.hab=c("red","blue"))

Continue reading Visualising a Classification in High Dimension

We will start, in our Data Science course, to discuss classification techniques (in the context of supervised models). Consider the following case, with 10 points, and two classes (red and blue)

> clr1 <- c(rgb(1,0,0,1),rgb(0,0,1,1)) > clr2 <- c(rgb(1,0,0,.2),rgb(0,0,1,.2)) > x <- c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) > y <- c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) > z <- c(1,1,1,1,1,0,0,1,0,0) > df <- data.frame(x,y,z) > plot(x,y,pch=19,cex=2,col=clr1[z+1])

To get a prediction, i.e. a partition of the space in two parts, consider some logistic regression

> reg=glm(z~x+y,data=df,family=binomial) > summary(reg) Call: glm(formula = z ~ x + y, family = binomial, data = df) Deviance Residuals: Min 1Q Median 3Q Max -1.6593 -0.4400 0.2564 0.5830 1.5374 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.706 1.999 -0.854 0.393 x -5.489 5.360 -1.024 0.306 y 8.568 5.515 1.554 0.120 (Dispersion parameter for binomial family taken to be 1) Null deviance: 13.4602 on 9 degrees of freedom Residual deviance: 8.1445 on 7 degrees of freedom AIC: 14.144 Number of Fisher Scoring iterations: 5

Given some point, the predicted class is obtained using

> pred_1 <- function(x,y){ + predict(reg,newdata=data.frame(x=x, + y=y),type="response")>.5 + }

(here, the predicted class is simply the one that is the most likely). To visualize it use

> x_grid<-seq(0,1,length=101) > y_grid<-seq(0,1,length=101) > z_grid <- outer(x_grid,y_grid,pred_1) > image(x_grid,y_grid,z_grid,col=clr2) > points(x,y,pch=19,cex=2,col=clr1[z+1])

Since the logistic regression is a (generalized) linear model, the line that separate the two regions is a straight line.

Continue reading Supervised Classification, Logistic and Multinomial

Most of the time, when we introduce binomial models, such as the logistic or probit models, we discuss only Bernoulli variables, . This year (actually also the year before), I discuss extensions to multinomial regressions, where is a function on some simplex. The multinomial logistic model was mention here. The idea is to consider, for instance with three possible classes

the following model

and

Tomorrow, around noon, I will be giving a talk on * predictive modeling* for actuaries. In the introduction, I will get back shortly on the idea that a

it is natural to use *least square* ideas. In order to illustrate all those concepts, we will use a simple dataset, with the sex, the height and the weight of a person, as well as declared weight.

Davis=read.table( "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")

Since there is a typo in this dataset, we have to invert to figures

Davis[12,c(2,3)]=Davis[12,c(3,2)]

but it’s not a big deal. The variable of interest, here, is someone’s weight

attach(Davis) Y=weight*2.204622

(here in pounds). We will use explanatory variables such as the sex of that person, or his/her height

X=Davis$height / 30.48

(in inches). So, we will start with the (standard) linear model, just to make sure that we all talk about the same thing.

The goal will be to use (possible) explanatory variable to improve our prediction. We will start with the standard linear model, but we will see that nonlinear models can also easily be obtained,

Non linearities will be discussed. But those models are Gaussian (as mentioned above). And homoscedastic. So we will see how generalized linear models can be used to model the mean and the variance, at the same time. For instance, with a Poisson regression (below), the variance will increase with the expected value.

After this general introduction, we will spend some time on 0-1 variables. We will see how to use a logistic regression, and also discuss more generally which kind of models can be used for classification. ROC curves will be presented, and explained.

Then, we will also see an alternative to the logistic model, namely classification trees and CART techniques

We will also discuss random forrests, bagging and boosting techniques

A pdf version of the slides can be downloaded.

I will give, with Qichun Xu, a joint webinar for the *Reinsurance Council* and the *Futurism Council* of the Society of Actuaries, on **Perspectives of Predictive Modeling with Case Studies **in a few days. The slides of my talk are now available (I do recommand to open the pdf version of the slides with Acrobat, since there are animated pictures in the slides that could not be visualized below for instance). The Society of Actuaries asked specifically for a *powerpoint* document, so I will use screenshots of the slides for the webinar. I do encourage to open and read the pdf file for a better quality… Sorry for the inconvenience. I will upload soon lines of codes to reproduce most of the graphs. All comments and remarks are welcome.