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)

# Tag Archives: GLM

# Simple Distributions for Mixtures?

The idea of GLMs is that given some covariates, has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that has a similar distribution… so there is no reason to test for a Gamma model for 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 has a joint Gaussien distribution, then both marginals are Gaussian, but also . So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for . The econometric interpretation is that with a standard Gaussian linear model, if is normally distributed, not only the conditional distribution is Gaussian but also the non-conditional distribution of .

> 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 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 is continuous. What if we consider a finite mixture here, i.e. takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for . But in practice, would we really reject the Gaussian assumption, for ? 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 . 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* that is Poisson distributed (that’s a Poisson regression) and also that is Poisson distributed. That simply comes from the fact that

while

and because of the conditional Poisson distribution, then

Thus,

So 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 .

More generally, it is very difficult to have a distribution family for that is also the distribution of the non-conditional variable . In the context of a finite mixture ( 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 with some given family, it is never a good idea to check if the non-conditional distribution 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()) > 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.

# 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 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

# 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…

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

# 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

Continue reading Visualising a Classification in High Dimension, part 2

# Classification with Categorical Variables (the fuzzy side)

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)

# Regression Models, It’s Not Only About Interpretation

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

# On Some Alternatives to Regression Models

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)

# Visualising a Classification in High Dimension

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