Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Actuariat IARD - ACT2040
Partie 2 - régression logistiqueet arbres de régression (Y ∈ {0, 1})
Arthur Charpentiercharpentier.arthur@uqam.ca
http ://freakonometrics.hypotheses.org/

Automne 20131

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Modèlisation d’une variable dichotomiqueRéférences : Frees (2010), chapitre 11 (p 305-342), Greene (2012), sections 17.2
et 17.3 (p 863-715), de Jong & Heller (2008), chapitre 7
Remarque : la régression logistique est un cas particulier des modèles GLM,
avec une loi binomiale et une fonction de lien logit.
logit = function(formula, lien="logit", data=NULL) {glm(formula,family=binomial(link=lien),data) }Modèles introduits par Gaddum (1933), Bliss (1935), ou Berkson (1944).

2

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La loi binomiale B(π)
Soient Yi des variables aléatoires prenant les valeurs yi ∈ {0, 1}, avec probabilité
1 − πi et πi (respectivement), i.e.
P(Yi = yi ) = πiyi [1 − πi ]1−yi , pour yi ∈ {0, 1}.
avec πi ∈ [0, 1]. Aussi, P(Yi = 1) = πi et P(Yi = 0) = 1 − πi . Alors E(Yi ) = πi etVar(Yi ) = π[1 − πi ].

En statistique, on suppose que πi = π ∀i ∈ {1, 2, · · · , n}, et on cherche à estimer π

3Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Maximum de vraisemblance et méthode des moments
La fonction de vraisemblance, pour un échantillon y1 , · · · , yn s’écrit
L(π; y) =nY

P(Yi = yi ) =

i=1

n
Y

π yi [1 − π]1−yi

i=1

et la log-vraisemblance est
log L(π; y) =nXyi log[π] + (1 − yi ) log[1 − π]i=1

La condition du premier ordre estn∂ log L(π; y) X yi
1 − yi
=
−
=0
∂π
π
1−π
i=1

4

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Intervalle de confiance pour π
On a un intervalle de confiance asymptotique car on sait que l’estimateur devraisemblance est asymptotiquement efficace,√Ln(π − πb) → N (0, I −1 (π))
où I(π) est l’information de Fisher, i.e.
I(π) =

1π[1 − π]d’où un intervalle de confiance approché (à 95%) pour π de la forme


p
1.96
π
b[1 − π
b] .
π
b± √
n
On peut aussi construire un intervalle de confiance, par le théorème central
limite, car π
b = X. On sait que√ X − E(X) L→ N (0, 1),np
Var(X)
5

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

avec ici X = πb, E(X) = π et Var(X) = π(1 − π), i.e. un intervalle de confiance
est obtenu par l’approximation
√

np

π
b−π
π
b[1 − πb]L

→ N (0, 1),d’où un intervalle de confiance (à 95%) pour π de la forme


p
1.96
π
b± √
π
b[1 − πb] .
nAvec la méthode du score

√

n

π
b−π
≤ 1.96π[1 − π]i.e.(en élevant au carré)p
1.96 1.962 + 4nb
π [1 − πb]2nb
π + 1.96±2(n + 1.962 )2(n + 1.962 )2

6

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Mise en oeuvre pratique
Considérons l’échantillon suivant, {y1 , · · · , yn }> set.seed(1)> n=20> (Y=sample(0:1,size=n,replace=TRUE))[1] 0 0 1 1 0 1 1 1 1 0 0 0 1 0 1 0 1 1 0 1

Comme Yi ∼ B(π), avec π = E(Y ), l’estimateur de la méthde des moments est
πb = y, soit ici> mean(X)[1] 0.55

On peut faire un test de H0 : π = π? contre H1 : π 6= π? (par exemple 50%)
On peut utiliser le test de Student,
T =

√

np

π
b − π?
π? (1 − π? )7

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

qui suit, sous H0 une loi de Student à n degrés de liberté.

0.3
0.2
0.1
0.0

dt(u, df = n)

0.4

> (T=sqrt(n)*(pn-p0)/(sqrt(p0*(1-p0))))[1] 0.4472136> abs(T)<qt(1-alpha/2,df=n)[1] TRUE

●

−3

−2

−1

0

1

2

3

On est ici dans la région d’acceptation du test.
8Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

On peut aussi calculer la p-value, P(|T | > |tobs |),

0.3
0.2
0.1
0.0

dt(u, df = n)0.4

> 2*(1-pt(abs(T),df=n))[1] 0.6595265

●

−3

−2

−1

0

1

2

3

Le maximum de vraisemblance est important car il permet de faire des tests
statistiques.
9Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Test de Wald l’idée est d’étudier la différence entre πb et π? . Sous H0 ,(bπ − π? ) 2 L 2
T = n −1→ χ (1)I (π? )
Test du rapport de vraisemblance l’idée est d’étudier la différence entre
b et log L(θ? ). Sous H0 ,log L(θ)!
log L(θ? ) L 2
T = 2 log
→ χ (1)blog L(θ)
Test du score l’idée est d’étudier la différence entre

T =

n
X

1
∂ log fπ? (xi )
n i=1
∂π

∂ log L(π? )et 0. Sous H0 ,
∂π

!2
L

→ χ2 (1)10Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Mise en oeuvre pratique
Pour l’estimateur de maximum de vraisemblance, on peut faire les calculs, ou lefaire numériquement. On commence par tracer la fonction θ 7→ log L(θ)> p=seq(0,1,by=.01)> logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}> plot(p,Vectorize(logL)(p),type="l",col="red",lwd=2)

11

−20
−30
−40
−50

Vectorize(logL)(p)Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

0.0

0.2

0.4

0.6

0.8

1.0

p

On peut trouver numériquement le maximum de log L,> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}> pml=optim(fn=neglogL,par=p0,method="BFGS")
> pml
$par
[1] 0.5499996

12

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013$value[1] 13.76278

i.e. on retrouve ici πb = y.
Maintenant, on peut tester H0 : π = π? = 50% versus H1 : π 6= 50%. Pour le test
de Wald, on peut commencer par calculer nI(θ? ) ou ,> nx=sum(X==1)> f = expression(nx*log(p)+(n-nx)*log(1-p))> Df = D(f, "p")> Df2 = D(Df, "p")> p=p0=0.5> (IF=-eval(Df2))[1] 80On peut d’ailleurs comparer cette valeur avec la valeur théorique, car
I(π)−1 = π(1 − π)
> 1/(p0*(1-p0)/n)[1] 8013

−15.0

−14.0

●

−16.0

Vectorize(logL)(p)−13.0

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

0.0

0.2

0.4

0.6

0.8

1.0

p

La statistique du test de Wald est ici> pml=optim(fn=neglogL,par=p0,method="BFGS")$par> (T=(pml-p0)^2*IF)[1] 0.199997

que l’on va chercher à comparer avec le quantile de la loi du χ2 ,14

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

> T<qchisq(1-alpha,df=1)[1] TRUE

i.e. on est dans la région d’acceptation du test. De manière duale, on peut
calculer la p-value du test,> 1-pchisq(T,df=1)[1] 0.6547233

Donc on va accepter H0 .

15

1.5
1.0
0.5
0.0

dchisq(u, df = 1)

2.0

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

●

0

1

2

3

4

5

6

Pour le test du rapport de vraisemblance, T est ici> (T=2*(logL(pml)-logL(p0)))[1] 0.2003347

16

−15.0

−14.0

●

−16.0

Vectorize(logL)(p)−13.0

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

0.0

0.2

0.4

0.6

0.8

1.0

p

Là encore, on est dans la région d’acceptation,> T<qchisq(1-alpha,df=1)[1] TRUE

Enfin, on peut faire le test du score. Le score se calcule facilement,
17Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

> nx=sum(X==1)> f = expression(nx*log(p)+(n-nx)*log(1-p))> Df = D(f, "p")> p=p0> score=eval(Df)La statistique de test est alors> (T=score^2/IF)[1] 0.218

−15.0

−14.0

●

−16.0

Vectorize(logL)(p)−13.0

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

0.0

0.2

0.4

0.6

0.8

1.0

p

Là encore, on est dans la région d’acceptation,> T<qchisq(1-alpha,df=1)[1] TRUE

19

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Prise en compte de l’hétérogénéité
Supposons que πi soit fonction de variables explicatives X i , e.g.
πi = E(Yi |X i ).
Le modèle classique est le modèle linéaire, E(Yi |X i ) = X 0i β
On peut faire une estimation par moindres carrés pour estimer β,( n)Xb = argmin
β[Yi − X 0 β]2 , β ∈ Rki

i=1

Problème : πi ∈ [0, 1] et X 0i β ∈ R.
Problème : modèle non homoscédastique, en effet 1 − X 0 β avec probabilité πiiεi = −X 0 β avec probabilité 1 − πiii.e. Var(ε) = X 0i β · [1 − xX 0i β].
20Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Prise en compte de l’hétérogénéité
baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv",header=TRUE,sep=",")
avocat=(baseavocat$ATTORNEY==1)
sexe=rep("H",length(avocat)); sexe[baseavocat$CLMSEX==2]="F"marital=rep("M",length(avocat)); marital[baseavocat$MARITAL==2]="C"
marital[baseavocat$MARITAL==3]="V"; marital[baseavocat$MARITAL==4]="D"ceinture=(baseavocat$SEATBELT==1);
age=baseavocat$CLMAGE
cout=baseavocat$LOSS
base=data.frame(avocat,sexe,marital,ceinture,age,cout)

La base ici ressemble à> tail(base,3)
avocat sexe marital ceinture age cout
1338 FALSE
F
M
TRUE 39 0.099
1339
TRUE
F
C
TRUE 18 3.277
1340 FALSEF
C
TRUE 30 0.688

21Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Utilisation de la cote (odds)
πi ∈ [0, 1] alors que X 0i β ∈ R, donc on ne peut pas supposer
πi = E(Yi |X i ) = X 0i βπi
∈ [0, ∞].
oddsi =
1−πou

log(oddsi ) = log

πi
1−π


∈ R.

On appelle logit cette transormation,logit(πi ) = log

πi
1−π



22

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Utilisation de la cote (odds)
On va supposer que
logit(πi ) = X 0i βou
−1πi = logit(X 0i β) =exp[X 0i β].
01 + exp[X i β]Remarque : si πi est faible (πi ∼ 0) alors logit−1 ∼ exp, i.e.
πi ∼ exp(X 0i β)

23Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La fonction logistique
F (x) =

exp[x], ∀x ∈ R.
1 + exp[x]24Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La fonction logistique
F (x) =

exp[2x], ∀x ∈ R.
1 + exp[2x]25Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La fonction logistique
F (x) =

exp[−x], ∀x ∈ R.
1 + exp[−x]26Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La fonction logistique
F (x) =

exp[1+x], ∀x ∈ R.
1 + exp[1+x]27Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La régression logistique
La log-vraisemblance est ici
log L(β) =nXyi log(πi )+(1−yi ) log(1−πi ) =

i=1

n
X

yi log(πi (β))+(1−yi ) log(1−πi (β))

i=1

La résolution se fait en écrivant les conditions du premier ordre. Pour cela, onécrit le gradient ∇ log L(β) = [∂ log L(β)/∂βk ], oùnX∂ log L(β)
yi ∂πi (β) 1 − yi ∂πi (β)
==
−
∂βk
π
(β)
∂β
πi (β) ∂βkk
i=1 i

or compte tenu de la forme de πi (β),∂πi (β)= πi (β)[1 − πi (β)]Xk,i∂βkon obtientnX∂ log L(β)==
Xk,i [yi − πi (β)]
∂βk
i=1

28

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Newton Raphson - Descente de gradientRemarque Pour rechercher le zéro d’une fonction (réelle), i.e. f (x) = 0 : on sedonne x0 , et on posef (xk−1 )xk = xk−1 − 0
f (xk−1 )Remarque Pour rechercher le maximum d’un fonction dérivable, on recherche lezéro de la dérivée, i.e. h0 (x) = 0 : on se donne x0 , et on posef 0 (xk−1 )xk = xk−1 − 00
f (xk−1 )Il faut trouver numériquement le maximum de log L(β).
L’algorithme standard est
1. partir d’une valeur initiale β 0
2. poser β k = β k−1 − H(β k−1 )−1 ∇ log L(β k−1 )29Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013où ∇ log L(β) est le gradient, et H(β) la matrice Hessienne (on parle parfois de
Score de Fisher).
Ici, le terme générique de la matrice Hessienne estn∂ 2 log L(β) X=
Xk,i X`,i [yi − πi (β)]∂βk ∂β`i=1

PosonsΩ = [ωi,j ] = diag(bπi (1 − πbi ))alors le gradient est∂ log L(β)= X 0 (y − π)
∇ log L(β) =∂βalors que la matrice Hessienne estimation∂ 2 log L(β)0
H(β) =
=
−X
ΩX
0
∂β∂β
30

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Aussi, l’algorithme de Newton Raphson s’écrit
β k = (X 0 ΩX)−1 X 0 ΩZ où Z = Xβ k−1 + X 0 Ω−1 y − π),On reconnait une régression pondérée itérée,β k = argmin {(Z − Xβ)0 Ω(Z − Xβ).}D’un point de vue computationnel, l’estimation se fait de la manière suivante :X=cbind(rep(1,nrow(base)),base$cout)
Y=base$avocat*1

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,
31Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

for(s in 1:20){
pi=exp(X%*%beta)/(1+exp(X%*%beta))
gradient=t(X)%*%(Y-pi)
omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi))hesienne=-t(X)%*%omega%*%X
beta=beta-solve(hesienne)%*%gradient}Ici, on obtient les valeurs suivantes
[,1]
[,2]
[1,] 0.50034948 0.001821554[2,] -0.05689403 0.015434935[3,] -0.19941087 0.055570172[4,] -0.45629828 0.140490153[5,] -0.73627462 0.253154397[6,] -0.85830164 0.309712845[7,] -0.87192186 0.316615292[8,] -0.87207800 0.316697581

i.e. on converge (et vite).
32Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Comportement asymptotique (et intervalle de confiance)Pb→
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
> pi=exp(X%*%beta)/(1+exp(X%*%beta))+ omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi))+ hesienne=-t(X)%*%omega%*%X+ sqrt(-diag(solve(hesienne)))[1] 0.09128379 0.02902411

33

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

On retrouve toutes ces valeurs en utilisant> reglogit=logit(Y~0+X)
Coefficients:
Estimate Std. Error z value Pr(>|z|)X1 -0.87208
0.09128 -9.553<2e-16 ***X2 0.31670
0.02902 10.912<2e-16 ***--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1857.6
Residual deviance: 1609.4
AIC: 1613.4

on 1340
on 1338

degrees of freedom
degrees of freedom

Number of Fisher Scoring iterations: 8
Message d’avis :glm.fit: fitted probabilities numerically 0 or 1 occurred34

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Test d’hypothèse et significativité
Considérons un test H0 : βk = 0 contre H1 : βk 6= 0 (classique dans le modèle
linéaire).
Compte tenu du comportement asymptotiquement Gaussien des estimateurs, ona automatiquement un test de significativité des coefficients : sous H0 ,βbkq

∼ N (0, 1).

d βbk )
Var(
Remarque : ce test correspond à un test de Wald (d’où la notation z sous R)

35

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Régression linéaire sur une variable catégorielle
Supposons que l’on régresse sur X, une variable binaire, i.e X ∈ {0, 1}, avec
P(Y = 1|X) = β0 + β1 X
qui prend les valeurs β0 si x = 0 (modalité de référence) et β0 + β1 si x = 1.
La log-vraisemblance devientX
X
log L(β0 , β1 ) =yi log[β0 ] +yi log[β0 + β1 ]i,xi =0+i,xi =1X(1 − yi ) log[1 − β0 ] +i,xi =0X(1 − yi ) log[1 − β0 + β1 ]i,xi =1

et les conditions du premier ordre donnent
X
X
X
yi
yi
yi
xi =0
i =0
i =1
c0 = xX
c1 = xX
β
et β
− X
1
1
1
xi =0

xi =0

xi =0

36

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Régression linéaire sur une variable catégorielle
Sur nos données, on peut considérer le sexe comme variable explicative.
> xtabs(~avocat + sexe, data=base)sexe
avocatF
H
FALSE 390 265
TRUE 352 333

Si on fait la régression linéaire
> reglm = lm(avocat ~ sexe, data=base)> summary(reglm)
Coefficients:
Estimate Std. Error t value Pr(>|t|)(Intercept) 0.47439
0.01830
25.92 < 2e-16 ***sexeH
0.08246
0.02740
3.01 0.00266 **37

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Régression linéaire sur une variable catégorielle
Ces paramètres correspondent aux fréquences empiriques,> attach(base)> sum((avocat==1)&(sexe=="F"))/sum((sexe=="F"))[1] 0.4743935> sum((avocat==1)&(sexe=="H"))/sum((sexe=="H"))+ sum((avocat==1)&(sexe=="F"))/sum((sexe=="F"))[1] 0.08246266

38

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Régression logistique sur une variable catégorielle
Avec un modèle logistique, on obtient les mêmes prédictions,> reglogit = logit(avocat ~ sexe, data=base)> summary(reglogit)
Coefficients:
Estimate Std. Error z value Pr(>|z|)(Intercept) -0.10252
0.07352 -1.394 0.16319
sexeH
0.330930.110372.998 0.00271 **--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1> exp(reglogit$coefficients[1])/+ (1+exp(reglogit$coefficients[1]))(Intercept)
0.4743935> exp(sum(reglogit$coefficients[1:2]))/+ (1+exp(sum(reglogit$coefficients[1:2])))[1] 0.5568562

39

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Modèle alternatif, le modèle probit
Soit H une fonction de répartition d’une variable rélle, H : R 7→ [0, 1]. On
suppose ici
πi = H(X 0i β) ou X 0i β = H −1 (πi )– si H est la fonction de répartition de la loi logistique, on retrouve le modèle
logit,– si H est la fonction de répartition de la loi N (0, 1), on obtient le modèle probit,
Remarque : probit signifie probability unit.
On suppose qu’il existe une variable latente Yi? non observée, telle que
πi = P(Yi = 1) = P(Yi? > 0)On suppose que Yi? = X 0i β + i , où i est une erreur de loi H. Alors
πi = P(Yi = 1) = P(Yi? > 0) = P(i > −X 0i β) = 1 − H(−X 0i β)40Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013alors que1 − πi = P(Yi = 0) = P(Yi? ≤ 0) = P(i ≤ −X 0i β) = H(−X 0i β)
Si H = Φ, alors H est symmétrique par rapport à 0, et
πi = P(Yi = 1) = 1 − Φ(−X 0i β) = Φ(X 0i β),ou encore ηi = X 0i β = Φ−1 (πi ).
ex, on a une loi logistique, qui est aussi symmétrique par rapport
Si H(x) =1 + exà 0, etπi = P(Yi = 1) = 1 − H(−X 0i β) = H(X 0i β),ou encore ηi = X 0i β = H −1 (πi ).
Remarque : pourquoi prendre un seuil nul ? sinon le modèle ne serait pas
identificable. De même si on ne fixe pas la variance des résidus. On prend donc
(arbitrairement) s = 0 et σ 2 = 1.
41Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La fonction logistique
Z x

 2
1
z
Φ(x) = √exp −dz, ∀x ∈ R.
2
2π −∞

42

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

La fonction logistique
√ Z x


2
3z
3
√
exp − 2 dz, ∀x ∈ R.
Φ(x) =2π
π 2π −∞

43

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Qualité de l’ajustement
Pour juger de la validité du modèle, on peut considérer plusieurs quantités.
On peut définir la déviance, qui compare la log-vraisemblance du modèle et cellle
du modèle parfait,nX
D = −2[Yi log[bπi ] + (1 − Yi ) log[1 − π
bi ]]
i=1

Plus la variance est faible, meilleur est le modèle .
On peut aussi considérer la statistique de Pearson
χ2 =

n
X
(yi − π
bi )2
i=1

π
bi [1 − πbi ]44Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Étude des résidus
A l’aide de la relation précédante, on peut définir
εbi = pYi − πbi
πbi [1 − πbi ]On parlera de résidus de Pearson.
Il existe aussi les résidus de déviance.
1εbi = ±[Yi log[bπi ] + (1 − Yi ) log[1 − πbi ]] 2
avec un signe positif si Yi ≥ πbi , négatif sinon.

45Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Visualisation des résidus> reglogit = logit(avocat ~ cout+sexe, data=base)> E = residuals(reglogit,type=‘‘response’’)46Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Visualisation des résidus> reglogit = logit(avocat ~ cout+sexe, data=base)> E = residuals(reglogit,type=‘‘deviance’’)47Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Prédiction avec un modèle logistiquebSous R, la prédiction sera (par défaut) une prédiction du score, i.e. x0 β.
> predict(reglogit,newdata=data.frame(cout=10,sexe="H"))12.530417

Il faut alors utiliser la loi logistique pour prévoir π,> exp(predict(reglogit,newdata=data.frame(cout=10,sexe="H")))/+ (1+exp(predict(reglogit,newdata=data.frame(cout=10,sexe="H"))))
1
0.9262468

ou utiliser directement> predict(reglogit,newdata=data.frame(cout=10,sexe="H"),type="response")
1
0.9262468

48

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Prédiction avec un modèle logistique49

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Prédiction avec un modèle logistique

Pour juger de la qualité du modèle, plusieurs outils existent.

50Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Courbe ROC, receiver operating characteristic
On dispose d’observations Yi prenant les valeurs 0 (on dira négatives) ou 1 (ondira positives).
On a construit un modèle qui prédit πi . En se fixant un seuil s, si π
bi ≤ s, on
prédit Ybi = 0, et si π
bi > s, on prédit Ybi = 1.
Le modèle sera bon si les positifs sont prédits positifs, et les négatifs sont prédits
négatifs.
Le choix du seuil s permettra de minimiser soit les faux positifs, soit les fauxnégatifs.
Y =1 Y =0
Ybi = 1
Ybi = 0

TP

FP

FN

TN

où TP sont les True Positive, FP les False Positive, etc.
51Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Courbe ROC, receiver operating characteristic
S=predict(reglogit,type="response")plot(0:1,0:1,xlab="False Positive Rate",ylab="True Positive Rate",cex=.5)
for(s in seq(0,1,by=.01)){
Ps=(S>s)*1
FP=sum((Ps==1)*(avocat==0))/sum(avocat==0)
TP=sum((Ps==1)*(avocat==1))/sum(avocat==1)
points(FP,TP,cex=.5,col="red")}52

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Courbe ROC, receiver operating characteristic53

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Courbe ROC, receiver operating characteristic
En reliant les points, on obtient la courbe ROC,
FP=TP=rep(NA,101)plot(0:1,0:1,xlab="False Positive Rate",ylab="True Positive Rate",cex=.5)
for(s in seq(0,1,by=.01)){
Ps=(S>s)*1
FP[1+s*100]=sum((Ps==1)*(avocat==0))/sum(avocat==0)TP[1+s*100]=sum((Ps==1)*(avocat==1))/sum(avocat==1)}
lines(c(FP),c(TP),type="s",col="red")54

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Courbe ROC, receiver operating characteristic

55

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Courbe ROC, receiver operating characteristic
Mais il existe aussi des librairies de fonctions dédiées à cette analyse,
library(ROCR)
pred=prediction(S,avocat)
perf=performance(pred,"tpr", "fpr")
plot(perf,colorize = TRUE)

56

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Courbe ROC, receiver operating characteristic57

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression
On va utiliser la méthode CART (Classification And Regression Trees),
introduite par Breiman (1984). Pour aller plus loin, Berk (2008), chapitre 3.
Un arbre de décision est une règle de classification basé sur des tests associés aux
attributs organisés de manière arborescente.
Considérons - pour commencer - une variable continue X1 (ou X), et une
variable à expliquer Y qui prend les valeurs {0, 1}.

58Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression - indice de Gini
Étant donne une partition S = {S1 , · · · , Sk },

 X
k
X
|Si |
Si
|Si | Sjgini(S) =
· 1−
=
·
|S|
S
|S|
S
i=1
i6=j

Dans une régression de Y sur X (partitionné en 2 classes, notées A et B),
X nx X nx,y nx,y1−gini(Y |X) = −nnxnxx∈{A,B}y∈{0,1}

Dans une régression de Y sur X (partitionné en 2 classes, notées A et B),


X
X
nx
nx,y
nx,y
gini(Y |X) = −
1−
n
nx
nx
x∈{A,B,C}y∈{0,1}59

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression - indice de Gini
Le code est ici> u=sort(unique(B$age))> gini=(NA,length(u))> for(i in 2:(length(u))){+ I=B$age<(u[i]-.5)+ ft1=sum(B$Y[I==TRUE])/nrow(B)+ ff1=sum(B$Y[I==FALSE])/nrow(B)+ ft0=sum(1-B$Y[I==TRUE])/nrow(B)+ ff0=sum(1-B$Y[I==FALSE])/nrow(B)+ ft=ft0+ft1+ f0=ft0+ff0+ gini[i] = -((ft1+ft0)*(ft1/(ft1+ft0)*(1-ft1/(ft1+ft0))++ft0/(ft1+ft0)*(1-ft0/(ft1+ft0))) ++(ff1+ff0)*(ff1/(ff1+ff0)*(1-ff1/(ff1+ff0))++
ff0/(ff1+ff0)*(1-ff0/(ff1+ff0))))+ }60Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression - entropie

entropie(S) = −

k
X
|Si |

|S|
i=1


· log

Si
S



Dans une régression de Y sur X (partitionné en 2 classes, notées A et B),
X nx X nx,ynx,yentropie(Y |X) = −
log
n
nx
nx
x∈{A,B}y∈{0,1}

Dans une régression de Y sur X (partitionné en 2 classes, notées A, B et C),


X
X
nx
nx,y
nx,y
entropie(Y |X) = −
log
n
nx
nx
x∈{A,B,C}y∈{0,1}61

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression - entropie> u=sort(unique(B$age))> entropie=rep(NA,length(u))> for(i in 2:(length(u))){+ I=B$age<(u[i]-.5)+ ft1=sum(B$Y[I==TRUE])/nrow(B)+ ff1=sum(B$Y[I==FALSE])/nrow(B)+ ft0=sum(1-B$Y[I==TRUE])/nrow(B)+ ff0=sum(1-B$Y[I==FALSE])/nrow(B)+ ft=ft0+ft1+ f0=ft0+ff0+ entropie[i] = -((ft1+ft0)*(ft1/(ft1+ft0)*log(ft1/(ft1+ft0))++ft0/(ft1+ft0)*log(ft0/(ft1+ft0))) ++(ff1+ff0)*(ff1/(ff1+ff0)*log(ff1/(ff1+ff0))++
ff0/(ff1+ff0)*log(ff0/(ff1+ff0)))
+ )}

62

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression - critère du χ2
On peut aussi utiliser un test d’indépendance basée sur la distance du χ2 . Posons
n⊥
x,y =

nx · nynet la distance du χ2 quand X est partitionné en 2 classes, notées A et B),
χ2 (Y |X) =X

2
X [nx,y − n⊥
x,y ]
n⊥
x,y

x∈{A,B} y∈{0,1}alors que quand X (partitionné en 2 classes, notées A, B et C),
2

χ (Y |X) =X

2
X [nx,y − n⊥
x,y ]
n⊥
x,y

x∈{A,B,C} y∈{0,1}63Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression - critère du χ2> u=sort(unique(B$age))> gini=entropie=chideux=rep(NA,length(u))> for(i in 2:(length(u))){+ I=B$age<(u[i]-.5)+ ft1=sum(B$Y[I==TRUE])/nrow(B)+ ff1=sum(B$Y[I==FALSE])/nrow(B)+ ft0=sum(1-B$Y[I==TRUE])/nrow(B)+ ff0=sum(1-B$Y[I==FALSE])/nrow(B)+ M=matrix(c(ft0,ff0,ft1,ff1),2,2)+ ft=ft0+ft1+ f0=ft0+ff0+ Mind=matrix(c(ft*f0,f0*(1-ft),(1-f0)*ft,(1-f0)*(1-ft)),2,2)+ Q=sum(nrow(B)*(M-Mind)^2/Mind)+ chideux[i] = Q+ }64

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression (2 variables)
Considérons deux variables continues X1,i et X2,i .
> library(tree)> tree(Y~X1+X2)

65

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression (2 variables)
Une animation est présentée sur le blog.
Le principe est le même qu’avec une seule variable, avec un choix de variable (enplus du choix du seuil).
Pour trouver le premier noeud, on calcule l’indice de Gini pour tous les seuilspossibles, de X1 et de X2 .

66

Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression (2 variables)67Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression (2 variables)

68Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Arbre de régression (2 variables)
Une fois le premier noeud établi (ici sur X2 , en x2 ), le second est obtenu– soit en coupant sur X1 , pour les X2 < x2
– soit en coupant sur X1 , pour les X2 > x2
– soit en coupant sur X2 , pour les X2 < x2
– soit en coupant sur X2 , pour les X2 < x269Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Considérons une variable qualitative X1,i .
La variable qualitative prend des modaliés (non ordonnées) {x1,1 , · · · , x1,I }.
On peut alors calculer des indices de Gini (e.g.) pour toute partition {J, J C } de I.
70Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013

Couper les arbres
En pratique, sur les grosses bases de données, lire (et interpréter) des arbres peutêtre difficile,> sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",+ header=TRUE,sep=";")> sinistres=sinistre[sinistre$garantie=="1RC",]> 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)> base = merge(contrat,nombre)> Y = base$nbre>0> X1 = base$ageconducteur

71Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013> library(tree)> arbre=tree(Y~X1,split="gini")> plot(arbre)> text(arbre,cex=.8)

On peut demander à constituer des classes suffisamment larges
72Arthur CHARPENTIER - ACT2040 - Actuariat IARD - Automne 2013> arbre=tree(Y~X1,split="gini", mincut = 5000)> plot(arbre)> text(arbre,cex=.8)73