To content | To menu | To search

Tuesday, May 15 2012

Finding Waldo, a flag on the moon and multiple choice tests, with R

I have to admit, first, that finding Waldo has been a difficult task. And I did not succeed. Neither could I correctly spot his shirt (because actually, it was what I was looking for). You know, that red-and-white striped shirt. I guess it should have been possible to look for Waldo's face (assuming that his face does not change) but I still have problems with size factor (and resolution issues too). The problem is not that simple. At the conference, a price was offered for writing an algorithm in Matlab. And one can even find Mathematica codes online. But most of the those algorithms are based on the idea that we look for similarities with Waldo's face, as described in problem 3 on's webpage. You can find papers on that problem, e.g. Friencly & Kwan (2009) (based on statistical techniques, but Waldo is here a pretext to discuss other issues actually), or more recently (but more complex) Garg et al. (2011) on matching people in images of crowds.

What about codes in R ? On, some ideas can be found (and thank Robert Hijmans for his help on his package). So let us try here to do something, on our own. Consider the following picture,

With the following code (based on the following file) it is possible to import the picture, and to extract the colors (based on an RGB decomposition),
> library(raster)
> waldo=brick(system.file("DepartmentStoreW.grd",
+ package="raster"))
> waldo
class       : RasterBrick
dimensions  : 768, 1024, 786432, 3 (nrow,ncol,ncell,nlayer)
resolution  : 1, 1  (x, y)
extent      : 0, 1024, 0, 768  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values      : C:\R\win-library\raster\DepartmentStoreW.grd
min values  : 0 0 0
max values  : 255 255 255
My strategy is simple: try to spot areas with white and red stripes (horizontal stripes). Note that here, I ran the code on a Windows machine, the package is not working well on Mac. In order to get a better understanding of what could be done, let us start with something much more simple. Like the picture below, with Waldo (and Waldo only). Here, it is possible to extract the three colors (red, green and blue),
> plot(waldo,useRaster=FALSE)


It is possible to extract the red zones (already on the graph above, since red is a primary color), as well as the white ones (green zones on the graphs means a white region on the picture, on the left)
# white component
white = min(waldo[[1]] , waldo[[2]] , waldo[[3]])>220
focalswhite = focal(white, w=3, fun=mean)
# red component
red = (waldo[[1]]>150)&(max(  waldo[[2]] , waldo[[3]])<90)
focalsred = focal(red, w=3, fun=mean)
i.e. here we have the graphs below, with the white regions, and the red ones,

From those two parts, it has been possible to extract the red-and-white stripes from the picture, i.e. some regions that were red above, and white below (or the reverse),
# striped component
striped = red; n=length(values(striped)); h=5
focalsstriped = focal(striped, w=3, fun=mean)
So here, we can easily spot Waldo, i.e. the guy with the red-white stripes (with two different sets of thresholds for the RGB decomposition)

 Let us try somthing slightly more complicated, with a zoom on the large picture of the department store (since, to be honest, I know where Waldo is...).

Here again, we can spot the white part (on the left) and the red one (on the right), with some thresholds for the RGB decomposition

Note that we can try to be (much) more selective, playing with threshold. Here, it is not very convincing: I cannot clearly identify the region where Waldo might be (the two graphs below were obtained playing with thresholds)

 And if we look at the overall pictures, it is worst. Here are the white zones, and the red ones,

and again, playing with RGB thresholds, I cannot spot Waldo,

Maybe I was a bit optimistic, or ambitious. Let us try something more simple, like finding a flag on the moon. Consider the picture below on the left, and let us see if we can spot an American flag,

Again, on the left, let us identify white areas, and on the right, red ones

Then as before, let us look for horizontal stripes

Waouh, I did it ! That's one small step for man, one giant leap for R-coders ! Or least for me... So, why might it be interesting to identify areas on pictures ? I mean, I am not Chloe O'Brian, I don't have to spot flags in a crowd, neither Waldo, nor some terrorists (that might wear striped shirts). This might be fun if you want to give grades for your exams automatically. Consider the two following scans, the template, and a filled copy,

A first step is to identify regions where we expect to find some "red" part (I assume here that students have to use a red pencil). Let us start to check on the template and the filled form if we can identify red areas,
exam = stack("C:\\Users\\exam-blank.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=3, fun=mean)
exam = stack("C:\\Users\\exam-filled.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=3, fun=mean)

First, we have to identify areas where students have to fill the blanks. So in the template, identify black boxes, and get the coordinates (here manually)
exam = stack("C:\\Users\\exam-blank.png")
black = max(  exam[[1]] ,exam[[2]] , exam[[3]])<50
focalsblack = focal(black, w=3, fun=mean)

The blue points above are where we look for students' answers. Then, we have to define the vector of correct answers,
points(CORRECTX, CORRECTY,pch=16,col="red",cex=1.3)
Now, let us get back on red areas in the form filled by the student, identified earlier,
exam = stack("C:\\Users\\exam-filled.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=5, fun=mean)

Here, we simply have to compare what the student answered with areas where we expect to find some red in,
points(CORRECTX, CORRECTY,pch=1,
Crosses on the graph on the right below are the answers identified as correct (here 13),
> icorrect=values(red)[(750-CORRECTY)*
+ 610+(CORRECTX)]
> points(CORRECTX[icorrect], CORRECTY[icorrect],
+ pch=4,col="black",cex=1.5,lwd=1.5)
> sum(icorrect)
[1] 13

In the case there are negative points for non-correct answer, we can count how many incorrect answers we had. Here 4.
> iuncorrect=values(red)[(750-UNCORRECTY)*610+
> sum(iuncorrect)
[1] 4
So I have not been able to find Waldo, but I least, that will probably save me hours next time I have to mark exams...

Thursday, March 22 2012

Do we appreciate sunbathing in Spring ?

We are currently experiencing an extremely hot month in Montréal (and more generally in North America). Looking at people having a beer, and starting the first barbecue of the year, I was wondering: if we asked people if global warming was a good or a bad thing, what do you think they will answer ? Wearing a T-shirt in Montréal in March is nice, trust me ! So how can we study, from a quantitative point of view, depending on the time of the year, what people think of global warming ?

A few month ago, I went quickly through

score.sentiment = function(sentences, pos.words,
neg.words, .progress='none')
scores = laply(sentences, 
function(sentence, pos.words, neg.words) { sentence = gsub('[[:punct:]]', '', sentence) sentence = gsub('[[:cntrl:]]', '', sentence) sentence = gsub('\\d+', '', sentence) sentence = tolower(sentence) word.list = strsplit(sentence, '\\s+') words = unlist(word.list) pos.matches = match(words, pos.words) neg.matches = match(words, neg.words) pos.matches = ! neg.matches = ! score = sum(pos.matches) - sum(neg.matches) return(score) }, pos.words, neg.words, .progress=.progress ) scores.df = data.frame(score=scores, text=sentences) return(scores.df) }

hu.liu.pos = scan("positive-words.txt", what="character",
hu.liu.neg = scan('negative-words.txt', what='character',

> score.sentiment("It's awesome I am so happy,
thank you all",
+ hu.liu.pos,hu.liu.neg)$score
[1] 3

> score.sentiment("I'm desperate, life is a nightmare,
I want to die",
+ hu.liu.pos,hu.liu.neg)$score
[1] -3

But one can easy see a big problem with this methodology. What if the sentence included negations ? E.g.

> score.sentiment("I'm no longer desperate, life is
not a nightmare anymore I don't want to die",
+ hu.liu.pos,hu.liu.neg)$score
[1] -3

Here the sentence is negative, extremely negative, if we look only at the score. But it should be the opposite. I simple idea is to change (slightly) the function, so that once a negation is found in the sentence, we take the opposite of the score. Hence, we just add at the end of the function


Here we obtain

> score.sentiment.neg("I'm no longer desperate,
life is not a nightmare anymore I don't want to die",
+ hu.liu.pos,hu.liu.neg)$score
[1] 3

But does it really work ? Let us focus on Tweets,


Consider the following tweet-extractions, based on two words, a negative word, and the negation of a positive word,

> tweets=searchTwitter('not happy',n=1000)
> NH.text= lapply(tweets, function(t) t$getText() )
> NH.scores = score.sentiment(NH.text,
+ hu.liu.pos,hu.liu.neg)
> tweets=searchTwitter('unhappy',n=1000)
> UH.text= lapply(tweets, function(t) t$getText() )
> UH.scores = score.sentiment(UH.text,
+ hu.liu.pos,hu.liu.neg)

> plot(density(NH.scores$score,bw=.8),col="red")
> lines(density(UH.scores$score,bw=.8),col="blue")

> UH.scores = score.sentiment.neg(UH.text,
+ hu.liu.pos,hu.liu.neg)
> NH.scores = score.sentiment.neg(NH.text,
+ hu.liu.pos,hu.liu.neg)
> plot(density(NH.scores$score,bw=.8),col="red") > lines(density(UH.scores$score,bw=.8),col="blue")

> w.tweets=searchTwitter("snow",since= LISTEDATE[k],
+ until= LISTEDATE[k+1],geocode="40,-100,2000mi")
> W.text= lapply(w.tweets, function(t) t$getText() )
> W.scores = score.sentiment.neg(W.text,
+ hu.liu.pos,hu.liu.neg, .progress='text')
> M[k]=mean(W.scores$score)
We obtain here the following score function, over three years, on Twitter,

Well, we have to admit that the pattern is not that obvious. There might me small (local) bump in the winter, but it is not obvious...

Let us get back to the point used to introduce this post. If we study what people "feel" when they mention global warming, let us run the same code, again in North America
> w.tweets=searchTwitter("global warming",since= LISTEDATE[k],
+ until= LISTEDATE[k+1],geocode="40,-100,2000mi")
Actually, I was expecting a nice cycle, with positive scores in Spring, and perhaps negative scores during heat waves, in the middle of the Summer...

What we simply observe is that global warming was related to "negative words" on Twitter a few years ago, but we have reached a neutral position nowadays.

And to be honest, I do not really know how to interpret that: is there a problem with the technique I use (obviously, I use here a very simple scoring function, even after integrating a minor correction to take into consideration negations) or is there something going one that can be interpreted ?

Monday, March 19 2012

Simulation d'un processus de Lévy, et discrétisation

Avec @renaudjf, on discutait l'autre jour de la simulation d'un processus de Lévy. Et on se posait la question d'un algorithme optimal pour combiner un processus de Poisson (ou un process Poisson composé) avec un processus de Wiener (avec éventuellement un drift, voire une diffusion plus générale). En fait, pour générer des processus de Poisson, j'ai toujours eu l'habitude de simuler les durées entre sauts (avec des lois exponentielles, indépendantes, comme dans des vieux billets). Jean François me suggérait d'utiliser une propriété d'uniformité des sauts sur un intervalle de temps donné, conditionnellement aux nombres de sauts.

Commençons par la première piste. On peut générer un processus de Wiener, éventuellement avec un drift, et à coté, on peut générer les lois exponentielles  (qui vont correspondre aux durées entre sauts), et éventuellement aussi des amplitudes de sauts (e.g. des pertes qui suivent des lois exponentielles). On a ici

où . On commence par générer  en notant que

où les incréments  sont Gaussiens (centrés et de variance ) et indépendants les uns des autres. Quant aux durées entre sauts, les  , ce sont des lois exponentielles indépendantes, de moyenne . Voilà le code qui permet de générer les trois composantes,


Le hic est que pour le processus de Wiener, on a du discrétiser, alors que pour le processus de Poisson composé, non. Pourtant, il va bien se ramener sur une même échelle de temps. Une première piste est de créer vraiment la fonction 


et pour faire un dessin ensuite, c'est un jeu d'enfant


Une autre possibilité est d'utiliser une propriété d'uniformité du processus de Poisson que j'évoquais en introduction. Car le processus de Poisson vérifie une propriété remarquable: si  est la date où survient le ième saut, , alors conditionnellement au fait que , les variables  correspondent aux statistiques d'ordres de  variables indépendantes, uniformément distribuées sur , i.e.

Cette propriété se trouve dans Wolff (1982). L'idée de la démonstration est relativement simple. Commençons avec un (unique) saut, alors pour ,

i.e. on retrouve la fonction de répartition d'une loi uniforme sur . On itère ensuite avec 2 sauts, 3 sauts, etc.

La traduction en R de cette idée est tout simplement (car on se place sur )


Ensuite, une stratégie est de discrétiser le processus de Poisson, avec le même pas de temps que le processus de Wiener,


On retrouve la même trajectoire qu'auparavant



Sauf qu'on a eu de la chance. Avec cette procédure, il ne faut pas que l'on ait deux sauts dans le même intervalle de temps ! Bon, il est vrai qu'une caractérisation du processus de Poisson est que

et donc on doit avoir très peu de chance d'avoir deux sauts au même instant d'autant plus que le pas de temps est petit. Mais "peu de chance" ne veut pas dire nul, et si on génère des milliers de trajectoires, la probabilité d'avoir une fois un soucis n'est pas négligeable.

Jean-François a eu l'idée brillante de proposer de tirer non pas des lois uniformes sur , mais des lois uniformes discrètes, dans



sans remise afin d'éviter d'avoir deux sauts au même moment. L'idée était séduisante, mais il me restait à etre convaincu que les durées entre sauts... continuaient à être (presque) exponentielles.

Pour ça, on peut faire quelques tests. Par exemple, générer quelques simulations pour avoir une centaines de sauts (et donc une centaine de durées entre sauts), puis faire un test de loi exponentielle (de moyenne )

for(ns in 1:20){

On fait ici 20 boucles car on avait fixé


et j'avais dit que je voulais une centaine d'observations pour faire un test de loi (ce qui est purement arbitraire, on en conviendra). On peut ensuite faire un test d'ajustement de loi exponentielle, 


Si on répète un grand nombre de fois, en changeant le pas de temps (ou le nombre de subdivisons de l'intervalle de temps), on notera qu'effectivement, avec un grand pas de temps (à gauche ci-dessous) on rejette souvent, voire presque toujours, l'hypothèse de loi exponentielle. Mais qu'assez vite, c'est une hypothèse qui n'est pas invraisemblable,

Je ne sais pas entre ces deux fonctions ce qui serait le plus rapide, mais on a deux jolis algorithmes pour générer les processus de Lévy. Non ?

Thursday, March 1 2012

Elections, sondages et intervalles de confiance

(ou, pour parapher le titre d'un article paru l'an passé dans Le Monde "Sondages : et si l'on prenait sérieusement en compte les marges d'erreur ?"). Il y a plusieurs mois, suite aux questions d' et @adelaigue , je m'étais plongé dans la problématique des intervalles de confiance quand trois possibilités sont offertes à des personnes sondées (dans un billet en français, et quelques jours plus tard un billet plus complet en anglais, avec quelques typos qu'il faudrait reprendre...).

L'idée était que quand on a le choix entre deux candidats, l'intervalle de confiance associé aux probabilités de voter pour chacun des candidats est de la forme désigne la fonction de répartition de la loi normale centrée réduite (valide compte tenue de l'approximation de la loi binomiale par une loi normale). C'est d'ailleurs cette formule que l'on retrouve dans Le Monde

> n=1000
> p=20/100
> alpha=5/100
> qnorm(1-alpha/2)*sqrt(p*(1-p)/n)
[1] 0.0247918

Si on regarde ce qui est communément fait, c'est de prendre la borne supérieure de cet intervalle de confiance,

(ce qui peut être très conservateur). On peut alors faire toutes sortes de tests comme, contre C'est assez classique. Sauf que si on a trois candidats, on peut se demander si le précédant test change, ou pas... Et cet intervalle de confiance permet de dire quoi que ce soir. Pour Le Monde la réponse semble claire, puisque tous les exemples évoquent trois candidats, puis l'article évoque ensuite (toujours) des intervalles de confiance construits dans un cas binomial... Le problème est un peu le même que celui des ellipsoïdes de confiance dans le cas Gaussien (les régions de confiances sont alors assez éloignés de pavés obtenus en s'intéressant uniquement à des intervalles de confiance obtenus de manière univariée, comme évoqué dans un précédant billet).

Donc formellement, on a maintenant non plus et (avec mais, et (avec Se placer ainsi sur la bordure du simplexe en dimension 3 permet de faire des jolis dessins. En effet, la bordure du simplexe forme un triangle équilatéral,

Or les triangles équilatéraux sont des figures très particulières: si on considère un point intérieur, la somme des distances à chacun des cotés est constante. C'est le théorème de Viviani.

(ce qui se démontre soit en considérant un triangle translaté, soit en travaillant sur les aires, comme le montre l'animation ci-dessus). On peut alors représenter un triplet dans le triangle,

Les probabilités étant alors les distances à chacun des côtés. Plus on est proche d'un sommet, plus les intentions de votes pour ce candidat sont importantes. Chercher lequel parmi deux candidats (sur trois) aura le plus de voix revient par exemple à s'interroger dans quelle partie du triangle on se trouve,

Bref, avoir le choix parmi trois, et non plus deux candidats est un peu plus compliqué, mais on peut toujours faire des dessins... Et des calculs. Si on utilise des intervalles de confiance obtenus dans le cas binomial, on aurait quelque chose du genre (ce qui revient à utiliser en régression des hypercubes pour les régions de confiance des estimateurs des paramètres, alors qu'on peut avoir des ellipsoïdes de taille beaucoup plus petite, comme détaillé dans un précédant billet)

En fait, on peut faire un peu mieux probablement... Si on néglige la corrélation qui existe entre les estimateurs, et que l'on continue à utiliser une approximation Gaussienne,on aurait

Numériquement, avec l'échantillon suivant

> alpha
[1]  20 20  50

on obtient graphiquement la distribution suivante pour le triplet

Si maintenant on corrige, en tenant compte du fait que les estimateurs des fréquences son négativement corrélés (comme expliqué dans un vieux billet)

Toujours sur le même exemple numérique, on obtient graphiquement la distribution suivante pour le triplet,

Plus généralement, on peut supposer qu'il y a candidats, et que les probabilités de voter pour chacun d'entre eux est, où désigne la frontière (supérieure) du simplexe de, i.e.


Un sondage auprès de n personnes peut être vu comme un tirage d'une loi multinomiale Si est le nombre de personnes qui ont déclarer voter pour chacun des candidats (on demande aux sondés de choisir un seul candidats parmi les, i.e. et alors

que l'on peut réécrire

avec Histoire de retrouver une loi de la famille exponentielle. Au lieu d'avoir une approche fréquentiste, on peut - pour changer - proposer un estimateur bayésien. Car on dispose de lois simples sur les simplexes: les loi de Dirichlet. Lois qui pourront servir de loi a priori pour les probabilités. En plus, c'est la loi conjuguée de la loi multinomiale... On supposera ici comme loi a priori

pour et
L'estimateur de Bayes obtenu en considérant une fonction de perte quadratique sera alors

On retrouve une grandeur qui est linéaire en la moyenne empirique, car on est dans un cas de crédibilité linéaire, à la Bühlmann pour reprendre une terminologie chère aux actuaires (on est dans des familles exponentielles, avec une version multivariée du modèle binomial-beta). Et plus généralement, la loi a posteriori pour est une loi de Dirichlet ! Sur l'exemple précédant, si on utilise une loi de Dirichlet

>  D=ddirichlet(p, alpha)

on obtient graphiquement la distribution suivante pour le triplet ,

Comme on a la loi de notre estimateur, on peut faire des tests, ou estimer des probabilités. Par exemple, si on a observé les résultats suivants à un sondage,

> alpha
[1]  40 100  50

et qu'on se demande quelle peut être la probabilité que, on peut utiliser trois modèles (ou hypothèses),

  • une approximation gaussienne, et supposer les estimateurs comme étant indépendants
  • une approximation gaussienne, et prendre en compte la corrélation entre les estimateurs
  • utiliser un modèle bayésien, et utiliser une loi de Dirichlet
>  SIGMA=matrix(0,3,3)
>  diag(SIGMA)=alpha*(sum(alpha)-alpha)/sum(alpha)^2
>  RI=rmnorm(ns, mean=(alpha/sum(alpha))[1:2],
+ varcov=(SIGMA/sum(alpha))[1:2,1:2]) > SIGMA=-alpha%*%t(alpha)/sum(alpha)^2 > diag(SIGMA)=alpha*(sum(alpha)-alpha)/sum(alpha)^2 > RG=rmnorm(ns, mean=(alpha/sum(alpha))[1:2],
+ varcov=(SIGMA/sum(alpha))[1:2,1:2]) > RD=rdirichlet(ns, alpha) > mean(1-(RI[,1]+RI[,2])>=RI[,1]) [1] 0.7759848 > mean(1-(RG[,1]+RG[,2])>=RG[,1]) [1] 0.8548025 > mean(RD[,3]>=RD[,1]) [1] 0.8554395

On note que les deux dernières approches sont relativement proches, mais que négliger la corrélation entre les estimateurs conduit à une forte sous-estimation de la probabilité.

On peut bien entendu aller au delà de la dimension 3, mais on perdra alors la représentation graphique. Par exemple dans l'article de Le Monde, on est en réalité en dimension 4 car en plus des 3 candidats qui semblent intéressant, soit une autre réponse possible est proposée, soit on autoriser la non-réponse. Si on regarde par exemple le premier sondage évoqué,

on peut se demander quelle peut être la probabilité que le dernier candidat dans le sondage soit qualifié pour le second tour, i.e. Numériquement, cela se calcule très très bien (et nul besoin de rappeler ces histoires sur les intervalles de confiances comme le fait constamment l'article, car ils ne servent à rien...)

>  alpha=c(.24,.23,.21)
>  alpha=c(alpha,1-sum(alpha))*1347
>  SIGMA=matrix(0,4,4)
>  diag(SIGMA)=alpha*(sum(alpha)-alpha)/sum(alpha)^2
>  RI=rmnorm(ns, mean=(alpha/sum(alpha))[1:3],
+ varcov=(SIGMA/sum(alpha))[1:3,1:3])
>  SIGMA=-alpha%*%t(alpha)/sum(alpha)^2
>  diag(SIGMA)=alpha*(sum(alpha)-alpha)/sum(alpha)^2
>  RG=rmnorm(ns, mean=(alpha/sum(alpha))[1:3],
+ varcov=(SIGMA/sum(alpha))[1:3,1:3])
>  RD=rdirichlet(ns, alpha)
> mean(RI[,3]>=apply(cbind(RI[,1],RI[,2]),1,min))
[1] 0.1231242
> mean(RG[,3]>=apply(cbind(RG[,1],RG[,2]),1,min))
[1] 0.162147
> mean(RD[,3]>=apply(cbind(RD[,1],RD[,2]),1,min))
[1] 0.1604497

Là encore, on observe clairement qu'oublier ces histoires de corrélations entre les estimateurs n'est pas réaliste, la probabilité que le troisième candidat soit au second tour est 35% plus élevée que celle calculée (trop) rapidement. Bref, les intervalles de confiance avec plus de deux candidats, c'est plus compliqué... mais on peut malgré tout dire plein de choses !

Friday, December 2 2011

ACT2040: lire une sortie de régression

Histoire que les choses soient bien claires, prenons 2 minutes pour revenir sur la lecture d'une sortie de régression,

>  reg=glm(nbre~carburant+zone+ageconducteur+
offset(log(exposition)), + data=baseFREQ,family=poisson(link="log")) > summary(reg)   Call: glm(formula = nbre ~ carburant + zone + ageconducteur + offset(log(exposition)), family = poisson(link = "log"), data = baseFREQ)   Deviance Residuals: Min 1Q Median 3Q Max -0.5372 -0.3643 -0.2731 -0.1503 4.5994   Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.360655 0.105711 -22.331 < 2e-16 *** carburantE -0.283956 0.043117 -6.586 4.53e-11 *** zoneB 0.122891 0.108713 1.130 0.258 zoneC 0.224469 0.088310 2.542 0.011 * zoneD 0.411243 0.086804 4.738 2.16e-06 *** zoneE 0.532070 0.088226 6.031 1.63e-09 *** ageconducteur -0.007031 0.001550 -4.537 5.70e-06 *** --- Signif. codes: 0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1   (Dispersion parameter for poisson family taken to be 1)   Null deviance: 13659 on 49999 degrees of freedom Residual deviance: 13538 on 49993 degrees of freedom AIC: 17788   Number of Fisher Scoring iterations: 6

L'écriture formelle de ce modèle serait quelque chose de la forme correspond au nombre de sinistres observés pour l'assuré , correspond à l'exposition (i.e. au temps pendant lequel l'assuré a été observé dans la base), est la première variable, ici le carburant (qui est une variable qualitative prenant 2 modalités), est la seconde variable, ici la zone géographique (qui est encore qualitative, avec 5 zones), et enfin est la troisième et dernière variable explicative, ici l'âge du conducteur principal (qui est pris ici comme une variable continue). Le modèle s'écrit, si l'on disjoncte les variables qualitatives s'écrit

Il y alors 6 paramètres à estimer. C'est ce qui est renvoyé dans la sortie de la régression,

               Estimate Std. Error z value Pr(>|z|)
(Intercept)   -2.360655   0.105711 -22.331  < 2e-16 ***
carburantE    -0.283956   0.043117  -6.586 4.53e-11 ***
zoneB          0.122891   0.108713   1.130    0.258
zoneC          0.224469   0.088310   2.542    0.011 *
zoneD          0.411243   0.086804   4.738 2.16e-06 ***
zoneE          0.532070   0.088226   6.031 1.63e-09 ***
ageconducteur -0.007031   0.001550  -4.537 5.70e-06 ***
Dans ce tableau, on a pour chaque paramètre, une estimation, par exemple vaut -2.36, alors que vaut 0.2244. On a ensuite une estimation de l'écart-type de ces estimateurs, par exemple vaut 0.086.  On peut alors utiliser un test de significative basé sur une hypothèse de normalité de ces estimateurs. On note que l'âge du conducteur est significatif dans ce modèle. Pour les facteurs, la significativité est jugée par rapport à la modalité de référence, qui est celle qui n'apparaît pas dans la sortie de la régression. Aussi, le carburant "Essence" est significatif par rapport à la modalité de référence qui est le Diesel. En revanche, la zone "B" n'est pas significative par rapport à la zone de référence qui est la zone "A", ce qui suggèrerait de regrouper les zone "A" et "B" (comme discuté dans un ancien billet). Notons que l'on peut d'ailleurs aller un peu plus loin, en utilisant des estimateurs plus robustes de l'écart-type de ces estimateurs (proposés par Cameron & Trivedi (1998)),
> library(sandwich)
> cov.reg<-vcovHC (reg, type="HC0")
> std.err<-sqrt(diag(cov.reg))
> estimation<-cbind(estimate=reg$coefficients, std.err,
+        pvalue=round(2*(1-pnorm(abs(reg$coefficients/std.err))),5),
+        lower=reg$coefficients-1.96*std.err,
+        upper=reg$coefficients+1.96*std.err)
> estimation
estimate     std.err  pvalue       lower        upper
(Intercept)   -2.360654805 0.110541249 0.00000 -2.57731565 -2.143993957
carburantE    -0.283955585 0.044492586 0.00000 -0.37116105 -0.196750117
zoneB          0.122890716 0.111122177 0.26877 -0.09490875  0.340690183
zoneC          0.224468748 0.091032422 0.01367  0.04604520  0.402892295
zoneD          0.411242820 0.089449063 0.00000  0.23592266  0.586562984
zoneE          0.532070373 0.091675282 0.00000  0.35238682  0.711753926
ageconducteur -0.007030691 0.001664597 0.00002 -0.01029330 -0.003768081
On retrouve ici des grandeurs du même ordre que tout à l'heure, avec un intervalle de confiance à 95% pour les estimateurs. Si 0 appartient à l'intervalle de confiance, on dira que le paramètre n'est pas significatif. A partir de ces estimations, il est facile de faire une prédiction, par exemple pour un assuré de 40 ans, résidant dans la zone D et conduisant un véhicule diesel, sa espérance annuelle d'accident est de 10.7%
>  predict(reg,newdata=data.frame(carburant="D",
+  zone="D",ageconducteur=40,exposition=1),type="response")
>  cbind(c(1,0,0,0,1,0,40),reg$coefficients)
[,1]         [,2]
(Intercept)      1 -2.360654805
carburantE       0 -0.283955585
zoneB            0  0.122890716
zoneC            0  0.224468748
zoneD            1  0.411242820
zoneE            0  0.532070373
ageconducteur   40 -0.007030691
>  exp(t(c(1,0,0,0,1,0,40))%*%reg$coefficients)
[1,] 0.1074597
Notons que dans ce cas, la probabilité de ne pas avoir d'accident est
>  lambda=exp(t(c(1,0,0,0,1,0,40))%*%reg$coefficients)
>  dpois(0,lambda)
[1] 0.8981127
>  exp(-lambda)
[1,] 0.8981127
Parmi les autres interprétations possibles de ces valeurs, notons qu'à zone et à âge identique, un assuré conduisant un véhicule essence à une fréquence de sinistres 25% plus faible qu'un assuré conduisant un véhicule diesel,
>  exp(reg$coefficients[2])
En bas, on a la déviance du modèle. Dans le cas d'une régression de Poisson, on sait que
soit ici
>  baseFREQ$pred=predict(reg,type="response")
>  Z=log(baseFREQ$pred^baseFREQ$nbre)
>  2*(sum(log(baseFREQ$nbre^baseFREQ$nbre)-baseFREQ$nbre)-
+  sum(Z- baseFREQ$pred))
[1] 13538.09
> deviance(reg)
[1] 13538.09
Mais la déviance n'est pas forcément pertinente pour comparer des modèles ayant des nombres différentes de paramètres, et donc la sortie inclue le critère d'Akaike.

Monday, October 3 2011

ACT2040, introduction aux modèles linéaires généralisés

On commencera ce mardi les GLM, après avoir introduit les lois exponentielles (qui ont du être revues en démonstration vendredi dernier). La notation utilisée sera que la loi (densité ou fonction de probabilité) de sera de la forme
Pour un complément plus exhaustif, je renvoie au chapitre en ligne.
  • Le modèle linéaire (Gaussien)
Le modèle de base est le modèle Gaussien que l'on avait revu au dernier cours,
> X=c(1,2,3,4)
> Y=c(1,2,5,6)
> base=data.frame(X,Y)
> reg1=lm(Y~1+X,data=base)
> nbase=data.frame(X=seq(0,5,by=.1))
> Y1=predict(reg1,newdata=nbase)
Pour une prédiction (unique), on obtient la prédiction suivante

Le code pour une telle représentation est le suivant
> plot(X,Y,pch=3,cex=1.5,lwd=2,xlab="",ylab="")
> lines(nbase$X,Y1,col="red",lwd=2)
> u=2
> mu=predict(reg1)[2]
> sigma=summary(reg1)$sigma
> y=seq(0,7,.05)
> loi=dnorm(y,mu,sigma)
> segments(u,y,loi+u,y,col="light green")
> lines(loi+u,y)
> abline(v=u,lty=2)
> points(X[2],Y[2],pch=3,cex=1.5,lwd=2)
> points(X[2],predict(reg1)[2],pch=19,col="red")
> arrows(u-.2,qnorm(.05,mu,sigma),
+ u-.2,qnorm(.95,mu,sigma),length=0.1,code=3,col="blue")
On peut multiplier les prédictions, en se basant sur l'hypothèse d'homoscédasticité (la variance sera alors constante)

Mais on peut aller plus loin
  • Le modèle linéaire généralisé
Plusieurs modèles peuvent etre estimés, en changeant la loi de la variable à expliquer, et la fonction lien,
> reg2=glm(Y~1+X,data=base,family=poisson(link="identity"))
> Y2=predict(reg2,newdata=nbase,type="response")
> reg3=glm(Y~1+X,data=base,family=poisson(link="log"))
> Y3=predict(reg3,newdata=nbase,type="response")
> reg4=glm(Y~1+X,data=base,family=gaussian(link="log"))
> Y4=predict(reg4,newdata=nbase,type="response")
> sigma=sqrt(summary(reg4)$dispersion)
Pour le modèle Poissonnien avec un lien identité, on obtient

On obtient ainsi une variance qui augmente avec la prédiction,

Pour une régression de Poisson avec un lien logarithmique,

i.e. pour nos quatre prédictions

On peut comparer avec une prédiction d'un modèle Gaussien avec un lien logarithmique,

i.e. pour les quatre prédictions

Friday, November 5 2010

Pretty R code in the blog

David Smith (alias @revodavid, see also on the Revolutions blog, here) pointed out that my R code was not easy to read (not only due to my computing skills, but mainly because of the typography I use). He suggested that I use the Pretty R tool (here). And I will...

So, just to answer quickly to a question I received by email (a few weeks ago, sorry for the delay), here is the code to get the following nice plot

library(evd); data(lossalae) 
x <- lossalae$Loss; y <- lossalae$ALAE
xhist <- hist(log(x), plot=FALSE)
yhist <- hist(log(y), plot=FALSE)
top <- max(c(xhist$counts, yhist$counts))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE),
c(3,1), c(1,3), TRUE)
plot(x, y, xlab="", ylab="",log="xy",col="red")
barplot(xhist$counts, axes=FALSE, ylim=c(0, top),
space=0,col="light green")
barplot(yhist$counts, axes=FALSE, xlim=c(0, top),
space=0, horiz=TRUE,col="light blue")
I guess it is more convenient to read it, isn't it...

Tuesday, November 2 2010

Statistique de l'assurance STT6705V, partie 9

Attention, étant donné que le changement d'heure a eu lieu en France samedi dernier, mais qu'il n'aura lieu au Canada que ce samedi, la transmission sera légèrement perturbée.. Le cours durera deux heures et pas trois: de 11:30-13:30 à Montréal (on commence une heure plus tard) et 16:30-18:30 à Rennes (car la salle n'était pas libre plus tôt).
Nous finirons la partie sur le provisionnement (je parlerais des projets) et sinon je renvoie ici pour des notes de cours détaillées sur ce que nous avons vu.
COMPLÉMENT: un paragraphe de dernière minute... je mets des lignes de code, ici. Il s'agit de ligne que nous utiliserons, et que nous commenterons pendant ce cours..

Wednesday, October 27 2010

A million ? what are the odds...

50 days ago, I published a post, here, on forecasting techniques. I was wondering what could be the probability to have, by the end of this year, one million pages viewed (from Google Analytics) on this blog. Well, initially, it was on my blog at the Université de Rennes 1 (, but since I transfered the blog, I had to revise my code. Initially, I had that kind of graphs,

and when I look at the cumulative distribution of the number of pages viewed on January first, I had

while for the distribution of the time I should read this million (the dual problem), I obtained

and I said that I should have around 35% chance to reach the million pages viewed by the end of this year.
Here is the updated graph, with the blog à Université de Rennes 1 (still in black) and the one here (in blue, where I add the two blogs together).

Actually, I decided to look at the evolution of the probability to reach the million by New Year's Eve...

The code looks like that,
> base=read.table("",
+ sep=";",header=TRUE)
> X2=cumsum(base$nombre)
> X=X1+X2
> kt=which(D==as.Date("01/06/2010","%d/%m/%Y"))
> D0=as.Date("08/11/2008","%d/%m/%Y")
> D=D0+1:length(X1)
> P=rep(NA,(length(X)-kt)+1)
> for(h in 0:(length(X)-kt)){
+ model  <- arima(X[1:(kt+h)],c(7 ,   # partie AR
+                     1,    # partie I
+                     7),method="CSS")   # partie MA
+ forecast <- predict(model,200)
+ u=max(D[1:kt+h])+1:300
+ k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
+ (P[h+1]=1-pnorm(1000000,forecast$pred[k],forecast$se[k]))
+ }
It has been a bit tricky, since I wanted an automatic fit of the ARIMA process, meaning that I had to assess a priori the orders of the ARIMA process. And I had numerical problems, since we got non stationary AR part at least at one period of time considered.... So finally I used here the CSS method which uses conditional-sum-of-squares to find starting values in the optimization procedure.

Actually, if we consider a classical descritption of traders, it looks like I act as a trader (dealing with millions and forgetting about real people): it is the same here, I do not know what a million means, I cannot imagine 250,000 visitors looking at that blog... But I can still do the maths. Anyway, a million is huge when I start to think about it... but perhaps I should not... I cannot possibility imagine that so many people might find interesting my mathematical lucubration*....
* initially I was looking for the analogous of "élucubration" in French, meaning "divagation, absurd theory" (the proper translation might be "rantings" (here) , "ravings" (here) or "wild imagining" (everywhere else here or there)). When I asked Google for a possible translation (here), I got "lucubration" which means "composed by night; that which is produced by meditation in retirement". Well, it was not initially what I intended to say, but since I usually work on my blog during the night, when I got awake by one of the girls, I decided to keep this word.... At least, I learnt something today, appart for the code mentioned above....

Friday, October 15 2010

Margin of error, and comparing proportions in the same sample

I recently tried to answer a simple question, asked by @adelaigue. Actually, I thought that the answer would be obvious... but it is a little bit more compexe than what I thought. In a recent pool about elections in Brazil, it was mentionned in a French newspapper that "Mme Rousseff, 62 ans, de 46,8% des intentions de vote et José Serra, 68 ans, de 42,7%" (i.e. proportions obtained from the survey). It is also mentioned that "la marge d'erreur du sondage est de 2,2% " i.e. the margin of error is 2.2%, which means (for the journalist) that there is a "grande probabilité que les 2 candidats soient à égalité" (there is a "large probability" to have equal proportions).
Usually, in sampling theory, we look at the margin of error of a single proportion. The idea is that the variance of \widehat{p}, obtained from a sample of size is
thus, the standard error is
The standard 95% confidence interval, derived from a Gaussian approximation of the binomial distribution is
The largest value is obtained when p is 1/2, and then we have a worst case confidence interval (an upper bound) which is
So with a margin of error means that Hence, with a 5% margin of error, it means that n=400. While 2.2% means that n=2000:
> 1/.022^2
[1] 2066.116   

Classically, we compare proportions between two samples: surveys at two different dates, surveys in different regions, surveys paid by two different newpapers, etc. But here, we wish to compare proportions within the same sample. This has been consider in an "old" paper published in 1993 in the American Statistician,

It contains nice figures to illustrate the difference between the standard approach,

and the one we would like to study here.

This point is mentioned in the book by Kish, survey sampling (thanks Benoit for the reference),

Let and denote empirical frequencies we have obtained from the sample, based on observations. Then since
we have
Thus, a natural margin of error on the difference between the two proportion is here
which is here 4 points
> n=2000
> p1=46.8/100
> p2=42.7/100
> 1.96*sqrt((p1+p2)-(p1-p2)^2)/sqrt(n)
[1] 0.04142327
Which is exactly the difference we have here ! Hence, the probability of reaching such a value is quite small (2%)
> s=sqrt(p1*(1-p1)/n+p2*(1-p2)/n+2*p1*p2/n)
> (p1-p2)/s
[1] 1.939972
> 1-pnorm(p1-p2,mean=0,sd=sqrt((p1+p2)-(p1-p2)^2)/sqrt(n))
[1] 0.02619152

Actually, we can compare the three margin of errors we have so far,
  • the upper bound
  • the "average one"
  • the more accurate one we just obtained,
> p=seq(0,.5,by=.01)
> ic1=rep(1.96/sqrt(4*n),length(p))
> ic2=1.96*sqrt(p*(1-p))/sqrt(n)
> delta=.01
> ic31=1.96*sqrt(2*p-delta^2)/sqrt(n)
> delta=.2
> ic32=1.96*sqrt(2*p-delta^2)/sqrt(n)
> plot(p,ic32,type="l",col="blue")
> lines(p,ic31,col="red")
> lines(p,ic2)
> lines(p,ic1,lty=2)
So on the graph below, the dotted line is the standard upper bound, the plain line in black being a more accurate one when the probability is (the x-axis). The red line is the true margin of error with a large difference between candidates (20 points) and the blue line with a small difference (1 point).

Remark: an alternative is to consider a chi-square test, comparering two multinomial distributions, with probabilities and where is the average proportion, i.e. 44.75%. Then
> p=(p1+p2)/2
> (x2=n*((p1-p)^2/p+(p2-p)^2/p))
[1] 3.756425
> 1-pchisq(x2,df=1)
[1] 0.05260495
Under the null hypothesis, should have a chi-square distribution, with one degree of freedom (since the average is fixed here). Here the probability to reach that level is around 5% (which can be compared with the 2% we add before).

So finally, I would think that here, stating that there is a "large probability" is not correct....

Thursday, October 7 2010

Studying joint effects in a regression

We've seen in the previous post (here)  how important the *-cartesian product to model joint effected in the regression. Consider the case of two explanatory variates, one continuous (, the age of the driver) and one qualitative (, gasoline versus diesel).

  • The additive model
Assume here that
Then, given (the exposure, assumed to be constant) and
Thus, there is a multplicative effect of the qualitative variate.
> reg=glm(nbre~bs(ageconducteur)+carburant+offset(exposition),
+     data=sinistres,family="poisson")
> ageD=data.frame(ageconducteur=seq(17,90),carburant="D",exposition=1)
> ageE=data.frame(ageconducteur=seq(17,90),carburant="E",exposition=1)
> yD=predict(reg,newdata=ageD,type="response")
> yE=predict(reg,newdata=ageE,type="response")
> lines(ageD$ageconducteur,yD,col="blue",lwd=2)
> lines(ageE$ageconducteur,yE,col="red",lwd=2)

On the graph below, we can see that the ratio
is constant (and independent of the age
> plot(ageD$ageconducteur,yD/yE)

  • The nonadditive model
In order to take into accound a more complex (non constant) interaction between the two explanatory variates, consider the following product model,
 > reg=glm(nbre~bs(ageconducteur)*carburant+offset(exposition),
+     data=sinistres,family="poisson")
> ageD=data.frame(ageconducteur=seq(17,90),carburant="D",exposition=1)
> ageE=data.frame(ageconducteur=seq(17,90),carburant="E",exposition=1)
> yD=predict(reg,newdata=ageD,type="response")
> yE=predict(reg,newdata=ageE,type="response")
> lines(ageD$ageconducteur,yD,col="blue",lwd=2)
> lines(ageE$ageconducteur,yE,col="red",lwd=2)

Here, the ratio
is not constant any longer,

  • Mixing additive and nonadditive
It is also possible to consider a model in between: we believe that there is no interaction for young people (say), while there is for older ones. Assume that the beak occurs at age 50,
> reg=glm(nbre~bs(ageconducteur*(ageconducteur<50))+
+     bs(ageconducteur*(ageconducteur>=50))*carburant+offset(exposition),
+     data=sinistres,family="poisson")
> ageD=data.frame(ageconducteur=seq(17,90,by=.1),carburant="D",exposition=1)
> ageE=data.frame(ageconducteur=seq(17,90,by=.1),carburant="E",exposition=1)
> yD=predict(reg,newdata=ageD,type="response")
> yE=predict(reg,newdata=ageE,type="response")
> lines(ageD$ageconducteur,yD,col="blue",lwd=2)
> lines(ageE$ageconducteur,yE,col="red",lwd=2)

Here, the ratio
is constant for young people, while it will change for older ones,

Regression and age class structure

Following my lecture last week on continuous variables in insurance ratemaking (and regression models), a student of mine asked me how  to run (in R) a regression on classes. Consider the following dataset,
>  sinistreUdM <- read.table("",
+  header=TRUE,sep=";")
>  sinistres=sinistreUdM[sinistreUdM$garantie=="1RC",]
>  contratUdM <- read.table("",
+  header=TRUE,sep=";")
>  T=table(sinistres$nocontrat)
>  T1=as.numeric(names(T))
>  T2=as.numeric(T)
>  nombre1 = data.frame(nocontrat=T1,nbre=T2)
>  I = contratUdM$nocontrat%in%T1
>  T1=contratUdM$nocontrat[I==FALSE]
>  nombre2 = data.frame(nocontrat=T1,nbre=0)
>  nombre=rbind(nombre1,nombre2)
>  sinistres = merge(contratUdM,nombre)
Assume that we have the following graph,
> library(tree)
> library(splines)
> arbre=tree((nbre>0)~ageconducteur,data=sinistres,split="gini",mincut = 1)
> age=data.frame(ageconducteur=18:90)
> y=predict(arbre,age)
> plot(age$ageconducteur,y,ylim=c(0,.1),xlab="Age of the driver",
+ ylab="Probability to have - at least - a claim")
> library(splines)
> reg=glm((nbre>0)~bs(ageconducteur,7),data=sinistres,family="binomial")
> y=predict(reg,age,type="response")
> lines(age$ageconducteur,y,col="blue")
> arbre=tree((nbre>0)~ageconducteur,data=sinistres,split="gini",mincut = 6000)
> y=predict(arbre,age)
> lines(age$ageconducteur,y,col="red",type="s",lwd=2)

Dots are empirical average probabilities per age, and here I ran a tree regression (to see which classes we might have, in red) and a logistic regression, with splines (in blue). The only way I found to get values of breaks it to use
> knot=substr(as.character((arbre$frame)$splits),2,
+ nchar(as.character((arbre$frame)$splits)))
> knot=unique(as.numeric(knot))
> abline(v=knot,lty=2,col="grey")
(there is probably a better way, but anyway...)
  • A step function as conditional expected value
With our tree regression, we have constant values for each class. Actually, we have the following predicted values,

or with numerical values,
> arbre
node), split, n, deviance, yval
      * denotes terminal node
 1) root 50000 2366.0 0.04980 
   2) ageconducteur < 29.5 6535  350.8 0.05692 *
   3) ageconducteur > 29.5 43465 2015.0 0.04873 
     6) ageconducteur < 41.5 14832  655.2 0.04632 
      12) ageconducteur < 36.5 8468  377.5 0.04676 *
      13) ageconducteur > 36.5 6364  277.7 0.04573 *
     7) ageconducteur > 41.5 28633 1359.0 0.04998 
      14) ageconducteur < 53.5 14837  734.5 0.05223 
        28) ageconducteur < 48.5 8574  418.3 0.05143 *
        29) ageconducteur > 48.5 6263  316.2 0.05333 *
      15) ageconducteur > 53.5 13796  624.8 0.04755 
        30) ageconducteur < 62.5 7411  312.6 0.04412 *
        31) ageconducteur > 62.5 6385  312.0 0.05153 *

The analogous in a glm regression is to run
> (brk=c(17,sort(knot),120))
[1]  17.0  29.5  36.5  41.5  48.5  53.5  62.5 120.0
> reg=glm((nbre>0)~cut(ageconducteur,brk),data=sinistres,family="binomial")
> y=predict(reg,age,type="response")
> y=predict(reg,age,type="response")
> lines(age$ageconducteur,y,col="purple",lwd=2)

So using function cut() we can generate classes, and predict probabilities to have claims per age class. But here, our prediction is constant per class, which might not be appropriate, as suggested in the regression with splines.
  • A piecewise "linear" function as conditional expected value
An alternative is to try to have a piecewise linear function (linear, at least, in the score). A first idea could be to run a regression per subclass
> (brk=c(17,sort(knot),120))
[1]  17.0  29.5  36.5  41.5  48.5  53.5  62.5 120.0
> classes=cut(sinistres$ageconducteur,brk)
> age=data.frame(ageconducteur=18:90)
> for(i in 1:length(levels(classes))){
+ reg=glm((nbre>0)~ageconducteur,data=sinistres,family="binomial",
+     subset=(classes==levels(classes)[i]))
+ I=((age>brk[i])& age<=brk[i+1])==TRUE
+ y=predict(reg,newdata=data.frame(ageconducteur=age[I]),type="response")
+ lines(age[I],y,col="purple",lwd=2)
+ }

An alternative will be to use the * cartesian product, between the age, and the class   
> reg=glm((nbre>0)~ageconducteur*cut(ageconducteur,brk),data=sinistres,family="binomial")
> reg
Call:  glm(formula = (nbre > 0) ~ ageconducteur * cut(ageconducteur,      brk), family = "binomial", data = sinistres)
              cut(ageconducteur, brk)(29.5,36.5] 
              cut(ageconducteur, brk)(36.5,41.5] 
              cut(ageconducteur, brk)(41.5,48.5] 
              cut(ageconducteur, brk)(48.5,53.5] 
              cut(ageconducteur, brk)(53.5,62.5] 
               cut(ageconducteur, brk)(62.5,120] 
ageconducteur:cut(ageconducteur, brk)(29.5,36.5] 
ageconducteur:cut(ageconducteur, brk)(36.5,41.5] 
ageconducteur:cut(ageconducteur, brk)(41.5,48.5] 
ageconducteur:cut(ageconducteur, brk)(48.5,53.5] 
ageconducteur:cut(ageconducteur, brk)(53.5,62.5] 
 ageconducteur:cut(ageconducteur, brk)(62.5,120] 
Degrees of Freedom: 49999 Total (i.e. Null);  49986 Residual
Null Deviance:      19790
Residual Deviance: 19760        AIC: 19790
> y=predict(reg,newdata=data.frame(ageconducteur=age),type="response")
> lines(age,y,col="purple",lwd=2)

Finally, note that it is also possible to use library(segmented)
> library(segmented)
> Y=(sinistres$nbre>0)
> X=sinistres$ageconducteur
> reg=glm(Y~X,family="binomial")
> slopes=segmented(reg,seg.Z=~X,psi=list(X=c(30,60)))
> slopes
Call: segmented.glm(obj = reg, seg.Z = ~X, psi = list(X = c(30, 60)))
Meaningful coefficients of the linear terms:
(Intercept)            X         U1.X         U2.X 
   -1.11619     -0.06738      0.06846      0.01084 
Estimated Break-Point(s) psi1.X psi2.X : 27.88 74.00
Degrees of Freedom: 49999 Total (i.e. Null);  49994 Residual
Null Deviance:     19790
Residual Deviance: 19770      AIC: 19780
> slope(slopes)
            Est.  St.Err. t value CI(95%).l CI(95%).u
slope1 -0.067380 0.020880 -3.2270 -0.108300 -0.026460
slope2  0.001081 0.002232  0.4843 -0.003294  0.005456
slope3  0.011920 0.015160  0.7863 -0.017800  0.041640

Monday, October 4 2010

Statistique de l'assurance STT6705V, polycopié (1)

Les notes de cours sont en ligne ici, avec les 2 premières parties évoquées en cours, à savoir la tarification du risque de masse (a priori), et le provisionnement pour sinistres à payer (que l'on commencera mercredi). La dernière partie sur la construction des tables de mortalité prospectives arrivera plus tard (je suis en train de rajouter une section d'assurance vie permettant de rappeler des notations qui sont utilisées par la suite). Il s'agit d'ébauches de notes, toutes les remarques sont les bienvenues.
Sinon pour revoir le dernier cours, c'est ici et (les liens ont été réparés pour les séances passées).

Friday, October 1 2010

Too large datasets for regression ? What about subsampling....

(almost) recently, a classmate working in an insurance company told me he had too large datasets to run simple regressions (GLM, which involves optimization issues), and that they were thinking of a reward for the one who will write the best R-code (at least the fastest). My first idea was to use subsampling techniques, saying that 10 regressions on 100,000 observations can take less time than a regression on 1,000,000 observations. And perhaps provide also better results...

  • Time to run a regression, as a function of the number of observations
Here, I generate a dataset as follows
and we fit
where is a spline function (just to make it as general as possible, since in insurance ratemaking, we include continuous variates that do not influence claims frequency linearly in the score). Yes, there might be also useless variables, including one of them which is strongly correlated with one that has an impact in the regression. The code to generate the dataset is simply
> n=10000
> X1=rexp(n)
> X2=sample(c("A","B","C"),size=n,replace=TRUE)
> X3=runif(n)
> Z=rmnorm(n,c(0,0),matrix(c(1,0.8,.8,1),2,2))
> X4=Z[,1]
> X5=Z[,2]
> X6=X1^2
> E=runif(n)
> lambda=.2*X5-4*dbeta(X3,2,5)+X1+
> Y=rpois(n,exp(lambda))
> base=data.frame(Y,X1,X2,X3,X4,X5,X6,E)
We would like the study the time it takes to run a regression, as a function of the size (i.e. the number of lines of the dataset. 
> system.time( glm(Y~bs(X1)+X2+X3+X4+
+ X5+X6+offset(log(E)),family=poisson,
+ data=base) )
utilisateur     système      écoulé
0.25        0.00        0.25 
Here, the time I look at is the last one. But so far, it was rather simple, but it is not the best model I can get. Let us use a stepwise (backward) variable selection,
> system.time( step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=poisson,
+ data=base)) )
Start:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))
Step:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + offset(log(E))
Df Deviance    AIC
<none>        2236.0 2882.1
- X5      1   2240.1 2884.2
- X4      1   2244.1 2888.2
- X3      1   4783.2 5427.3
- X2      2   5311.4 5953.5
- bs(X1)  3   6273.7 6913.8
utilisateur     système      écoulé
1.82        0.03        1.86 
Finally, from the first regression, we have points in black (based on 200 simulated datasets), and with a stepwise procedure, we have the points in red.

i.e. it might look linear (proportional), but if it was linear, then on a log-log scale, we should have also straigh lines, with slope 1,

Actually, it looks like a convex function.

The interpretation of that convexity might lead to misinterpretation. On the graph below on the left, on a dataset two times bigger than the previous one (black point) will be less than two times longer to run, while on the right, it  will be more than two timess longer,

Convexity can simply be interpreted as "too large datasets take time, and too small too...". Which is a first step: it should be interesting, in some cases, to run several regressions on smaller datasets....
  • Running 100 regressions on 100 lines, or running 1 regression on 10,000 lines ?
Here, we have datasets with,000 lines. The questions is how long will it take if we subdived into subsamples (of equal size), and run regressions ?
> nk=trunc(n/k)rep(1:k,each=nk); nt=nk*k
> base=data.frame(Y[1:nt],X1[1:nt],
+ X2[1:nt],X3[1:nt],X4[1:nt],X5[1:nt],
+ X6[1:nt],E[1:nt],classe)
> system.time( for(j in 1:k){
+  glm(Y~bs(X1)+X2+X3+X4+X5+
+ X6+offset(log(E)),family=poisson
+ ,data=base,subset=classe==j) })
utilisateur     système      écoulé
1.31        0.00        1.31
> system.time( for(j in 1:k){
+      step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=
+ poisson,data=base,subset=classe==j)) })
Start:  AIC=183.97
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))
  Df Deviance    AIC
<none>        117.15 213.04
- X2      2   250.15 342.04
- X3      1   251.00 344.89
- X4      1   420.63 514.53
- bs(X1)  3   626.84 716.74
utilisateur     système      écoulé
11.97        0.03       12.31 
On the graph below, we have the time (y-axis, here on a log scale) it took to run regression on samples of size, as function of (x-axis), including the time it took to run the regression on a dataset of size which is the concentration of dots on the left (i.e., both on the 6 regressors - in black - and with a strepwise procedure - in red. One has to keep in mind that I did not remove the printing option in the stepwise procedure, so it might be difficult to compare the two clouds (black vs. red). Nevertheless, we clearly see that if we run regression on samples of size, when is not too large, i.e. less than 10 or 15, it is not longer than the regression on,000 lines.

So here we see that running 100 regressions on 2,000 lines is longer than running 1 regression on 200,000 lines... But maybe we are not comparing things that are actually comparable: what if it takes a bit longer, but we strongely improve the quality of our estimators ?
  • What about the quality of the output ?
Here, we consider only one dataset, with,000 lines (just to make it run a bit faster). And subsets. Recall that the generated dataset is from
and we fit
Here, we plot here and a confidence interval, defined as
The lightblue segment is the initial estimator, while the blue one is obtained from the stepwise procedure. The grey area represent the estimation on the overall sample, while the segments on the right are the estimators (each on samples of size

We can see that we have much more volatility on those estimators, but the average (horizontal doted lines) are not so bad... The true value (i.e. the one used to generate the dataset is the dotter black horizontal line).
And if we repeat that on 1,000 simulated dataset, we obtaind the following distribution for (blue line), so we have an unbiased estimator of our parameter (the verticular line being here the true value), here including a stepwise procedure,

But if we add the the red curve is the average of the the previous one being now the clear blue line in the back, we see that taking average of estimators on subsamples is not bad at all, on the contrary,

and for those who think that the stepwise procedure is a mistake, here is what we get without it,

So what we can see is that running 20 regressions can take (a little) more time (from what we've seen earlier) than running only one on the whole dataset.... but it provides better estimates. So the tradeoff is not that simple, and maybe running severall regressions on huge datasets can be a proper alternative.

Wednesday, September 29 2010

Detecting distributions with infinite mean

In a post I published a few month ago (in French, here, based on some old paper, there), I mentioned a statistical procedure to test if the underlying distribution of an i.i.d. sample had a finite mean (based on extreme value results). Since I just used it on a small dataset (yes, with real data), I decided to post the R code, since it is rather simple to use. But instead of working on that dataset, let us see what happens on simulated samples. Consider observations generated from a Pareto distribution
(here, as a start)
> a=2
> X=runif(200)^(-1/a)
Here, we will use the package developped by Mathieu Ribatet,
> library(RFA)
Le chargement a nécessité le package : tcltk
Chargement de Tcl/Tk... terminé
Message d'avis :
le package 'RFA' a été compilé avec la version R 2.10.1
  • Using Generalized Pareto Distribution (and LR test)
A first idea is to fit a GPD distribution on observations that exceed a threshold >1.
Since we would like to test (against the assumption that the expected value is finite, i.e., it is natural to consider the likelihood ratio
Under the null hypothesis, the distribution of should be chi square distribution with one degree of freedom. As mentioned here, the significance level is attained with a higher accuracy by employing Bartlett correction (there). But let  us make it as simple as possible for the blog, and use the chi-square distribution to derive the p-value.
Since it is rather difficult to select an appropriate threshold, it can be natural (as in Hill's estimator) to consider, and thus, to fit a GPD on the largest values. And then to plot all that on a graph (like the Hill plot)
> Xs=rev(sort(X))
> s=0;G=rep(NA,length(Xs)-14);Gsd=G;LR=G;pLR=G
> for(i in length(X):15){
+ s=s+1
+ FG=fitgpd(X,Xs[i],method="mle")
+ FGc=fitgpd(X,Xs[i],method="mle",shape=1)
+ G[s]=FG$estimate[2]
+ Gsd[s]=FG$std.err[2]
+ FGc$fixed
+ LR[s]=FGc$deviance-FG$deviance
+ pLR[s]=1-pchisq(LR[s],df=1)
+ }
Here we keep the estimated value of the tail index, and the associated standard deviation of the estimator, to draw some confidence interval (assuming that the maximum likelihood estimate is Gaussian, which is correct only when n is extremely large). We also calculate the deviance of the model, the deviance of the constrained model (, and the difference, which is the log likelihood ratio. Then we calculate the p-value (since under the likelihood ratio statistics has a chi-square distribution).
If, we have the following graph, with on top the p-value (which is almost null here), the estimation of the tail index he largest values (and a confidence interval for the estimator),
If (finite mean, but infinite variance), we have
i.e. for those two models, we clearly reject the assumption of infinite mean (even if gets too close from 1, we should consider thresholds large enough). On the other hand, if (i.e. infinite mean) we clearly accept the assumption of infinite mean (whatever the threshold),
  • Using Hill's estimator
An alternative could be to use Hill's estimator (with Alexander McNeil's package). See here for more details on that estimator. The test is simply based on the confidence interval derived from the (asymptotic) normal distribution of Hill's estimator,
> library(evir)
> Xs=rev(sort(X))
> HILL=hill(X)
> G=rev(HILL$y)
> Gsd=rev(G/sqrt(HILL$x))
> pLR=1-pnorm(rep(1,length(G)),mean=G,sd=Gsd)
Again, if, we clearly rejct the assumption of infinite mean,
and similarly, if (finite mean, but infinite variance)
Here the test is more robust than the one based on the GPD. And if (i.e. infinite mean), again we accept,
Note that if, it is still possible to run the procedure, and hopefully, it suggests that the underlying distribution has infinite mean,
(here without any doubt). Now you need to wait a few days to see some practical applications of the idea (there was on in the paper mentioned above actually, on business interruption insurance losses).

- page 1 of 2