Tag Archives: ACT2040

Segmentation en tarification, compléments

Dans le premier cours d’actuariat IARD, nous avons vu l’importance de la ségmentation, et son implication sur le calcul des primes (passer d’une espérance mathématique à une espérance conditionnelle). Pour aller un peu plus loin, quelques compléments,

pour une lecture plus économique de la problématique de la segmentation en assurance

ou pour une lecture plus légale

Sinon, plusieurs articles de vulgarisation peuvent être lu sur internet,

La première démo aura lieu lundi, en salle informatique. Karim sera une introduction au langage R, à la manipulation des variables (qualitatives et quantitatives). Je mettrais en ligne les transparents en fin de semaine, et les codes seront mis en ligne dans le courant de la semaine prochaine.

Comme annoncé hier, il n’y aura pas cours mercredi prochain. Le mercredi suivant, nous verrons la modélisation des variables indicatrices, i.e. la régression logistique, et les arbres de régression. On supposera que le modèle linéaire aura été vu, je mets un lien vers les transparents du cours ACT6420 de la session passée, notes de cours transparents1 et transparents2. Il est aussi possible de relire Frees (2010), chapitres 3, 4, 5 et 6.

Pour commencer à pratiquer la régression logistique, on utilisera la petite fonction suivante

logit = function(formula, lien="logit", data=NULL) {
glm(formula,family=binomial(link=lien),data)
}

Sinon, la Casualty Actuarial Society a mis en ligne plusieurs documents en ligne sur les arbres de régression (qui sont peu abordés dans les livres mentionnés auparavant),

pour une comparaison de toutes les méthodes

Les transparents seront mis en ligne en fin de semaine prochaine. A suivre donc…

Actuariat IARD

Cet hiver (même si la nouvelle ne sera officielle qu’à la rentrée), je devrais donner le cours ACT2040, actuariat IARD. Le plan de cours sera bientôt en ligne, mais je peux déjà dire que le cours sera basé sur le Tome 2 du livre écrit avec Michel Denuit il y a quelques années, mathématiques de l’assurance non-vie. Le cours est une suite du cours ACT6420 méthodes de prévisions, donné cet automne (qui est un prérequis indiqué sur le site du registrariat http://websysinfo.uqam.ca/…): je partirais donc du fait que le modèle linéaire de régression est connu (et compris) et que tout le monde sait utiliser R, et lire des sorties de régression. Mais les premières démonstrations reviendront sur l’utilisation de R, et sur l’analyse de la variance, que l’on n’a pas vraiment eu le temps d’aborder dans le cours de régression. Pour des références sur R, je conseille

  • “R pour les débutants” d’Emmanuel Paradis, (PDF)
  • “Introduction à la programmation en S” par Vincent Goulet, (PDF)

pour les documents en français, ou pour des documents plus complets, mais en anglais

  • “R for Beginners” d’Emmanuel Paradis (PDF),
  • “An Introduction to R” par Longhow Lam (PDF)
  • “The R language — a short companion” par Marc Vandemeulebroecke (PDF),
  • “The R Guide” par Jason Owen (PDF),
  • “Econometrics in R” par Grant Farnsworth (PDF) pour aller plus loin sur les régressions,
  • “Practical Regression and Anova using R” by Julian Faraway (PDF) sur le meme sujet
  • “Statistics with R and S-Plus” d’Hugo Quené (PDF)
  • “Statistical Computing and Graphics Course Notes” par Frank Harrell, (PDF).
  • “Using R for Data Analysis and Graphics – Introduction, Examples and Commentary” par John Maindonald (PDF).

Sinon, les transparents du premier cours sont en ligne ici, et le plan de cours est

et je mettrais bientôt en ligne des liens vers des bases de données que l’on utilisera tout au long du cours, ou en démonstration.

Pour les références, je citerais deux livres sur lesquels je m’appuierai beaucoup car je les connais presque par cœur. Ils sont disponibles à la Coopuqam

ACT2040: examen final (suite et fin)

Comme annoncé, l’examen terminal a lieu ce matin, et il est basé sur des données mises en ligne sur le blog depuis une dizaine de jours. L’énoncé est un peu longue, avec une bonne trentaine de pages de sorties informatiques, et un énoncécomprenant une trentaine de questions. Des éléments de correction sont en ligne ci-dessous (plutôt que de taper une correction en pdf, autant mettre le code en ligne).

On suppose importée les bases (sinon, le code est ici). Pour la première partie – sur la tarification (a priori) – la fréquence annuelle de sinistres est tout simplement

> mean(BASEN$nombre)/mean(BASEN$exposition)
[1] 0.05241692
> exp(glm(nombre~1,offset=log(exposition),data=BASEN,
+ family=poisson(link="log"))$coefficients)
(Intercept)
0.05241692

(je vais essayer de donner à chaque question deux méthodes pour obtenir le résultat). Si on suppose que le nombre de sinistres par an suit une loi de Poisson, alors la probabilité d’avoir au moins un accident est

> lambda=mean(BASEN$nombre)/mean(BASEN$exposition)
> 1-ppois(0,lambda)
[1] 0.05106684
> 1-exp(-lambda)
[1] 0.05106684

Si on cherche à distinguer par sexe du conducteur principal

> mean(BASEN$nombre[BASEN$sexeconducteur=="H"])/
+ mean(BASEN$exposition[BASEN$sexeconducteur=="H"])
[1] 0.05218143
> mean(BASEN$nombre[BASEN$sexeconducteur=="F"])/
+ mean(BASEN$exposition[BASEN$sexeconducteur=="F"])
[1] 0.05457013
> exp(glm(nombre~0+sexeconducteur,offset=log(exposition),
+ data=BASEN,family=poisson(link="log"))$coefficients)
sexeconducteurF sexeconducteurH
0.05457013      0.05218143

Ces deux valeurs ne sont pas significativement différentes. Si on regarde la sortie de la régression

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
sexeconducteurF -2.90827    0.05464  -53.23   <2e-16 ***
sexeconducteurH -2.95303    0.01848 -159.82   <2e-16 ***

les deux coefficients sont quasiment identiques (en enlevant la constante). Si on rajoute la constante, la modalité de référence devientle sexe féminin

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -2.90827    0.05464 -53.230   <2e-16 ***
sexeconducteurH -0.04476    0.05768  -0.776    0.438

et la modalité homme n’est alors plus significative. Bref, le sexe du conducteur n’est pas une variable explicative.

Ensuite, parmi les variables qui semble pouvoir etre utilisée pour modéliser le nombre de sinistre, on pourra retenir l’age du conducteur (décroissant jusqu’à 30 ans, puis croissant), l’age du véhicule (décroissant), l’ancienneté du permis de conduire (décroissant les 10 premières années puis croissant – ou indépendant), éventuellement certains type de véhicules (le coupé cabriolet semble avoir – en moyenne – plus d’accident que la Jeep), la situation familiale, le lieu d’habitation (en particulier rural semble avoir moins d’accident qu’urbain), le statut d’habitation, le type de paiement de la prime (mensuel semblant plus risqué qu’annuel), le poids du véhicule, ou la marque (ou plutot certaines marques, comme Peugeot versus Burstner Mobil). Bref, toutes, sauf le sexe du conducteur.

Pour un conducteur de 45 ans qui a son permis de conduire depuis seulement 5 ans, les trois modèles donnent la fréquence annuelle suivante

> nbase=data.frame(ageconducteur=45,agepermis=5,exposition=1)
> exp(t(c(1,5))%*%reg6$coefficients)
[,1]
[1,] 0.03395972
> predict(reg6,newdata=nbase,type="response")
1
0.03395972
> exp(t(c(1,45))%*%reg7$coefficients)
[,1]
[1,] 0.04100356
> predict(reg7,newdata=nbase,type="response")
1
0.04100356
> exp(t(c(1,5,45))%*%reg8$coefficients)
[,1]
[1,] 0.05395398
> predict(reg8,newdata=nbase,type="response")
1
0.05395398

Les trois modèles sont respectivement un modèle de Poison avec un lien log, un modèle binomiale négative avec un lien log, et un modèle binomial, avec un lien logit. Les deux premiers sont très proches: ils ont la meme fonction de lien, et comme nous l’avions vu à l’examen intra, le modèle binomiale négative est un modèle de Poisson avec une variable latente non observée qui suivrait une loi Gamma. Donc pour ces deux modèles, si les classes sont homogènes, et que par classe, la loi de Poisson est adaptée, il est normal d’avoir des coefficients très proches. Pour le modèle binomial avec un lien logistique, il faut se souvenir que la fonction lien est

http://freakonometrics.blog.free.fr/public/perso4/act2040-1.gif

car ici la probabilité d’avoir un accident est faible. Donc les deux fonctions liens sont proches. De plus, pour une loi de Poisson

http://freakonometrics.blog.free.fr/public/perso4/act2040-2.gif
http://freakonometrics.blog.free.fr/public/perso4/act2040-3.gif

Bref, pour notre loi de Poisson, on a presque une loi binomiale de paramètre http://freakonometrics.blog.free.fr/public/perso4/act2040-4.gif si http://freakonometrics.blog.free.fr/public/perso4/act2040-4.gif est petit.

Pour la prédiction suivant les trois modèles, on utilise classiquement

> exp(t(c(1,40,0,0,1,0,1,0))%*%reg9$coefficients)
           [,1]
[1,] 0.04351037

Comme on a accès aux sorties, pour les trois modèles, les trois prédictions sont très proches (mais on s’en doutait d’après la question précédante)

> predict(reg9,newdata=nbase,type="response")
         1 
0.04351037 
> predict(reg10,newdata=nbase,type="response")
         1 
0.04351056 
> predict(reg11,newdata=nbase,type="response")
         1 
0.04352741

Pour la méthode de biais minimal, il s’agit du principe de balancement, qui coïncide avec la régression de Poisson sur les classes. Pour les régressions stepwise, dans un cas, on élimine (étape par étape) des variables, et dans l’autre cas, des modalités de variables (pour les variables qualitatives).

Les deux modèles semblent valides, au sens où toutes les variables sont significatives. Le découpage en classes pour l’age du véhicule semble justifié par l’arbre de régression, qui suggèrait de découper à 8 et 14 ans. Pour les prédictions, on a pour le premier modèle

> nbase15=data.frame(ageconducteur=40,agepermis=20,
+ habitation="urbain",locataire=TRUE,paymensuel=FALSE,
+ marie=TRUEFORD=TRUE,RENAULT=FALSE,PEUGEOT=FALSE,
+ exposition=1)
> predict(reg15,newdata=nbase15,type="response")
         1 
0.03398312

alors que le second donne

> nbase16=data.frame(ageconducteur=40,agevehF="NEUF",rural=FALSE,
+ locataire=TRUE,paymensuel=FALSE,marie=TRUE,
+ jeune=FALSE,ancien=FALSE,exposition=1)
> predict(reg16,newdata=nbase16,type="response")
         1 
0.05208656

Moralité, nos deux modèles (qui semblent autant valides l’un que l’autre) donnent des prédictions relativement différentes. Si l’assuré change de voiture, cela n’a pas d’influence pour le second modèle (le véhicule est toujours dans la catégorie neuf). Pour le premier, il perd le rabais qu’il avait en conduisant un véhicule Ford, et donc la fréquence annuelle de sinistre augmente de 23.5%

> 1/exp(reg15$coefficients["FORDTRUE"])
FORDTRUE 
1.235329

Bref, la prédiction devient

> predict(reg15,newdata=nbase15,type="response")/
+ exp(reg15$coefficients["FORDTRUE"])
         1 
0.04198033

On peut d’ailleurs le vérifier simplement,

> nbase15b=data.frame(ageconducteur=40,agepermis=20,
+ habitation="urbain",locataire=TRUE,paymensuel=FALSE,
+ marie=TRUE,FORD=FALSE,RENAULT=FALSE,PEUGEOT=FALSE,
+ exposition=1)
> predict(reg15,newdata=nbase15b,type="response")
         1 
0.04198033

Si on passe maintenant à la modélisation de la charge de sinistres, le coût moyen par sinistre serait

> mean(BASEY$cout)
[1] 1468.522

Cette fois, contrairement à la modélisation de la fréquence, il est plus difficile de trouver des variables explicatives. A priori, au vu des graphiques, aucune variable ne semblerait convenir…
Si on regarde la prédiction pour le premier modèle obtenu,

> nbase7=data.frame(VOLKSWAGEN=FALSE,agepermis=20,
+ sexeconducteur="H",usage="TOUS_DEPLACEMENTS")
> predict(regr7,newdata=nbase7,type="response")
       1 
1388.622

alors que pour le second,

> exp(predict(regr8,newdata=nbase7))*
+ exp(summary(regr8)$sigma^2/2)
       1 
1359.552

Si c’est la femme de l’assuré qui avait pris l’assurance, cela n’aurait rien changé pour le second modèle. En revanche, pour le premier, l’assuré perdrait le rabais qu’il avait en étant un homme, et donc le cout moyen attendu augmenterait de 23%,

> 1/exp(regr7$coefficients["sexeconducteurH"])
sexeconducteurH 
       1.237832

i.e.

> predict(regr7,newdata=nbase7,type="response")/
+ exp(regr7$coefficients["sexeconducteurH"])
       1 
1718.88

La encore, on retrouve le meme résultat directement

> nbase7b=data.frame(VOLKSWAGEN=FALSE,agepermis=20,
+ sexeconducteur="F",usage="TOUS_DEPLACEMENTS")
> predict(regr7,newdata=nbase7b,type="response")
       1 
1718.88

Passons maintenant aux calculs de la prime pure. Sans segmentation, la prime pure serait

> sum(BASEN$nombre)/sum(BASEN$exposition)*
+ mean(BASEY$cout)
[1] 76.9754

Les prédictions des modèles pour la fréquence donnaient

> (n1=predict(reg15,newdata=nbase15,type="response"))
         1 
0.03398312 
> (n2=predict(reg16,newdata=nbase16,type="response"))
         1 
0.05208656

alors que pour les couts moyens par sinistre, nous avions

> (c1=predict(regr7,newdata=nbase7,type="response"))
       1 
1388.622 
> (c2=exp(predict(regr8,newdata=nbase7))
+ *exp(summary(regr8)$sigma^2/2))
       1 
1359.552

Les quatre primes pures obtenues sont alors

> c(n1,n2)%*%t(c(c1,c2))
            1        1
[1,] 47.18972 46.20183
[2,] 72.32855 70.81440

Passons maintenant à la seconde partie – sur le provisionnement des sinistres à payer. On commence par travailler sur le triangle

TRIANGLE2010
          D1      D2      D3      D4      D5      D6      D7      D8      D9
A2002 357848 1124788 1735330 2218270 2745596 3319994 3466336 3606286 3833515
A2003 352118 1236139 2170033 3353322 3799067 4120063 4647867 4914039 5339085
A2004 290507 1292306 2218525 3235179 3985995 4132918 4628910 4909315      NA
A2005 310608 1418858 2195047 3757447 4029929 4381982 4588268      NA      NA
A2006 443160 1136350 2128333 2897821 3402672 3873311      NA      NA      NA
A2007 396132 1333217 2180715 2985752 3691712      NA      NA      NA      NA
A2008 440832 1288463 2419861 3483130      NA      NA      NA      NA      NA
A2009 359480 1421128 2864498      NA      NA      NA      NA      NA      NA
A2010 376686 1363294      NA      NA      NA      NA      NA      NA      NA

Pour obtenir le montant payé en 2008 pour les sinistres survenus en 2004, on utilise la différence entre

> TRIANGLE2010["A2004",c("D4","D5")]
           D4      D5
A2004 3235179 3985995

i.e.

> TRIANGLE2010["A2004","D5"]-TRIANGLE2010["A2004","D4"]
[1] 750816

Pour le montant payé au total en 2005, on utilise

> CUMULPAIEMENTS2010=TRIANGLE2010
> INCREMENTPAIEMENTS2010=CUMULPAIEMENTS2010
> INCREMENTPAIEMENTS2010[,2:9]=CUMULPAIEMENTS2010[,2:9]-
+ CUMULPAIEMENTS2010[,1:8]
> sum(diag(as.matrix(INCREMENTPAIEMENTS2010[1:4,4:1])))
[1] 2729241

On attaque maintenant la méthode Chain Ladder. Les coefficients de développement sont

> lambda=rep(NA,8)
> for(k in 1:8){
+ lambda[k]=(sum(TRIANGLE2010[1:(9-k),k+1])/
+ sum(TRIANGLE2010[1:(9-k),k]))}
> lambda[1:5]
[1] 3.474193 1.704149 1.460866 1.161765 1.095763

On peut alors obtenu ce qui sera payé en 2011 pour les sinistres survenus en 2010,

> TRIANGLE2010["A2010","D1"]*(lambda[1]-1)
[1] 931993.8

alors que pour le montant payé en 2011 pour les sinistres survenus en 2009,

> TRIANGLE2010["A2009","D2"]*(lambda[2]-1)
[1] 1000686

Pour le montant payé en 2010 et 2011 pour les sinistres survenus en 2010,

> TRIANGLE2010["A2010","D1"]*lambda[1]*lambda[2]
[1] 2230186

Pour les sinistres survenus en 2008, la charge total estimée sera

> TRIANGLE2010["A2008","D3"]*prod(lambda[3:8])
[1] 5531128

ce qui fait un montant de provision

> TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1)
[1] 3111267

On rajoute maintenant une diagonale dans notre triangle, i.e. une nouvelle année de paiements (en 2011).
Le montant payé en 2011 pour les sinistres survenus en 2010 est

> TRIANGLE2011["A2010","D2"]-TRIANGLE2011["A2010","D1"]
[1] 986608

à comparer avec

> TRIANGLE2010["A2010","D1"]*(lambda[1]-1)
[1] 931993.8

prédit par la méthode Chain-Ladder, soit une erreur de l’ordre de 5%,

> (TRIANGLE2011["A2010","D2"]-
+ TRIANGLE2011["A2010","D1"])/
+ (TRIANGLE2010["A2010","D1"]*(lambda[1]-1))-1
[1] 0.05859927

Si on compare maintenant les coefficients de transition de la méthode Chain Ladder,

> lambda2011=rep(NA,9)
> for(k in 2:9){
+ lambda2011[k]=(sum(TRIANGLE2011[1:(10-k),k+1])/
+ sum(TRIANGLE2011[1:(10-k),k]))}
> lambda2011[1:5]
[1]       NA 1.747333 1.457413 1.173852 1.103824
> lambda[1:5]
[1] 3.474193 1.704149 1.460866 1.161765 1.095763

Le coefficient permettant de passer de la colonne 2 à la colonne 3 augmente (légèrement). Si on regarde la prédiction du montant de paiements effectué en 2012 pour les sinistres survenus en 2009,

> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)
[1] 1310258

qui sont à comparer avec les

> TRIANGLE2010["A2009","D2"]*lambda[2]*(lambda[3]-1)
[1] 1116132

estimés un an auparavant. La variation entre les estimations est de l’ordre de 17%,

> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)/
+ (TRIANGLE2010["A2009","D2"]*lambda[2]*(lambda[3]-1))-1
[1] 0.1739278

Entre 2010 et 2012, on peut estimer que l’on payera

> TRIANGLE2011["A2010","D2"]*lambda2011[2]
[1] 2382128

pour les sinistres survenus en 2010.
Si on se focalise sur les sinistres survenus en 2008, la prédiction de la charge finale est ici

> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])
[1] 5660771

à comparer avec la prédiction faire en 2010,

> TRIANGLE2010["A2008","D3"]*prod(lambda[3:8])
[1] 5531128

i.e. notre nouvelle prédiction est 2% plus élevée,

> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])/
+ (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))-1
[1] 0.02343867

Le provisions pour sinistres à payer pour les sinistres survenus en 2008 est estimé, fin 2011, à

> TRIANGLE2011["A2008","D4"]*(prod(lambda2011[4:9])-1)
[1] 2177641

alors que fin 2010, un montant de provision avait été constitué, à hauteur de

> TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1)
[1] 3111267

On notera que la différence

> TRIANGLE2011["A2008","D4"]*(prod(lambda2011[4:9])-1)-
+ TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1)
[1] -933626.7

vient des paiements effectués en 2011, et du changement dans l’estimation de charge ultime

> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])-
+ (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))
[1] 129642.3
> TRIANGLE2011["A2008","D4"]-TRIANGLE2011["A2008","D3"]
[1] 1063269
> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])-
+ (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))-
+ (TRIANGLE2011["A2008","D4"]-TRIANGLE2011["A2008","D3"])
[1] -933626.7

On attaque enfin la dernière partie, avec la régression de Poisson. La prédiction du montant de paiements faits en 2012 pour les sinistres survenus en 2009 est ici

> Y=as.vector(as.matrix(INCREMENTPAIEMENTS2011))
> ANNEE=rep(2002:2010,10)
> DEVEL=rep(1:10,each=9)
> baseTriangle=data.frame(Y,A=as.factor(ANNEE),
+ D=as.factor(DEVEL))
> reg=glm(Y~A+D,data=baseTriangle,famil=poisson(link="log"))
> nbase=data.frame(A="2009",D="4")
> predict(reg,newdata=nbase,type="response")
      1 
1310258

qui correspond à ce que donnait la méthode Chain Ladder

> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)
[1] 1310258

(mais on savait que les deux modèles prédisaient la même chose). Pour les sinistres survenus en 2005, le montant de provisions à constituer est

> nbase=data.frame(A="2005",D=c("8","9","10"))
> predict(reg,newdata=nbase,type="response")
       1        2        3 
247190.0 370179.3  92268.5 
> sum(predict(reg,newdata=nbase,type="response"))
[1] 709637.8

qui, là encore, correspond à ce que prédisait la méthode Chain Ladder.

> TRIANGLE2011["A2005","D7"]*(prod(lambda2011[7:9])-1)
[1] 709637.8

Simulations et incréments négatifs

Pour faire suite (rapidement) à mon dernier billet, le cas le plus fréquent pour observer des incréments négatifs est lorsque l’on fait des simulations pour obtenir la distribution du montant de paiements futurs. Rappelons que pour générer un pseudo-triangle, on utilise les résidus de Pearson, et on pose

http://freakonometrics.blog.free.fr/public/perso4/pseudo-triangle.gifOn voit que l’on risque d’avoir un incrément négatif si le résidu est négatif, et plus grand (en valeur absolue) que la racine carrée de la prédiction (par notre modèle Poissonnien, en l’occurrence). Or dans le triangle vu en cours, on a

> source("https://perso.univ-rennes1.fr/
arthur.charpentier/bases.R")
> INC=PAID
> INC[,2:6]=PAID[,2:6]-PAID[,1:5]
> Y=as.vector(INC)
> D=rep(1:6,each=6)
> A=rep(2001:2006,6)
> base=data.frame(Y,D,A)
> reg=glm(Y~as.factor(D)+as.factor(A),
+     data=base,family=poisson(link="log"))
> Yp=predict(reg,type="response",
+ newdata=base)
> erreurs=residuals(reg,"pearson")
> min(sqrt(Yp[is.na(Y)==FALSE]))
[1] 2.868171
> min(erreurs)
[1] -2.213449

autrement dit, on n’aura jamais d’incrément négatif lors de nos boucle, si l’on génère les résidus par bootstrap. Sinon, avec des lois paramétriques non-bornées inférieurement (loi normale, loi de Student), il est tout a fait possible d’obtenir des incréments négatifs. La méthode pour éviter le problème des incréments négatifs dans les boucles… est probablement de ne pas prendre en compte les scénarios où des incréments négatifs ont été obtenus, i.e.

R=rep(NA,10000)
for(s in 1:10000){
serreur=sample(erreurs,
size=36,replace=TRUE)
E=matrix(serreur,6,6)
sY=matrix(Yp,6,6)+E*sqrt(matrix(Yp,6,6))
if(min(sY[is.na(Y)==FALSE])>=0){
sbase=data.frame(sY=as.vector(sY),D,A)
sbase$sY[is.na(Y)==TRUE]=NA
sreg=glm(sY~as.factor(D)+as.factor(A),
data=sbase,family=poisson(link="log"))
sYp=predict(sreg,type="response",
newdata=sbase)
R[s]=sum(sYp[is.na(Y)==TRUE])}
}

Une stratégie un peu plus propre pourrait être de mettre des 0 dès qu’on obtient un incrément négatif… Mais encore une fois, c’est de la cuisine.

Incréments négatifs dans les triangles de paiements

Considérons le triangle d’incréments de payements suivants,

> PAID
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 4372 4411 4428 4435 4456
[2,] 3367 4659 4696 4720 4730   NA
[3,] 3871 5345 5338 5420   NA   NA
[4,] 4239 5917 6020   NA   NA   NA
[5,] 4929 6794   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA
> INC=PAID
> INC[,2:6]=PAID[,2:6]-PAID[,1:5]
> INC
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163   39   17    7   21
[2,] 3367 1292   37   24   10   NA
[3,] 3871 1474   -7   82   NA   NA
[4,] 4239 1678  103   NA   NA   NA
[5,] 4929 1865   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA

Comme on peut le voir, sur la troisième ligne, on a un incrément de paiementnégatif. A priori, ce n’est pas forcément gênant. En particulier, on peut faire tourner la méthode Chain Ladder,

> lambda=rep(NA,5)
> for(k in 1:5){
+ lambda[k]=(sum(PAID[1:(6-k),k+1])/
+             sum(PAID[1:(6-k),k]))}
> lambda
[1] 1.380933 1.008476 1.008515 1.001858 1.004735
> PROJECTION=PAID
> for(k in 1:5){
+ PROJECTION[((7-k):6),k+1]=
+ PROJECTION[((7-k):6),k]*lambda[k]}
> sum(PROJECTION[,6]-
+ diag(PROJECTION[,6:1]))
[1] 2469.703

et la sortie coïncide avec la fonction R,

> MackChainLadder(PAID)
MackChainLadder(Triangle = PAID)

Latest Dev.To.Date Ultimate    IBNR Mack.S.E CV(IBNR)
1  4,456       1.000    4,456     0.0    0.000      NaN
2  4,730       0.995    4,752    22.4    0.146  0.00652
3  5,420       0.993    5,456    35.8    2.405  0.06721
4  6,020       0.985    6,111    91.3   41.679  0.45629
5  6,794       0.977    6,956   161.5   71.620  0.44334
6  5,217       0.707    7,376 2,158.6   95.750  0.04436

Totals
Latest:            32,637.00
Dev:                    0.93
Ultimate:          35,106.70
IBNR:               2,469.70
Mack S.E.:            146.62
CV(IBNR):  0.059366227164502
Message d'avis :
In Mack.S.E(CL[["Models"]], FullTriangle, est.sigma = est.sigma,
'loglinear' model to estimate sigma_n doesn't appear appropriate
p-value > 5.
est.sigma will be overwritten to 'Mack'.
Mack's estimation method will be used instead.

On notera le message d’avis, laissant entendre qu’il peut y avoir un petit soucis. En fait, le soucis apparaît clairement si on souhaite faire une régression de Poisson. Car autant les valeurs non-entières ne sont pas trop gênantes, autant les valeurs négatives le sont !

> Y=as.vector(INC)
> D=rep(1:6,each=6)
> A=rep(2001:2006,6)
> base=data.frame(Y,D,A)
> reg=glm(Y~as.factor(D)+as.factor(A),
+     data=base,family=poisson(link="log"))
Erreur dans eval(expr, envir, enclos) :
les valeurs négatives sont interdites pour la famille poisson

Il faut alors trouver une solution si on a des incréments négatifs. Et aucune solution ne sera intellectuellement satisfaisante, car avec une valeur négative, il n’y a aucune légitimité pour continuer à utiliser un modèle Poissonnien. Donc on va bricoler…

  • Prendre de l’argent à droite et à gauche

La première solution consiste à dire que cet incrément n’a pas de raison d’être. On va donc aller chercher de l’argent dans la colonne avant (sur la même ligne, on va supposer qu’il s’agit d’un problème sur les cadences de paiements) ou sur la colonne après. Au lieu d’avoir

[3,] 3871 1474   -7   82   NA   NA

On peut renflouer en prenant à gauche,

[3,] 3871 1467   0    82   NA   NA

ou à droite

[3,] 3871 1474   0    75   NA   NA

En faisant ces opérations, on change seulement localement la courbe de paiements pour la troisième année de survenance. On peut se demander si le fait de prendre à gauche a un impact sur l’estimation du montant de provisions, ou sur l’incertitude associée. Pour ça, on peut utiliser les fonctions suivantes,

library(ChainLadder)
CL=function(T){
sum(
MackChainLadder(T)$FullTriangle[,ncol(T)]
-diag(T[,ncol(T):1]))}
CL.SE=function(T){
MackChainLadder(T)$Total.Mack.S.E}
CumInc=function(T){
m=T
for(i in 2:nrow(T)){m[,i]=apply(T[,1:i],1,sum)}
return(m)}

On peut alors regarder, si on prend à droite ou à gauche ce qui se passe, ou ce qui se passe si on renfloue non pas à 0 mais à une valeur strictement positive (e.g. 1),

VCL=rep(NA,8)
for(j in 1:8){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+7,j-1-7)
VCL[j]=CL(CumInc(INCd))}
plot(1:8,VCL,type="b",col="blue",xlim=c(0,9),ylim=c(2462,2474))
VCL=rep(NA,10)
for(j in 1:10){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+9,j-1-9)
VCL[j]=CL(CumInc(INCd))}
lines(0:9,VCL,type="b",pch=0,col="red")

avec la courbe bleu si on renfloue à 0, et rouge si on renfloue à 1, pour l’estimation du montant total de provisions,

VCL=rep(NA,8)
for(j in 1:8){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+7,j-1-7)
VCL[j]=CL.SE(CumInc(INCd))}
plot(1:8,VCL,type="b",col="blue",xlim=c(0,9),ylim=c(130,145))
VCL=rep(NA,10)
for(j in 1:10){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+9,j-1-9)
VCL[j]=CL.SE(CumInc(INCd))}
lines(0:9,VCL,type="b",pch=0,col="red")

avec cette fois l’impact des transferts sur l’écart-type. En abscisse, on a (à peu de choses près) le montant enlevé sur la colonne de droite. Bref, les transferts venant de la gauche ou de la droite sont une solution, mais le choix d’où vient l’argent ne sera pas neutre, ni sur l’estimation, ni sur la variance de l’estimation (et l’intervalle de confiance).

  • Jouer à faire des translations…

Une autre piste pourrait être de noter que, dans le modèle linéaire, si on translate nos données (vers le haut), on ne change pas la pente, et la constante est juste augmenté (de la taille de la translation).

> lm(dist~speed,data=cars)

Call:
lm(formula = dist ~ speed, data = cars)

Coefficients:
(Intercept)        speed
-17.579        3.932

> lm((dist+10)~speed,data=cars)

Call:
lm(formula = (dist + 10) ~ speed, data = cars)

Coefficients:
(Intercept)        speed
-7.579        3.932

Autrement dit, la prédiction faite par notre modèle n’est pas modifiée par la translation, à condition d’opérer la translation ensuite sur la prédiction. Sur le dessin ci-dessous, on fait pareil, mais avec une régression de Poisson: la courbenoire est sur les données brutes, la courbe bleu est obtenue sur les points translatés vers le haut, et la rouge si on translate la prédiction sur les points translatés,
http://freakonometrics.blog.free.fr/public/perso4/glm-translation.gif
Si la translation est trop importante, la prédiction finale (en rouge) se retrouve davantage éloignée de la vraie valeur (en noir).
Pareillement, on peut translater nos paiements vers le haut dans les triangles, de manière à avoir des paiements positifs. On peut commencer par translater tout le triangle, par exemple en translatant juste assez pour que tous les incréments soient positifs ou nuls,

> k=7
> Y=as.vector(INC)+k
> D=rep(1:6,each=6)
> A=rep(2001:2006,6)
> base=data.frame(Y,D,A)
> reg=glm(Y~as.factor(D)+as.factor(A),
+ data=base,family=poisson(link="log"))
> Yp=predict(reg,type="response",
+ newdata=base)-k
> sum(Yp[is.na(Y)==TRUE])
[1] 2508.620

On peut aussi se demander ce qui se passerait si on décalait davantage vers le haut. On l’a vu sur l’animation, plus la translation est importante, plus la prédiction se décale. Un stratégie pourrait être de faire une estimation pour plusieurs valeurs de translations – de telle sorte que tous les incréments soient positifs ou nuls – et ensuite d’extrapoler pour savoir ce qui se passerait si on ne translatait pas,

> K=7:20
> R=rep(NA,length(K))
> for(i in 1:length(K)){
+ k=K[i]
+ Y=as.vector(INC)+k
+ D=rep(1:6,each=6)
+ A=rep(2001:2006,6)
+ base=data.frame(Y,D,A)
+ reg=glm(Y~as.factor(D)+as.factor(A),
+ data=base,family=poisson(link="log"))
+ Yp=predict(reg,type="response",
+ newdata=base)-k
+ R[i]=sum(Yp[is.na(Y)==TRUE])}
> plot(K,R,xlim=c(0,20),ylim=c(2465,max(R)))
> abline(lm(R~K),col="blue")
> (yp=predict(lm(R~K),newdata=(K=0)))
1
2470.199
> points(0,yp,col="red",pch=19)

 

On a cette fois une estimation plus proche de celle que nous avions en utilisant directement, sur le triangle brut, la méthode chain ladder. Mais on peut aller un peu plus loin. Car dans le triangle, on a translaté toutes les valeurs, alors que seule une était négative. On pourrait par exemple translater seulement les valeurs de la troisième colonne,

> K=7:20
> R=rep(NA,length(K))
> for(i in 1:length(K)){
+ k=K[i]
+ Y=as.vector(INC)
+ D=rep(1:6,each=6)
+ A=rep(2001:2006,6)
+ Y[D==3]=Y[D==3]+k
+ base=data.frame(Y,D,A)
+ reg=glm(Y~as.factor(D)+as.factor(A),
+ data=base,family=poisson(link="log"))
+ Yp=predict(reg,type="response",
+ newdata=base)
+ Yp[D==3]=Yp[D==3]-k
+ R[i]=sum(Yp[is.na(Y)==TRUE])}
> predict(lm(R~K),newdata=(K=0))
1
2469.703

Tiens, on retombe exactement sur l’estimateur de la méthode chain ladder… Étonnant, non ?

Une petite précision. Dans la vraie vie, on ne devrait pas avoir d’incréments négatifs dans les triangles de paiements. Maintenant, il faut reconnaître qu’on en voit apparaître lorsque l’on bootstrappe les résidus, et que l’on génère de pseudo-triangles…
Si des praticiens ont des commentaires sur les incréments négatifs, les commentaires sont ouverts (et peu modérés), donc racontez nous comment vous faites, je suis preneur…

ACT2040: examen final (suite)

L’examen final pour le cours d’actuariat IARD aura lieu mardi prochain, et il sera basé sur des données mises en ligne sur le blog depuis une dizaine de jours. Je mets aujourd’hui en ligne les sorties qui seront fournies le jour de l’examen, et sur lesquelles des questions seront posées. J’encourage les étudiants à survoler ces sorties avant l’examen. En revanche, je ne répondrai à aucune question sur ce document ! Amenez vos calculatrices à l’examen.

Une région géographique n’est pas une variable continue

En relisant les devoirs maisons, je me suis rendu compte que certains avaient tenté de regrouper les régions (géographiques) par régions homogènes. Sauf que les régions étaient codées par un numéro (selon la codification officielle). Par exemple, dans une des bases, nous avions des assurés dans 4 zones géographiques, à savoir la région 82 (région Rhône-Alpes en rouge) la région 54 (région Poitou-Charentes en vert) la région 73 (région Midi-Pyrénées en bleu) et enfin la région 41 (région Lorraine en mauve).

> unique(baseFREQ$region) 
[1] 82 54 73 41

Une idée intéressante pour regrouper les régions pouvait être d’utiliser les arbres. Les régions étant des couleurs (on le voit bien sur la carte) et pas des variables quantitatives, il est normal de travailler sur des facteurs. D’ailleurs le code pour faire la carte est le suivant,

> library(maps) 
>  france<-map(database="france") 
>  dpt=c("Ain","Ardeche","Drome","Isere","Loire ","Rhone",  
+ "Savoie","Haute-Savoie","Charente","Charente-Maritime", 
+ "Deux-Sevres","Vienne","Ariege","Aveyron","Haute-Garonne",  
+ "Gers","Lot","Hautes-Pyrenees","Tarn","Tarn-et-Garonne", 
+ "Meurthe-et-Moselle","Meuse","Moselle","Vosges") 
>  couleur=c(rep(2,8),rep(3,4),rep(4,8),rep(6,4))  
>  match=match.map(france,dpt) 
>  color=couleur[match] 
>  map(database="france", fill=TRUE, col=color)

L’arbre sur les régions en tant que facteurs donne le découpage suivant

>  baseFREQ$fregion=as.factor(baseFREQ$region) 
>  ARBRE1=tree(nombre~fregion,data=baseFREQ,split="gini")  
>  plot(ARBRE1) 
>  text(ARBRE1)

Bon, R a la mauvaise idée de recoder les classes (mais il garde l’ordre, i.e. a correspond à la région 41, b à 54, c à 73 et d à 82). Visuellement, on retient qu’il est possible de considérer deux grandes régions, AC (i.e. 41 et 73) et BD (i.e. 54 et 72). L’intérêt des arbres sur des variables qualitatives, des facteurs, c’est que tous les regroupements sont possibles. En revanche, si on fait un arbre sur la région qui est lue en tant que nombre (quantitatif), on obtient

>  ARBRE2=tree(nombre~region,data=baseFREQ,split="gini") 
>  plot(ARBRE2) 
>  text(ARBRE2)

Il est alors impossible de regrouper dans une même classe deux régions séparées par un nombre, i.e. on ne peut regrouper 41 et 82 dans la même classe. R suggère de distinguer peut être trois régions, à savoir 82 (à droite), puis 73 (au centre) et enfin de mettre éventuellement 41 et 54 ensemble. Ce qui n’est pas la stratégie optimale quand on regroupe des facteurs.

ACT2040: examen final

L’examen final du cours ACT2040 aura lieu le 13 décembre, durant trois heures. L’examen sera en deux parties: une de tarification, et une de provisionnement. Les deux porteront sur la lecture et l’analyse critique de sorties informatiques. Pour la tarification, les bases de données utilisées sont en ligne ci-dessous

> BASEN=read.table("http://freakonometrics.free.fr/baseN.txt",header=TRUE,sep=";")
> BASEY=read.table("http://freakonometrics.free.fr/baseY.txt",header=TRUE,sep=";")
> head(BASEN)
ageconducteur agepermis sexeconducteur situationfamiliale  habitation zone
1            57        39              F             Celiba peri-urbain    8
2            54        35              H             Celiba      urbain    3
3            51        32              F             Celiba      urbain    1
4            53        35              H              Marie       rural    4
5            61        43              H              Marie      urbain    8
6            60        29              F              Marie peri-urbain    1
agevehicule proprietaire    payment  marque         poids     usage
1          12    locataire     Annuel  AUTRES     8.>3500kg PROMENADE
2          20     sans mrp Semestriel PEUGEOT 4.3100-3199kg PROMENADE
3           4     sans mrp     Annuel  RAPIDO     1.<2700kg PROMENADE
4           1     sans mrp     Annuel  AUTRES 3.3000-3099kg PROMENADE
5           1 proprietaire     Annuel    FIAT 6.3300-3399kg PROMENADE
6          10     sans mrp    Mensuel    FIAT     8.>3500kg PROMENADE
exposition nombre   voiture
1          1      0 Monospace
2          1      0   Berline
3          1      0  sans avp
4          1      0  sans avp
5          1      1 Monospace
6          1      0  sans avp

Parmi les variables, la description (sommaire) est la suivante,

  • ageconducteur: âge du conducteur principal du véhicule
  • agepermis: ancienneté du permis de conduire du conducteur principal du véhicule
  • sexeconducteur: sexe du conducteur principal (H ou F)
  • situationfamiliale: situation familiale du conducteur principal (“Celiba”, “Marie” ou “Veuf/Div”)
  • habitation: zone d’habitation du conducteur principal (“peri-urbain”, “rural” ou “urbain” )
  • zone: zone d’habitation (allant de 1 à 8)
  • agevehicule: age du véhicule
  • proprietaire: si le conducteur principal possède un contrat Habitation, son statut (“locataire” ou “proprietaire”)  Sinon “sans mrp”
  • payment:type de fractionnement de la prime d’assurance automobile (“Annuel”, “Mensuel” ou “Semestriel”)
  • marque: marque du véhicule
> levels(BASEN[,10])
[1] "ADRIA"       "AUTOSTAR"    "AUTRES"      "BURSTNER MOBIL"
[5] "CHALLENGER"  "CHAUSSON"    "CITROEN"     "FIAT"
[9] "FORD"        "HYMERMOBIL"  "MERCEDES"    "PEUGEOT"
[13] "PILOTE"     "RAPIDO"      "RENAULT"     "VOLKSWAGEN"
  • poids: classe de poids du véhicule
> levels(BASEN[,11])
[1] "1.<2700kg"    "2.2700-2999kg""3.3000-3099kg""4.3100-3199kg"
[5] "5.3200-3299kg""6.3300-3399kg""7.3400-3499kg""8.>3500kg"
  • usage: utilisation du véhicule principal (“PROMENADE” ou “TOUS_DEPLACEMENTS”)
  • exposition: exposition, en années
  • nombre: nombre d’accident responsabilité civile du conducteur principal, pendant l’année passée
  • cout: cout du sinistre
  • voiture: type de véhicule
> levels(BASEN[,15])
[1] "Berline"            "Break"              "Buggy"
[4] "Cabriolet"          "Combispace"         "Coup\xe9"
[7] "Coup\xe9 Cabriolet" "Jeep"               "Minibus"
[10] "Minispace"          "Monospace"          "sans avp"

Concernant la partie sur le provisionnement, des sorties associées au triangle suivant seront proposées,

> TRIANGLE=read.table("http://freakonometrics.blog.free.fr/public/code/triangleACT2040.txt")
> TRIANGLE
X1      X2      X3      X4      X5      X6      X7      X8      X9     X10
1  357848 1124788 1735330 2218270 2745596 3319994 3466336 3606286 3833515 3901463
2  352118 1236139 2170033 3353322 3799067 4120063 4647867 4914039 5339085      NA
3  290507 1292306 2218525 3235179 3985995 4132918 4628910 4909315      NA      NA
4  310608 1418858 2195047 3757447 4029929 4381982 4588268      NA      NA      NA
5  443160 1136350 2128333 2897821 3402672 3873311      NA      NA      NA      NA
6  396132 1333217 2180715 2985752 3691712      NA      NA      NA      NA      NA
7  440832 1288463 2419861 3483130      NA      NA      NA      NA      NA      NA
8  359480 1421128 2864498      NA      NA      NA      NA      NA      NA      NA
9  376686 1363294      NA      NA      NA      NA      NA      NA      NA      NA
10 344014      NA      NA      NA      NA      NA      NA      NA      NA      NA

Confidence interval for predictions with GLMs

Consider a (simple) Poisson regression https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss01.gif. Given a sample https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss02.gif where https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss03.gif, the goal is to derive a 95% confidence interval for https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss04.gif given https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss05.gif, where https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss04.gif is the prediction. Hence, we want to derive a confidence interval for the prediction, not the potential observation, i.e. the dot on the graph below

> r=glm(dist~speed,data=cars,family=poisson)
> P=predict(r,type="response",
+ newdata=data.frame(speed=seq(-1,35,by=.2)))
> plot(cars,xlim=c(0,31),ylim=c(0,170))
> abline(v=30,lty=2)
> lines(seq(-1,35,by=.2),P,lwd=2,col="red")
> P0=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> points(30,P1$fit,pch=4,lwd=3)

i.e.

Let https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif denote the maximum likelihood estimator of https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss07.gif. Then
https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss40.gif
where https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss101.gif is Fisher information of https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif (from standard maximum likelihood theory). Recall that
https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss13.gif
where computation of those values is based on the following calculations
http://freakonometrics.blog.fre<br /><br /> e.fr/public/latex/poiss21.gif
In the case of the log-Poisson regression
https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss36.gif
Let us get back to our initial problem.

  • confidence interval for the linear combination

A first idea to get a confidence interval for https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss49.gif is to get a confidence interval for https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss100.gif (by taking exponential values of bounds, since the exponential is a monotone function). Asymptotically, we know that
https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss40.gif

thus, an approximation for the variance matrix of https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif will be based on https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss45.gif, obtained by plugging estimators of the parameters.
Then, since https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif as an asymptotic multivariate distribution, any linear combination of the parameters will also be normal, i.e.
https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss47.gif has a normal distribution, centered on https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss49.gif, with variance https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss102.gif where https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/Poiss110.gif is the variance of https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif. All those quantities can be easily computed. First, we can get the variance of the estimators

> i1=sum(predict(reg,type="response"))
> i2=sum(cars$speed*predict(reg,type="response"))
> i3=sum(cars$speed^2*predict(reg,type="response"))
> I=matrix(c(i1,i2,i2,i3),2,2)
> V=solve(I)

Hence, if we compare with the output of the regression,

> summary(reg)$cov.unscaled
(Intercept)         speed
(Intercept)  0.0066870446 -3.474479e-04
speed       -0.0003474479  1.940302e-05
> V
[,1]          [,2]
[1,]  0.0066871228 -3.474515e-04
[2,] -0.0003474515  1.940318e-05

Based on those values, it is easy to derive the standard deviation for the linear combination,

> x=30
> P2=predict(r,type="link",se.fit=TRUE,
+ newdata=data.frame(speed=x))
> P2
$fit
1
5.046034

$se.fit
[1] 0.05747075

$residual.scale
[1] 1

> sqrt(V[1,1]+2*x*V[2,1]+x^2*V[2,2])
[1] 0.05747084
> sqrt(t(c(1,x))%*%V%*%c(1,x))
[,1]
[1,] 0.05747084

And once we have the standard deviation, and normality (at least asymptotically), confidence intervals are derived, and then, taking the exponential of the bounds, we get confidence interval

> segments(30,exp(P2$fit-1.96*P2$se.fit),
+ 30,exp(P2$fit+1.96*P2$se.fit),col="blue",lwd=3)

Based on that technique, confidence intervals are no longer centered on the prediction. But who cares ?

  • delta method

Actually, those who like to use “more or less” expressions for confidence intervals will not like non centered intervals. So, an alternative is to use the delta method. Instead of writing (again) something on the theory, we can use a package which computes that method,

> estmean=t(c(1,x))%*%coef(reg)
> var=t(c(1,x))%*%summary(reg)$cov.unscaled%*%c(1,x)
> library(msm)
> deltamethod (~ exp(x1), estmean, var)
[1] 8.931232
> P1=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> P1
$fit
1
155.4048

$se.fit
1
8.931232

$residual.scale
[1] 1

The delta method gives us (asymptotic) normality, so once we have a standard deviation, we get the confidence interval.

> segments(30,P1$fit-1.96*P1$se.fit,30,
+ P1$fit+1.96*P1$se.fit,col="blue",lwd=3)

Note that those quantities – obtained with two different approaches – are rather close here

> exp(P2$fit-1.96*P2$se.fit)
1
138.8495
> P1$fit-1.96*P1$se.fit
1
137.8996
> exp(P2$fit+1.96*P2$se.fit)
1
173.9341
> P1$fit+1.96*P1$se.fit
1
172.9101
  • bootstrap techniques

And a third method (but far from what I expect to teach on that course) is to use bootstrap techniques to about those results based on asymptotic normality (we have only 50 observations). The idea is to sample from out dataset, and to run a log-Poisson regression on those new samples, and to repeat a lot of time,

Régression sur des variables catégorielles

Petit complément par rapport au cours de mardi. On avait évoqué tout d’abord la lecture des sorties lorsque l’on régresse sur des variables catégorielles (des facteurs). Commençons par supprimer la constante de la régression

> reg0=glm(nbre~0+zone,offset=log(exposition),data=base, 
+ family=poisson(link="log"))
> summary(reg0)

Call:
glm(formula = nbre ~ 0 + zone, family = poisson(link = "log"), 
    data = base, offset = log(exposition))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5717  -0.3968  -0.2996  -0.1547  12.6722  

Coefficients:
      Estimate Std. Error z value Pr(>|z|)    
zoneB -2.54187    0.06287  -40.43   <2e-16 ***
zoneA -2.54912    0.05285  -48.23   <2e-16 ***
zoneC -2.38525    0.03753  -63.56   <2e-16 ***
zoneD -2.13454    0.03878  -55.05   <2e-16 ***
zoneE -2.00204    0.03965  -50.49   <2e-16 ***
zoneF -2.06932    0.11547  -17.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 50966  on 50000  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(reg0,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
        1         2         3         4         5 
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On voit que toutes les modalités sont présentes, et toutes sont significatives. Si on régresse sur la constante, il faudra supprimer une modalité pour rendre le modèle identifiable. On peut forcer pour que la modalité de référence soit la seconde,

> base$zone=relevel(base$zone,"B")
> regB=glm(nbre~zone,offset=log(exposition),data=base,
+ family=poisson(link="log"))
> summary(regB)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54187    0.06287 -40.431  < 2e-16 ***
zoneA       -0.00725    0.08213  -0.088 0.929661
zoneC        0.15662    0.07322   2.139 0.032432 *
zoneD        0.40733    0.07387   5.514 3.50e-08 ***
zoneE        0.53983    0.07433   7.263 3.80e-13 ***
zoneF        0.47255    0.13148   3.594 0.000325 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(regB,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
1         2         3         4         5
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On notera que les prédictions ne changent pas. On peut aussi choisir la première comme modalité de référence,

> base$zone=relevel(base$zone,"A")
> reg=glm(nbre~zone,offset=log(exposition),
> data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54912    0.05285 -48.232  < 2e-16 ***
zoneB        0.00725    0.08213   0.088 0.929661
zoneC        0.16387    0.06482   2.528 0.011471 *
zoneD        0.41458    0.06555   6.324 2.54e-10 ***
zoneE        0.54708    0.06607   8.280  < 2e-16 ***
zoneF        0.47980    0.12699   3.778 0.000158 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Le fait que la seconde modalité ne soit pas significative se lit par rapport à la modalité de référence (en l’occurrence la première): non significatif signifie alors non significativement différente. Autrement dit, on peut regrouper les modalités en une seule.

> base$zonesimple=base$zone
> base$zonesimple[base$zone%in%c("A","B")]="A"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54612    0.04046 -62.937  < 2e-16 ***
zonesimpleC  0.16087    0.05518   2.915  0.00355 **
zonesimpleD  0.41158    0.05604   7.345 2.06e-13 ***
zonesimpleE  0.54408    0.05665   9.605  < 2e-16 ***
zonesimpleF  0.47681    0.12235   3.897 9.74e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

On note qu’avec ce regroupement, les autres modalités sont sensiblement différentes. On peut aussi retenir la troisième comme modalité de référence

> base$zonesimple=relevel(base$zonesimple,"C")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.38525    0.03753 -63.557  < 2e-16 ***
zonesimpleA -0.16087    0.05518  -2.915  0.00355 **
zonesimpleD  0.25071    0.05396   4.646 3.39e-06 ***
zonesimpleE  0.38321    0.05460   7.019 2.24e-12 ***
zonesimpleF  0.31593    0.12142   2.602  0.00927 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Comme toutes les modalités semblent significatives, on peut tenter de prendre comme modalité de référence une des dernières (dont les estimations des coefficients donnent des résultats très proches)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.06932    0.11547 -17.921  < 2e-16 ***
zonesimpleC -0.31593    0.12142  -2.602  0.00927 **
zonesimpleA -0.47681    0.12235  -3.897 9.74e-05 ***
zonesimpleD -0.06522    0.12181  -0.535  0.59232
zonesimpleE  0.06727    0.12209   0.551  0.58161
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Au vue de cette dernière sortie, on peut tenter de fusionner toutes les dernières classes ensembles

> base$zonesimple[base$zone%in%c("D","E","F")]="F"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5660  -0.3959  -0.3004  -0.1547  12.5929

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.07182    0.02696 -76.853  < 2e-16 ***
zonesimpleC -0.31344    0.04621  -6.783 1.18e-11 ***
zonesimpleA -0.47431    0.04861  -9.757  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15698  on 49997  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Bon, formellement, regrouper deux modalités (i.e. décréter que deux variables sont non significatives simultanément) demande un peu plus qu’un test de Student, ou que deux tests de Student…. Si on remonte un peu en arrière, on aurait pu faire un test multiple avant de fusionner les trois modalités (un test de type Fisher, ou une ANOVA)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> library(car)
> linearHypothesis(reg,c("zonesimpleD=0","zonesimpleE=0"))
Linear hypothesis test

Hypothesis:
zonesimpleD = 0
zonesimpleE = 0

Model 1: restricted model
Model 2: nbre ~ zonesimple

Res.Df Df  Chisq Pr(>Chisq)
1  49997
2  49995  2 5.7073    0.05763 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Manifestement, on peut accepter l’hypothèse que ces trois catégories n’en font qu’une. La zone géographique peut alors se découper en trois grandes zones, et pas en six. On notera que cela correspond à ce que propose un arbre de régression

> library(tree)
> arbre=tree(nbre~zone+offset(log(exposition)),
+ data=base,split="gini")
> plot(arbre)
> text(arbre)

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 http://freakonometrics.blog.free.fr/public/perso4/Yi-ltx.gif sera de la forme

http://freakonometrics.blog.free.fr/public/perso4/loi-exponentielle.gif

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

ACT2040, exercices de tarification

Les exercices commenceront ce vendredi, et les énoncés (ainsi que pour quelques uns des pistes de correction) sont en ligne ici. Les exercices (et les corrections) ont été tapés pour ce cours par Jean-Philippe Boucher. La premiere série d’exercices sont des exercices de remise en route sur les calculs d’espérances, et d’espérance conditionnelle. On reviendra dessus mardi prochain.