Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Actuariat IARD - ACT2040
Partie 3 - régression Poissonienne
et biais minimal (Y ∈ N)
Arthur Charpentiercharpentier.arthur@uqam.ca
http ://freakonometrics.hypotheses.org/

Automne 20131

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Modélisation d’une variable de comptage
Références : Frees (2010), chapitre 12 (p 343-361) Greene (2012), section 18.3
(p 802-828) de Jong Heller (2008), chapitre 6, sur la régression de Poisson. Sur
les méthodes de biais minimal, de Jong Heller (2008), section 1.3. Sinon
Cameron, A.C. & Trivedi, P.K. (1998). Regression analysis of count data,
Cambridge University Press.
Denuit, M., Marechal, X., Pitrebois, S. & Walhin, J.-F. (2007) Actuarial
Modelling of Claim Counts : Risk Classification, Credibility and Bonus-Malus
Systems. Wiley.
Hilbe, J. M. (2007). Negative Binomial Regression, Cambridge University Press.
Remarque : la régression de Poisson est un cas particulier des modèles GLM,
avec une loi de Poisson et une fonction de lien logarithmique.
regPoisson = function(formula, lien="log", data=NULL) {glm(formula,family=poisson(link=log),data) }2

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Base pour les données de comptage> sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",+ header=TRUE,sep=";")> sinistres=sinistre[sinistre$garantie=="1RC",]> sinistres=sinistre[sinistres$cout>0,]> contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt",+ header=TRUE,sep=";")> T=table(sinistres$nocontrat)> T1=as.numeric(names(T))> T2=as.numeric(T)> nombre1 = data.frame(nocontrat=T1,nbre=T2)> I = contrat$nocontrat%in%T1> T1= contrat$nocontrat[I==FALSE]> nombre2 = data.frame(nocontrat=T1,nbre=0)> nombre=rbind(nombre1,nombre2)> baseFREQ = merge(contrat,nombre)3Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Base pour les données de comptage
La variable d’intérêt Y est ici nbre ou baseFREQ$nbre> head(baseFREQ,4)
nocontrat exposition zone puissance agevehicule ageconducteur bonus marque
1
27
0.87
C
7
0
56
50
12
2
115
0.72
D
5
0
45
50
12
3
121
0.05
C
6
0
37
55
12
4
142
0.90
C
10
10
42
50
12
carburant densite region nbre
1
D
93
13
0
2
E
54
13
0
3
D
11
13
0
4
D
93
13
0

Remarque : la base ici est une sous-base (particulière) de celle utilisée dans
http ://cran.r-project.org/doc/contrib/Charpentier_Dutang_actuariat_avec_R.pdf

Les conclusions risques d’être différentes aussi...
4Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La loi de Poisson
La loi de Poisson est connue comme la loi des petits nombres,−λ λP(N = k) = ekk!

, ∀k ∈ N.

Pour rappel (cf premier cours), E(N ) = Var(N ).

5Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La loi Poisson mélange
En présence de sur-dispersion E(N ) < Var(N ), on peut penser à une loi Poisson
mélange, i.e. il existe Θ, variable aléatoire positive, avec E(Θ) = 1, telle quek−λθ [λθ]P(N = k|Θ = θ) = ek!

Frequence pour la classe (1)Frequence pour la classe (2)Frequence du melange

0.10, ∀k ∈ N.

Frequence pour la classe (1)Frequence pour la classe (2)Frequence du melange
0.08

0.08

0.06
0.06

0.04
0.04

0.02
0.02

0.00

0.00
0

2
1

4
3

6
5

8
7

10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50
9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

0

2
1

4
3

6
5

8
7

9

10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50
11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

6

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La loi Binomiale Négative
La loi binomiale négative apparît dans le modèle Poisson mélange, lorsqueΘ ∼ G(α, α). Dans ce cas,ααexp (−αx)
α−1π(θ) = xΓ(α)E(Θ) = 1 et Var(Θ) =1.
αDans ce cas, la loi (non conditionnelle) de N est
Z ∞P(N = k) =P(N = k|Θ = θ)π(θ)dθ,
0

P(N = k) =Γ(k + α)Γ(k + 1)Γ(α)



11 + λ/αα 
1−

11 + λ/αk, ∀k ∈ N

7

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Le processus de Poisson
Pour rappel, (Nt )t≥0 est un processus de Poisson homogène (de paramètre λ) s’ilest à accroissements indépendants, et le nombre de sauts observés pendant lapériode [t, t + h] suit une loi P(λ · h).

Ns+1 − Ns ∼ P(λ) est indépendant de Nt+h − Nt ∼ P(λ · h).

8Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Exposition et durée d’observation
Soit Ni la frénquence annulisée de sinistre pour l’assuré i, et supposons
Ni ∼ P(λ).
Si l’assuré i a été observé pendant une période Ei , le nombre de sinistre observé
est Yi ∼ P(λ · Ei ).

L(λ, Y , E) =nY
e−λEi [λEi ]Yi

i=1log L(λ, Y , E) = −λ

n
X

Ei +

i=1

n
X

Yi !
Yi log[λEi ] − log

i=1

n
Y

!
Yi !

i=1

qui donne la condition du premier ordre
n
n
X
1X
∂
Ei +
Yi
log L(λ, Y , E) = −
∂λ
λ
i=1
i=1

9

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

qui s’annule pour
Pn
n
X
Y
Yi
Ei
i
i=1
b
λ = Pn=ωiavec ωi = Pn
Ei
i=1 Eii=1 Eii=1> N <- baseFREQ$nbre> E <- baseFREQ$exposition> (lambda=sum(N)/sum(E))[1] 0.07279295> weighted.mean(N/E,E)[1] 0.07279295> dpois(0:3,lambda)*100[1] 92.979 6.768 0.246 0.006

Remarque : pour Ei on parlera d’exposition ou d’années police.

10Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Fréquence annuelle, une variable tarifaire> X1 <- baseFREQ$carburant> sum(N[X1=="D"])[1] 998> sum(N[X1=="E"])[1] 926> tapply(N, X1, sum)DE
998 926

On a ici le nombre d’accidents par type de carburant.
> tapply(N, X1, sum)/tapply(E, X1, sum)DE
0.07971533 0.06656323

11

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Fréquence annuelle, une variable tarifaire
Notons que l’on peut comparer ces moyenne (pondérées)> library(weights)> wtd.t.test(x=(N/E)[X1=="D"], y=(N/E)[X1=="E"],+ weight=E[X1=="D"], weighty=E[X1=="E"],samedata=FALSE)$test[1] "Two Sample Weighted T-Test (Welch)"$coefficientst.valuedf
p.value
2.835459e+00 4.959206e+04 4.577838e-03$additional
Difference
Mean.xMean.yStd. Err
0.013152097 0.079715332 0.066563235 0.004638437

Avec une p-value de l’ordre de 0.004 on accepte l’hypothèse que les deuxmoyennes sont distinctes. On retrouve un résultat semblable avec une régression
linéraire
12Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Fréquence annuelle, une variable tarifaire> summary(lm(N/E~X1,weights=E))
Coefficients:
Estimate Std. Error t value Pr(>|t|)(Intercept) 0.079715
0.003386 23.540 < 2e-16 ***X1E-0.013152
0.004668 -2.818 0.00484 **--Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 0.3789 on 49998 degrees of freedom
Multiple R-squared: 0.0001588,
Adjusted R-squared: 0.0001388
F-statistic: 7.939 on 1 and 49998 DF, p-value: 0.00484

13

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Fréquence annuelle et tableau de contingence
Supposons que l’on prenne en compte ici deux classes de risques.
> N <- baseFREQ$nbre> E <- baseFREQ$exposition> X1 <- baseFREQ$carburant> X2 <- cut(baseFREQ$agevehicule,c(0,3,10,101),right=FALSE)> names1 <- levels(X1)> names2 <- levels(X2)> (POPULATION=table(X1,X2))
X2
X1 [0,3) [3,10) [10,101)
D 7492 10203
6718
E 6275
9417
9895

14

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Fréquence annuelle et tableau de contingence
Mais la population (i.e. le nombre d’assurés) ne correspond pas forcément au
nombre d’année-police.
> EXPOSITION=POPULATION> for(k in 1:nrow(EXPOSITION)){+ EXPOSITION[k,]=tapply(E[X1==names1[k]],+ X2[X1==names1[k]],sum)}> EXPOSITION
X2
X1
[0,3)
[3,10) [10,101)
D 3078.938 5653.109 3787.503
E 2735.014 5398.950 5777.619[0,3)[3,10)[10,101)D

3078.945653.113787.50E

2735.015398.95

5777.6215

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Fréquence annuelle et tableau de contingence
On peut ensuite compter les sinistres> SINISTRE=POPULATION> for(k in 1:nrow(SINISTRE)){+ SINISTRE[k,]=tapply(N[X1==names1[k]],+ X2[X1==names1[k]],sum)}> SINISTRE
X2
X1 [0,3) [3,10) [10,101)
D
242
479
277
E
163
368
395

[0,3)

[3,10)[10,101)

D

242

479

277

E

163

368

395

16

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Fréquence annuelle et tableau de contingenceet interpréter les valeurs précédantes en fréquence annuelle de sinistre,> (FREQUENCE=SINISTRE/EXPOSITION)
X2
X1
[0,3)
[3,10)[10,101)
D 0.07859854 0.08473214 0.07313526
E 0.05959749 0.06816140 0.06836727[0,3)[3,10)[10,101)D

0.07860.08470.0731E

0.0596

0.06820.068417

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Fréquence annuelle et tableau de contingence
Notons Yi,j le nombre de sinistres observé (non annualisé) lorsque la première
variable (i.e. X1 ) prend la valeur i (ici i prend deux valeurs, diesel ou essenceorginaire) et la seconde variable (i.e. C) prend la valeur j (ici j prend trois
valeurs, voiture neuve (0-3), récente (3-10) ou vieille (plus de 10 ans)). On notera
Ni,j la fréquence annualisée.
La matrice Y = [Yi,j ] est ici FREQUENCE. On suppose qu’il est possible demodéliser Y à l’aide d’un modèle multiplicatif à deux facteurs, associés à chaquedes des variables. On suppose que
Ni,j = Li · Cj .
On notera Ei,j l’exposition, i.e. EXPOSITION. L’estimation de L = (Li ) et deC = (Cj ) se fait généralement de trois manières : par moindres carrés, parminimisation d’une distance (e.g. du chi-deux) ou par un principe debalancement (ou méthode des marges).
18

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Méthode des marges, Bailey (1963)
Dans la méthode des marges (selon la terminologie de Bailey (1963),
formellement, on veut
X
X
X
Yi,j =
Ei,j Ni,j =Ei,j Li · Cj ,jjj

en somment sur la ligne i, pour tout i, ou sur la colonne j,
X
X
X
Yi,j =
Ei,j Ni,j =Ei,j Li · Cj ,iii

La première équation donneP
Li = Pj Yi,jj Ei,j Cj

et la secondePYi,jCj = P i.
i Ei,j Li
19

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Méthode des marges, Bailey (1963)
On résoud alors ce petit systeme de maniere itérative (car il n’y a pas de solution
analytique simple).
> (m=sum(SINISTRE)/sum(EXPOSITION))[1] 0.07279295> L<-matrix(NA,100,2);C<-matrix(NA,100,3)> L[1,]<-rep(m,2);colnames(L)=names1> C[1,]<-rep(m,3);colnames(C)=names2> for(j in 2:100){+ L[j,1]<-sum(SINISTRE[1,])/sum(EXPOSITION[1,]*C[j-1,])+ L[j,2]<-sum(SINISTRE[2,])/sum(EXPOSITION[2,]*C[j-1,])+ C[j,1]<-sum(SINISTRE[,1])/sum(EXPOSITION[,1]*L[j,])+ C[j,2]<-sum(SINISTRE[,2])/sum(EXPOSITION[,2]*L[j,])+ C[j,3]<-sum(SINISTRE[,3])/sum(EXPOSITION[,3]*L[j,])+ }20Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Méthode des marges, Bailey (1963)
On obtient ici les séquences suivantes> L[1:5,]
DE[1,] 0.07279295 0.07279295[2,] 1.09509693 0.91441875[3,] 1.09368399 0.91569373[4,] 1.09366449 0.91571133[5,] 1.09366422 0.91571157> C[1:5,][0,3)[3,10)[10,101)[1,] 0.07279295 0.07279295 0.07279295[2,] 0.06896336 0.07611702 0.07125554[3,] 0.06897350 0.07612457 0.07124032[4,] 0.06897364 0.07612468 0.07124011[5,] 0.06897364 0.07612468 0.07124011

ce qui donne les prédictions pour la fréquence21

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

> PREDICTION1=SINISTRE> PREDICTION1[1,]<-L[100,1]*C[100,]> PREDICTION1[2,]<-L[100,2]*C[100,]> PREDICTION1
X2
X1
[0,3)
[3,10)[10,101)
D 0.07543400 0.08325484 0.07791276
E 0.06315996 0.06970825 0.06523539[0,3)[3,10)[10,101)D

0.07540.08330.0779E

0.06320.06970.0652On notera que les marges sont identiques, par exemple pour la première ligne> sum(PREDICTION1[1,]*EXPOSITION[1,])[1] 998> sum(SINISTRE[1,])[1] 99822

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Remarque : Cette technique est équivalente à utiliser une régression log-Poisson
sur les deux variables qualitatives,> donnees <- data.frame(N,E,X1,X2)> regpoislog <- glm(N~X1+X2,offset=log(E),data=donnees,+ family=poisson(link="log"))> newdonnees <- data.frame(X1=factor(rep(names1,3)),E=rep(1,6),+ X2=factor(rep(names2,each=2)))> matrix(predict(regpoislog,newdata=newdonnees,+ type="response"),2,3)
[,1]
[,2]
[,3]
[1,] 0.07543401 0.08325484 0.07791276[2,] 0.06315996 0.06970825 0.06523539[0,3)[3,10)[10,101)D

0.07540.08330.0779E

0.06320.06970.065223

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Méthode des moindres carrésParmi les méthodes proches de celles évoquées auparavant sur la méthode desmarges, il est aussi possible d’utiliser une méthode par moindres carrés(pondérée). On va chercher à minimiser la somme des carrés des erreurs, i.e.
XD=
Ei,j (Nij − Li · Cj )2
i,j

P
∂D
La condition du premier ordre donne ici ∂L=−2
j Cj Ei,j (Ni.j − Li · Cj ) = 0
i
soit
P
P
C
E
N
j Cj Yi.jj j i,j i.j=P
Li = P
2
2
E
C
E
C
i,j
i,j
j
j
j
j
L’autre condition du premier ordre donnePP
i Li Ei,j Ni.j
i Li Yi.j
P
P
Cj =
=
2
2
E
L
E
L
i,j
i,j
i
i
i
i
On résoud alors ce petit systeme de maniere itérative (car il n’y a pas de solution
analytique simple).
24

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

> L=matrix(NA,100,2);C=matrix(NA,100,3)> L[1,]=c(1,1);colnames(L)=rownames(EXPOSITION)> C[1,]=c(1,1,1);colnames(C)=colnames(EXPOSITION)> for(j in 2:100){+ L[j,1]=sum(SINISTRE[1,]*C[j-1,])/sum(EXPOSITION[1,]*C[j-1,]^2)+ L[j,2]=sum(SINISTRE[2,]*C[j-1,])/sum(EXPOSITION[2,]*C[j-1,]^2)+ C[j,1]=sum(SINISTRE[,1]*L[j,])/sum(EXPOSITION[,1]*L[j,]^2)+ C[j,2]=sum(SINISTRE[,2]*L[j,])/sum(EXPOSITION[,2]*L[j,]^2)+ C[j,3]=sum(SINISTRE[,3]*L[j,])/sum(EXPOSITION[,3]*L[j,]^2)+ }> PREDICTION2[1,]=L[100,1]*C[100,]> PREDICTION2[2,]=L[100,2]*C[100,]> PREDICTION2
X2
X1
[0,3)
[3,10)[10,101)
D 0.07575130 0.08339880 0.07753709
E 0.06342566 0.06982882 0.06492089

25

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Méthode des moindres carrés
La prédiction pour chacune des classes est alors ici[0,3)[3,10)[10,101)D

0.07580.08340.0775E

0.06340.06980.064926

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Méthode du χ2
Parmi les méthodes proches de celles évoquées dans la section ?? sur la méthodedes marges, il est aussi possible d’utiliser une méthode basée sur la distance duchi-deux. On va chercher à minimiser
Q=

X Ei,j (Ni,j − Li · Cj )2
Li · Cji,jLà encore on utilise les conditions du premier ordre, et on obtient2 P  Ei,j Yi,j

 j
Li =  P

Cj

j Ei,j Cj

 12



et une expression du même genre pour Cj .
> L=matrix(NA,100,2);C=matrix(NA,100,3)> L[1,]=c(1,1);colnames(L)=rownames(EXPOSITION)

27

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

> C[1,]=c(1,1,1);colnames(C)=colnames(EXPOSITION)> for(j in 2:100){+ L[j,1]=sqrt(sum(EXPOSITION[1,]*FREQUENCE[1,]^2/C[j-1,])/sum(EXPOSITION[1,]*C[j-1,]))+ L[j,2]=sqrt(sum(EXPOSITION[2,]*FREQUENCE[2,]^2/C[j-1,])/sum(POPULATION[2,]*C[j-1,]))+ C[j,1]=sqrt(sum(EXPOSITION[,1]*FREQUENCE[,1]^2/L[j,])/sum(EXPOSITION[,1]*L[j,]))+ C[j,2]=sqrt(sum(EXPOSITION[,2]*FREQUENCE[,2]^2/L[j,])/sum(EXPOSITION[,2]*L[j,]))+ C[j,3]=sqrt(sum(EXPOSITION[,3]*FREQUENCE[,3]^2/L[j,])/sum(EXPOSITION[,3]*L[j,]))+ }> PREDICTION3=SINISTRE> PREDICTION3[1,]=L[100,1]*C[100,]> PREDICTION3[2,]=L[100,2]*C[100,]> PREDICTION3
X2
X1
[0,3)
[3,10)[10,101)
D 0.08556959 0.09529172 0.09358149
E 0.05254307 0.05851284 0.05746269

28

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

[0,3)

[3,10)[10,101)D

0.0856

0.0953

0.0936E

0.05250.05850.0575

Remarque : on le verra par la suite, cette méthode donne les mêmes résultats
qu’une régression Gaussienne, avec un lien logarithmiqueUne régression - standard - Gaussienne
Yi= β0 + β1 X1,i + β2 X2,i + εi
Ei(avec des poids liés à l’exposition) donnerait ici>>+

reggaus <- lm(N/E~X1+X2,weight=E,data=donnees)matrix(predict(reggaus,newdata=newdonnees,
type="response"),2,3)
[,1]
[,2]
[,3]

29

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

[1,] 0.07575523 0.08296676 0.07808159[2,] 0.06279835 0.07000988 0.06512471

1

2

3

1

0.07580.08300.0781

2

0.06280.07000.0651

On va ici chercher un modèle Gaussien avec un lien log,
Yi = Ei · exp [β0 + β1 X1,i + β2 X2,i ] + εi(malheureusement ici, l’algorithme ne converge pas...)> reggauslog <- glm(N~X1+X2,offset=log(E),data=donnees,+ family=gaussian(link="log"),start=coefficients(reggaus))
Erreur dans eval(expr, envir, enclos) :impossible de trouver des valeurs initiales correctes : priÃ¨re d’en fournir

30Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La régression de Poisson
L’idée est la même que pour la régression logistique : on cherche un modèle
linéaire pour la moyenne. En l’occurence,
Yi ∼ P(λi ) avec λi = exp[X 0i β].
Dans ce modèle, E(Yi |X i ) = Var(Yi |X i ) = λi = exp[X 0i β].
Remarque : on posera parfois θi = ηi = X 0i βLa log-vraisemblance est icinX
log L(β; Y ) =[Yi log(λi ) − λi − log(Yi !)]i=1

ou encorelog L(β; Y ) =nX

Yi · [X 0i β] − exp[X 0i β] − log(Yi !)

i=1

31

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Le gradient est icin∂ log L(β; Y ) X∇ log L(β; Y ) ==(Yi − exp[X 0i β])X 0i∂βi=1

alors que la matrice Hessienne s’écrit
nX∂ 2 log L(β; Y )H(β) ==−(Yi − exp[X 0i β])X i X 0i
0
∂β∂β
i=1

La recherche du maximum de log L(β; Y ) est obtenu (numériquement) parl’algorithme de Newton-Raphson,
1. partir d’une valeur initiale β 0
2. poser β k = β k−1 − H(β k−1 )−1 ∇ log L(β k−1 )où ∇ log L(β) est le gradient, et H(β) la matrice Hessienne (on parle parfois de
Score de Fisher).
32

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Par exemple, si on régresse sur l’âge du véhicule
Y <- baseFREQ$nbre
X1 <- baseFREQ$agevehicule
X=cbind(rep(1,length(X1)),X1)

on part d’une valeur initiale (e.g. une estimation classique de modèle linéaire)
beta=lm(Y~0+X)$coefficients

On fait ensuite une boucle (avec 50,000 lignes, l’algorithme du cours 2. nefonctionne pas)
> for(s in 1:20){+ gradient=t(X)%*%(Y-exp(X%*%beta))+ hessienne=matrix(0,ncol(X),ncol(X))+ for(i in 1:nrow(X)){+ hessienne=hessienne + as.numeric(exp(X[i,]%*%beta))* (X[i,]%*%t(X[i,]))}+ beta=beta+solve(hessienne)%*%gradient}33

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Ici, on obtient les valeurs suivantes
[,1]
[,2]
[1,] 0.034651 0.0005412[2,] -0.9318582 0.001040094[3,] -1.843778 0.002294[4,] -2.624283 0.005060941[5,] -3.144156 0.009455123[6,] -3.334736 0.01274805[7,] -3.354665 0.01331968[8,] -3.354849 0.01332686[19,] -3.354849 0.01332686

i.e. on converge (et vite là encore).

34

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La régression de PoissonPb→
On peut montrer que β
β et

√

L

b − β) → N (0, I(β)−1 ).
n(βNumériquement, la encore, on peut approcher I(β)−1 qui est la varianceb(asymptotique) de notre estimateur. Or I(β) = H(β), donc les écart-types de βsont
> sqrt(diag(solve(hessienne)))[1] 0.036478386 0.003782998

35

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

On retrouve toutes ces valeurs en utilisant> summary(regPoisson(N~X1))
Coefficients:
Estimate Std. Error z value Pr(>|z|)(Intercept) -3.354849
0.036478 -91.968 < 2e-16 ***X1
0.013327
0.003783
3.523 0.000427 ***--Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1(Dispersion parameter for poisson family taken to be 1)
Null deviance: 12931
Residual deviance: 12919
AIC: 16594

on 49999
on 49998

degrees of freedom
degrees of freedom

Number of Fisher Scoring iterations: 636

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Prise en compte de l’exposition (offset)
oups, on a oublié de prendre l’exposition dans notre modèle. On a ajusté unmodèle de la forme
Yi ∼ P(λi ) avec λi = exp[β0 + β1 X1,i ]
mais on voudrait
Yi ∼ P(λi · Ei ) avec λi = exp[β0 + β1 X1,i ]ou encore
Yi ∼ P(λ̃i ) avec λ̃i = Ei · exp[β0 + β1 X1,i ] = exp[β0 + β1 X1,i + log(Ei )]Aussi, l’exposition intervient comme une variable de la régression, mais en
prenant le logarithme de l’exposition, et en forçant le paramètre à être unitaire,i.e.
Yi ∼ P(λ̃i ) avec λ̃i = Ei · exp[β0 + β1 X1,i ] = exp[β0 + β1 X1,i + 1 log(Ei )]37

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

> summary(regPoisson(N~X1+offset(log(E))))
Coefficients:
Estimate Std. Error z value Pr(>|z|)(Intercept) -2.586986
0.038389 -67.388<2e-16 ***X1-0.004364
0.004103 -1.064
0.287--Signif. codes: 0 â***â 0.001 â**â 0.01 â*â 0.05 â.â 0.1 â â 1(Dispersion parameter for poisson family taken to be 1)
Null deviance: 12680
Residual deviance: 12679
AIC: 16354

on 49999
on 49998

degrees of freedom
degrees of freedom

Number of Fisher Scoring iterations: 638

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Régression de Poisson multiple
On peut (bien entendu) régresser sur plusieurs variables explicatives> model1=regPoisson(nbre~zone+as.factor(puissance)+agevehicule++ageconducteur+carburant+offset(log(exposition)),data=baseFREQ)> summary(model1)
Coefficients:(Intercept)
zoneB
zoneC
zoneD
zoneE
zoneF
as.factor(puissance)5
as.factor(puissance)6
as.factor(puissance)7
as.factor(puissance)8
as.factor(puissance)9

Estimate Std. Error z value Pr(>|z|)-2.546970
0.122723 -20.754 < 2e-16 ***0.011487
0.097559
0.118
0.90630.196208
0.077090
2.545
0.0109 *0.403382
0.078788
5.120 3.06e-07 ***0.594872
0.079207
7.510 5.90e-14 ***0.684673
0.143612
4.768 1.87e-06 ***0.135072
0.081393
1.659
0.0970 .
0.161305
0.079692
2.024
0.0430 *0.164168
0.079039
2.077
0.0378 *0.122254
0.110876
1.103
0.27020.1819780.123996
1.468
0.142239

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

as.factor(puissance)10 0.254358
0.119777
2.124
0.0337 *as.factor(puissance)11 0.001156
0.170163
0.007
0.9946as.factor(puissance)12 0.243677
0.223207
1.092
0.2750as.factor(puissance)13 0.513950
0.284159
1.809
0.0705 .
as.factor(puissance)14 0.582564
0.295482
1.972
0.0487 *as.factor(puissance)15 0.173748
0.383322
0.453
0.6504
agevehicule
0.001467
0.004191
0.350
0.7264ageconducteur-0.008844
0.001658 -5.335 9.58e-08 ***carburantE-0.201780
0.049265 -4.096 4.21e-05 ***--Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1(Dispersion parameter for poisson family taken to be 1)
Null deviance: 12680
Residual deviance: 12524
AIC: 16235

on 49999
on 49980

degrees of freedom
degrees of freedom

Number of Fisher Scoring iterations: 640

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Effets marginaux, et élasticité
Les effets marginaux de la variable k pour l’individu i sont donnés par∂E(Yi |X i )∂ exp[X 0i β]== exp[X 0i β] · βk
∂Xi,k
∂Xk
b · βbk
estimés pas exp[X 0i β]Par exemple, pour avoir l’effet marginal de la variable ageconducteur,> coef(model1)[19]
ageconducteur-0.008844096> baseFREQ[1:4,]
nocontrat exposition zone puissance agevehicule ageconducteur
1
27
0.87
C
7
0
56
2
115
0.72
D
5
0
45
3
121
0.05
C
6
0
37
4
142
0.90
C
10
10
42

...
41

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

> effet19=predict(model1,type="response")*coef(model1)[19]> effet19[1:4]
1
2
3
4-5.265678e-04 -4.690555e-04 -3.569767e-05 -6.846838e-04

On peut aussi calculer les effets marginaux moyens de la variable k, Y · βbk> mean(predict(model1,type="response"))*coef(model1)[19]
ageconducteur-0.0003403208

Autrement dit, en vieillissant d’un an, il y aura (en moyenne) 0.00034 accident demoins, par an, par assuré.
Ici, on utilise des changements en unité (∂Xi,k ), mais il est possible détudier
l’impact de changement en proportion. Au lieu de varier d’une utilité, on va
considérer un changement de 1%.

42Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Interprétation, suite
Ici, les zones prennent les valeurs A B C D E ou F, selon la densité en nombred’habitants par km2 de la commune de résidence (A ="1-50", B="50-100",
C="100-500", D="500-2,000", E="2,000-10,000", F="10,000+").
Dans la sortie, nous avons obtenu
Coefficients:
carburantE

Estimate Std. Error z value Pr(>|z|)-0.201780
0.049265 -4.096 4.21e-05 ***i.e.
> coefficients(model1)["carburantE"]
carburantE-0.2017798

Autrement dit, à caractéristiques identiques, un assuré conduisant un véhicule
essence a une fréquence de sinistres presque 20% plus faible qu’un assuré
conduisant un véhicule diesel,
43Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

> exp(coefficients(model1)["carburantE"])
carburantE
0.8172749

Si on regarde maintenant la zone, la zone A a disparu (c’est la zone de référence).
On notera que la zone B n’est pas significativement différente de la zone A.
Coefficients:
zoneB
zoneC
zoneD
zoneE
zoneF

Estimate Std. Error z value Pr(>|z|)0.011487
0.097559
0.118
0.90630.196208
0.077090
2.545
0.0109 *0.403382
0.078788
5.120 3.06e-07 ***0.594872
0.079207
7.510 5.90e-14 ***0.684673
0.143612
4.768 1.87e-06 ***Notons que l’on pourrait choisir une autre zone de référence> baseFREQ$zone=relevel(baseFREQ$zone,"C")

Si on refait la régression, on trouve que tous les zones sont disinctes de la zone C.
44Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

> model1=regPoisson(nbre~zone+as.factor(puissance)+agevehicule++ageconducteur+carburant+offset(log(exposition)),data=baseFREQ)> summary(model1)
Coefficients:
Estimate Std. Error z value Pr(>|z|)(Intercept)-2.350761
0.112384 -20.917 < 2e-16 ***zoneA-0.196208
0.077090 -2.545 0.010921 *zoneB-0.184722
0.086739 -2.130 0.033202 *zoneD
0.207174
0.064415
3.216 0.001299 **zoneE
0.398664
0.064569
6.174 6.65e-10 ***zoneF
0.488465
0.135765
3.598 0.000321 ***as.factor(puissance)5
0.135072
0.081393
1.659 0.097017 .
as.factor(puissance)60.161305
0.079692
2.024 0.042959 *as.factor(puissance)7
0.164168
0.079039
2.077 0.037798 *as.factor(puissance)80.122254
0.110876
1.103 0.270195
as.factor(puissance)90.1819780.123996
1.468 0.142209
as.factor(puissance)10 0.254358
0.119777
2.124 0.033704 *as.factor(puissance)11 0.001156
0.170163
0.007 0.994578
as.factor(puissance)12 0.243677
0.223207
1.092 0.274960
as.factor(puissance)13 0.513950
0.284159
1.809 0.070502 .

45

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

as.factor(puissance)14 0.582564
0.295482
1.972 0.048659 *as.factor(puissance)15 0.173748
0.383322
0.453 0.650356
agevehicule
0.001467
0.004191
0.350 0.726375
ageconducteur-0.008844
0.001658 -5.335 9.58e-08 ***carburantE-0.201780
0.049265 -4.096 4.21e-05 ***--Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1(Dispersion parameter for poisson family taken to be 1)
Null deviance: 12680
Residual deviance: 12524
AIC: 16235

on 49999
on 49980

degrees of freedom
degrees of freedom

Number of Fisher Scoring iterations: 6

Pareil pour la zone D. En revanche, si la modalité de référence devient la zone E,
on note que la zone F ne se distingue pas,
> baseFREQ$zone=relevel(baseFREQ$zone,"E")46

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

puisque l’on a
Coefficients:(Intercept)
zoneD
zoneC
zoneA
zoneB
zoneF

Estimate Std. Error z value Pr(>|z|)-1.952098
0.112128 -17.410 < 2e-16 ***-0.191490
0.066185 -2.893 0.00381 **-0.398664
0.064569 -6.174 6.65e-10 ***-0.594872
0.079207 -7.510 5.90e-14 ***-0.583385
0.088478 -6.594 4.29e-11 ***0.089801
0.135986
0.660 0.50902

On pourrait être tenté de regroupe A et B, et E et F. En effet, les tests deStudent suggèrent des regroupement (pour chacune des paires). Pour faire un
regroupement des deux paires, on fait un test de Ficher,> baseFREQ$zone=relevel(baseFREQ$zone,"C")> model1=regPoisson(nbre~zone+as.factor(puissance)+agevehicule++ageconducteur+carburant+offset(log(exposition)),data=baseFREQ)> library(car)> linearHypothesis(model1,c("zoneA=zoneB","zoneE=zoneF"))

47

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Linear hypothesis test
Hypothesis:zoneA - zoneB = 0
zoneE - zoneF = 0
Model 1: restricted model
Model 2: nbre ~ zone + as.factor(puissance) + agevehicule + ageconducteur +carburant + offset(log(exposition))
Res.Df Df Chisq Pr(>Chisq)
1 49982
2 49980 2 0.4498
0.7986

On peut accepter (avec une telle p-value) un regroupement. Construisons cettenouvelle variable> baseFREQ$zonesimple=baseFREQ$zone> baseFREQ$zonesimple[baseFREQ$zone%in%c("A","B")]="A"> baseFREQ$zonesimple[baseFREQ$zone%in%c("E","F")]="E"48

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

(pour des raisons techniques, il faut le nom existe déjà, i.e on ne peut pas créer laclasse AB).
> model1=regPoisson(nbre~zonesimple+as.factor(puissance)+agevehicule++ageconducteur+carburant+offset(log(exposition)),data=baseFREQ)> summary(model1)
Coefficients:(Intercept)
zonesimpleD
zonesimpleC
zonesimpleA

Estimate Std. Error z value Pr(>|z|)-1.943339
0.111282 -17.463 < 2e-16 ***-0.201446
0.064289 -3.133 0.00173 **-0.408524
0.062658 -6.520 7.03e-11 ***-0.599785
0.066128 -9.070 < 2e-16 ***49

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Lissage et modèles non-paramétriques
Dans le modèle où seul l’âge du conducteur intervient, on a
> model2=regPoisson(nbre~ageconducteur+offset(log(exposition)),data=baseFREQ)> summary(model2)Call:glm(formula = formula, family = poisson(link = log), data = data)
Deviance Residuals:
Min
1Q
Median-0.5063 -0.3504 -0.2601

3Q-0.1412Max
13.3387

Coefficients:
Estimate Std. Error z value Pr(>|z|)(Intercept)-2.163061
0.077380 -27.95 < 2e-16 ***ageconducteur -0.009891
0.001635-6.05 1.45e-09 ***--Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

50Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

(Dispersion parameter for poisson family taken to be 1)
Null deviance: 12680
Residual deviance: 12642
AIC: 16317

on 49999
on 49998

degrees of freedom
degrees of freedom

Number of Fisher Scoring iterations: 6

i.e. la fréquence annuelle de sinistre suit une loi de Poisson de paramètreb = exp[βb0 + βb1 x] pour un assuré d’âge x.
λage=18:80
lambda=predict(model2,newdata=data.frame(ageconducteur=age,exposition=1),type="response")plot(age,lambda,type="l",lwd=2,col="blue")abline(h=sum(N)/sum(E))

51

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Mais est-ce raisonnable de supposer une décroissance exponentielle en fonction del’âge ?
On peut tenter une régression paramétrique en prenant l’âge comme variable
catégorielle> model2b=regPoisson(nbre~as.factor(ageconducteur)+offset(log(exposition)),data=baseFREQ)> summary(model2b)52

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Call:glm(formula = formula, family = poisson(link = log), data = data)
Deviance Residuals:
Min
1Q
Median-0.6949 -0.3425 -0.2521

3Q-0.1381

Max
12.8912

Coefficients:(Intercept)
as.factor(ageconducteur)19
as.factor(ageconducteur)20
as.factor(ageconducteur)21
as.factor(ageconducteur)22
as.factor(ageconducteur)23
as.factor(ageconducteur)24
as.factor(ageconducteur)25

Estimate Std. Error z value Pr(>|z|)-14.2522
233.9921 -0.061
0.951
12.5148233.9923
0.053
0.957
12.8309
233.9922
0.055
0.956
12.1014233.9923
0.052
0.959
12.2352
233.9922
0.052
0.958
12.6185
233.9922
0.054
0.957
12.1043233.9922
0.052
0.959
12.0647233.9922
0.052
0.959

...
as.factor(ageconducteur)97-1.9218

1227.3574-0.0020.999

53

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

as.factor(ageconducteur)98
as.factor(ageconducteur)99
as.factor(ageconducteur)100-2.0504-1.7031-1.8760

2116.3379
926.5006
2116.3379-0.001-0.002-0.001

0.999
0.999
0.999

(Dispersion parameter for poisson family taken to be 1)
Null deviance: 12680
Residual deviance: 12517
AIC: 16352

on 49999
on 49918

degrees of freedom
degrees of freedom

Number of Fisher Scoring iterations: 1454

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

On peut utiliser des méthodes de lissage pour estimer s(·) telle que
b = exp[βb0 + sb(x)]λ55Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Lissage et modèles non-paramétriques56Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Lissage et modèles non-paramétriques57

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Le lissage par splines ?
Une première idée pourrait ètre de faire des modèles linéaires par morceaux etcontinus. Par exemple, on peut introduire une rupture en x = 30

K=30model2spline=regPoisson(nbre~bs(ageconducteur,knots=c(K),degree=1)+offset(log(exposition)),data=baseFREQ)
lambda=predict(model2spline,newdata=data.frame(ageconducteur=age,exposition=1),type="respo
plot(age,lambda,type="l",lwd=2,col="blue")abline(h=sum(N)/sum(E),col="red",lty=2)
abline(v=K)

58

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

ou encore en x = 25 et x = 35,
K=c(25,35)59

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Mais le lissage par des fonctions linéaires par morceaux garantie la continué de lafonction, pas de sa dérivée.
On peut introduire les blueb-splines de degré n, définies, sur une grille de mpoints ti (appelés noeuds) dans [0, 1], avec 0 ≤ t0 ≤ t1 ≤ . . . ≤ tm ≤ 1, par60

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

récurrence sur le degré inférieur par une relation de la forme1 si t 6 t < tjj+1
bj,0 (t) :=
0
sinon
et

t − tjtj+n+1 − t
bj,n (t) :=bj,n−1 (t) +bj+1,n−1 (t).
tj+n − tj
tj+n+1 − tj+1

Numériquement
B=function(x,j,n,K){b=0
a1=a2=0if(((K[j+n+1]>K[j+1])&(j+n<=length(K))&(n>0))==TRUE){a2=(K[j+n+1]-x)/(K[j+n+1]-K[j+1])*B(x,j+1,n-1,K) }if(((K[j+n]>K[j])&(n>0))==TRUE){a1=(x-K[j])/(K[j+n]-K[j])*B(x,j,n-1,K)}
if(n==0){ b=((x>K[j])&(x<=K[j+1]))*1 }if(n>0){b=a1+a2}
return(b)}61Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Pour les splines de degré 1, avec un seul noeud on a comme base de fonctions
u=seq(0,1,by=.01)
plot(u,B(u,1,1,c(0,.4,1,1)),lwd=2,col="red",type="l",ylim=c(0,1))
lines(u,B(u,2,1,c(0,.4,1,1)),lwd=2,col="blue")abline(v=c(0,.4,1),lty=2)

et les splines de degré 2, avec un seul noeud on a comme base de fonctions
62

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

u=seq(0,1,by=.01)
plot(u,B(u,1,2,c(0,0,.4,1,1,1)),lwd=2,col="red",type="l",ylim=c(0,1))
lines(u,B(u,2,2,c(0,0,.4,1,1,1)),lwd=2,col="blue")lines(u,B(u,3,2,c(0,0,.4,1,1,1)),lwd=2,col="green")
abline(v=c(0,.4,1),lty=2)

On peut uiliser des b-splines de degré 2 pour faire du lissage,
K=30

63

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

model2spline=regPoisson(nbre~bs(ageconducteur,knots=c(K),degree=2)+offset(log(exposition)),data=baseFREQ)
lambda=predict(model2spline,newdata=data.frame(ageconducteur=age,exposition=1),type="respo
plot(age,lambda,type="l",lwd=2,col="blue")abline(h=sum(N)/sum(E),col="red",lty=2)
abline(v=K)

avec un noeud, ou deux,64Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

K=c(25,35)
model2spline=regPoisson(nbre~bs(ageconducteur,knots=c(K),degree=2)+offset(log(exposition)),data=baseFREQ)

Pour trouver où mettre le noeud de manière optimale, on peut fixer un critère(classique) comme minimiser la somme des carrés des erreurs (de Pearson, e.g.)ou le critère d’Akaike
65Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

SSR=AKAIKE=matrix(NA,60)
for(i in 1:60){
model2spline=regPoisson(nbre~bs(ageconducteur,knots=18+i,degree=2)+offset(log(exposition))
SSR[i]=sum(residuals(model2spline,type="pearson")^2)
AKAIKE[i]=AIC(model2spline)}

Pour la somme des carrés des résidus, on a une rupture vers 29 ans,plot(18+1:60,SSR,type="b",col="blue")66Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013et pour le critère d’Akaike, on a une rupture vers 35 ans,plot(18+1:60,AKAIKE,type="b",col="red")67Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Remarque : Par défaut, R utilise des splines cubiques, de degré 3.

68