Après la loi de Bernoulli, la loi binomiale et la loi multinomiale, on attaque la régression logistique. Les slides sont en ligne (slides 2) et la vidéo aussi (slides 2)
Tag Archives: Logistique
Variables Catégorielles et Modèle Logistique
Petit complément, suite aux coquilles qu’il y avait dans les slides, sur une propriété de la régression logistique quand on régresse sur des variables catégorielles. On était sur la base des avocats
> avocat <- read.table("http://freakonometrics.free.fr/AutoBI.csv",header=TRUE,sep=",") > avocat$CLMSEX <- factor(avocat$CLMSEX, labels=c("M","F")) > avocat$MARITAL <- factor(avocat$MARITAL, labels=c("M","C","V","D")) > avocat=avocat[!is.na(avocat$CLMSEX),]
> attach(avocat)
> sum((ATTORNEY==2)&(CLMSEX=="F"),na.rm=TRUE)/sum(CLMSEX=="F",na.rm=TRUE)
[1] 0.5256065
> sum((ATTORNEY==2)&(CLMSEX=="M"),na.rm=TRUE)/sum(CLMSEX=="M",na.rm=TRUE)
[1] 0.4453925
Autrement dit, la proportion d’hommes qui se sont fait représentés par un avocat est de 44.539%, et la proportion de femmes 52.56%. On peut visualiser le tableau croisé ci-dessous
> tab=xtabs(~ATTORNEY+CLMSEX,data=avocat)
> require(vcd)
> mosaic(tab, shade=TRUE, legend=TRUE)
Quand on fait une régression linéaire (Gaussienne), on retrouve ces probabilités dans les valeurs de coefficients: pour les femmes, on a 44.539%, et pour les hommes, on a la différence avec les femmes (cette dernière modalité étant la modalité de référence ici)
> reglm = lm((ATTORNEY==2) ~ CLMSEX, data=avocat)
> summary(reglm)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.44539 0.02060 21.620 < 2e-16 ***
CLMSEXF 0.08021 0.02756 2.911 0.00367 **
C’est ce que l’on a ci-dessous
> sum(coefficients(reglm))
[1] 0.5256065
> coefficients(reglm)[1]
(Intercept)
0.4453925
On a un résultat similaire avec une régression logistique
> reglogit = glm((ATTORNEY==2) ~ CLMSEX, data=avocat,family=binomial)
> summary(reglogit)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.21930 0.08312 -2.639 0.00833 **
CLMSEXF 0.32182 0.11097 2.900 0.00373 **
même si les coefficients n’ont pas la même interprétation. En transformant ces coefficients, on retrouve très exactement les mêmes valeurs
> exp(sum(reglogit$coefficients[1:2])) /(1+exp(sum(reglogit$coefficients[1:2])))
[1] 0.5256065
> exp(reglogit$coefficients[1]) /(1+exp(reglogit$coefficients[1]))
(Intercept)
0.4453925
Désolé pour la typo.
Modèle de régression et interaction(s) entre facteurs
Dans un modèle de régression, on veut écrire
Quand on se limite à un modèle linéaire, on écrit
ou encore
Mais on de doute que l’on rate quelque chose… en particulier, on va rater toutes les interactions possibles. On peut croiser les variables, et supposer que
qui peut s’étendre d’avantage, à l’ordre 3,
voire davantage.
Supposons que nos variables soient ici qualitatives, et plus précisément binaires. Prenons un exemple simple, avec des données (classiques) en risque de crédit1. On peut trouver la base via
library(evtree) db=GermanCredit
ou encore directement
myVariableNames = c("checking_status","duration","credit_history",
"purpose","credit_amount","savings","employment","installment_rate",
"personal_status","other_parties","residence_since","property_magnitude",
"age","other_payment_plans","housing","existing_credits","job",
"num_dependents","telephone","foreign_worker","class")
GermanCredit = read.table(
"http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data",
header=FALSE,col.names=myVariableNames)
Retenons pour commencer trois variables explicatives,
db=data.frame(Y=GermanCredit$class-1,
X1=GermanCredit$checking_status%in%c("A12","A13"),
X2=GermanCredit$credit_history%in%c("A30","A31"),
X3=GermanCredit$savings%in%c("A61","A62"))
reg=glm(Y~X1+X2+X3,data=db,family=binomial)
summary(reg)
La régression sans interaction donne ici
Call:
glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5431 -0.8421 -0.6295 1.3994 1.9999
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.8544 0.1699 -10.915 < 2e-16 ***
X1TRUE 0.3363 0.1496 2.249 0.0245 *
X2TRUE 1.3462 0.2347 5.735 9.76e-09 ***
X3TRUE 1.0001 0.1787 5.596 2.19e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1221.7 on 999 degrees of freedom
Residual deviance: 1143.6 on 996 degrees of freedom
AIC: 1151.6
Number of Fisher Scoring iterations: 4
Il existe plusieurs interactions possibles ici (limitons nous aux paires). C’est ce que l’on observe quand on fait la régression
reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=db,family=binomial)
summary(reg)
Call:
glm(formula = Y ~ X1 + X2 + X3 + X1:X2 + X1:X3 + X2:X3, family = binomial,
data = db)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5369 -0.8281 -0.6439 1.3954 1.9638
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.77109 0.20070 -8.825 < 2e-16 ***
X1TRUE 0.30296 0.33737 0.898 0.369186
X2TRUE 0.88353 0.54255 1.628 0.103421
X3TRUE 0.87709 0.22583 3.884 0.000103 ***
X1TRUE:X2TRUE -0.37917 0.49343 -0.768 0.442225
X1TRUE:X3TRUE 0.09178 0.37278 0.246 0.805522
X2TRUE:X3TRUE 0.80923 0.58185 1.391 0.164293
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1221.7 on 999 degrees of freedom
Residual deviance: 1141.0 on 993 degrees of freedom
AIC: 1155
Number of Fisher Scoring iterations: 4
On peut faire un dessin pour visualiser les interactions : on a trois sommets (nos trois variables), et on visualiser les interactions
indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))
k=3
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",
xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)
ce qui donne ici, pour nos trois variables
Ce modèle pourrait sembler incomplet, car on ne regarde que les interactions entre les modalités, par paires. En fait, c’est parce qu’il manque (visuellement) les variables non-croisées. On peut les rajouter si on veut (au risque de surcharger le dessin)
cercle=function(c,r,cl) lines(c[1]+r*cos(seq(0,2*pi,length=501)),
c[2]+r*sin(seq(0,2*pi,length=501)),col=cl)
reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=db,family=binomial)
indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))
k=3
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)
soit ici
Si on change le ‘sens‘ de nos variables (en recodant a l’envers, en permutant les vrais et les faux), on obtient le graphique suivant
dbinv=db
dbinv[,2:k]=1-dbinv[,2:k]
reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=dbinv,family=binomial)
indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))
k=3
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)
qui peut alors être comparé au graphique précédant
Avec 5 variables, on augmente les interactions possibles… même si beaucoup risquent d’être non-significatifs. On peut déjà se focaliser sur les paires possibles d’interactions croisées. Pour simplifier le code, on va utiliser deux fonctions locales,
vrepeach=function(x,e){
v=NULL
for(i in 1:length(e)){v=c(v,rep(x[i],each=e[i]))}
return(v)}
vreplength=function(x,l){
v=NULL
for(i in 1:length(l)){v=c(v,x[l[i]:length(x)])}
return(v)}
et ensuite, on adapte le code précédant
indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)
ce qui donne un schéma plus complexe,
On peut aussi prendre juste 2 variables, prenant 3 et 4 modalités respectivement. On va extraire deux variables indicatrices pour la première (la modalité restante sera la modalité de référence) et trois pour la seconde,
db=data.frame(Y=GermanCredit$class-1,
X1=GermanCredit$checking_status=="A12",
X2=GermanCredit$checking_status=="A13",
X3=GermanCredit$checking_status=="A14",
X4=GermanCredit$employment%in%c("A72","A73"),
X5=GermanCredit$employment%in%c("A74","A75"))
k=5
indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
if(!is.na(coefficients(reg)[1+k+i])){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)
On voit que plusieurs interactions ne sont alors plus possibles, sur la partie gauche (les trois modalités de la même variable) et sur la partie droite
On peut d’ailleurs simplifier les graphs, en ne visualisant que les interactions significatives.
indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
if(!is.na(coefficients(reg)[1+k+i])){
if(summary(reg)$coefficients[1+k+i,4]<.1){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}}}
for(i in 1:k){
if(summary(reg)$coefficients[1+i]<.1){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)
soit ici
Ici, une seule interactions croisée est significative, et presque toutes les variables le sont. Et si on reprend le modèle avec 5 facteurs,
db=data.frame(Y=GermanCredit$class-1,X1=GermanCredit$checking_status%in%c("A12","A13"),
X2=GermanCredit$credit_history%in%c("A30","A31"),
X3=GermanCredit$savings%in%c("A61","A62"),
X4=GermanCredit$employment%in%c("A71","A72"),
X5=GermanCredit$other_payment_plans=="A143")
indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
if(!is.na(coefficients(reg)[1+k+i])){
if(summary(reg)$coefficients[1+k+i,4]<.1){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}}}
for(i in 1:k){
if(summary(reg)$coefficients[1+i]<.1){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)
on obtient
Je ne sais pas si mes graphiques sont pertinents, ou pas. Mais je trouve ça joli. En fait, je suis tombé un peu par hasard2 sur les Tables de Taguchi, développées par Gen’ichi Taguchi (田口 玄一). Le soucis est que je n’ai rien compris… Enfin, disons que je croyais comprendre, puis j’ai continué à faire des dessins… Si quelqu’un pourrait m’expliquer sur mon exemple les graphiques de Taguchi, je suis preneur ! car je doute que ce soit ce que je fais depuis tout à l’heure…
1. Cette base est largement utilisée dans le quatrième chapitre de Computational Actuarial Science with R, à paraître dans les mois à venir.
2.En l’occurence, le hasard est @Benavent qui a suscité ma curiosité ce matin en me parlant de ces tables, dont je n’avais alors jamais entendu parlé ! J’avais même lu rapidement Taniguchi (谷口 ジロー) et je ne voyais pas le rapport avec les statistiques….
Introduction à la régression logistique et aux arbres
Pour le second cours ACT2040, on va finir l’introduction (et les rappels de statistique inférentiel) puis attaque la première grosse section, sur la régression logistique et aux arbres de classification. base tirée du livre de Jed Frees, http://instruction.bus.wisc.edu/jfrees/…
> baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv",
+ header=TRUE,sep=",")
> tail(baseavocat)
CASENUM ATTORNEY CLMSEX MARITAL CLMINSUR SEATBELT CLMAGE LOSS
1335 34204 2 2 2 2 1 26 0.161
1336 34210 2 1 2 2 1 NA 0.576
1337 34220 1 2 1 2 1 46 3.705
1338 34223 2 2 1 2 1 39 0.099
1339 34245 1 2 2 1 1 18 3.277
1340 34253 2 2 2 2 1 30 0.688
On dispose d’une variable dichotomique indiquant si un assuré – suite à un accident de la route – a été représenté par un avocat (1 si oui, 2 si non). On connaît le sexe de l’assuré (1 pour les hommes et 2 pour les femmes), le statut marital (1 s’il est marié, 2 s’il est célibataire, 3 pour un veuf, et 4 pour un assuré divorcé). On sait aussi si l’assuré portait ou non une ceinture de sécurité lorsque l’accident s’est produit (1 si oui, 2 si non et 3 si l’information n’est pas connue). Enfin, une information pour savoir si le conducteur du véhicule était ou non assuré (1 si oui, 2 si non et 3 si l’information n’est pas connue). On va recoder un peu les données afin de les rendre plus claires à lire.
Les transparents sont en ligne sur le blog,
Des compléments théoriques sur les arbres peuvent se trouver en ligne http://genome.jouy.inra.fr/…, http://ensmp.fr/…, ou http://ujf-grenoble.fr/… (pour information, nous ne verrons que la méthode CART). Je peux renvoyer au livre (et au blog) de Stéphane Tuffery, ou (en anglais) au livre de Richard Berk, dont un résumé se trouve en ligne sur http://crim.upenn.edu/….
Régression logistique et arbres
Pour le cours de mercredi prochain, la base utilisée sera une base tirée du livre de Jed Frees, http://instruction.bus.wisc.edu/jfrees/…
> baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv",header=TRUE,sep=",") > tail(baseavocat) CASENUM ATTORNEY CLMSEX MARITAL CLMINSUR SEATBELT CLMAGE LOSS 1335 34204 2 2 2 2 1 26 0.161 1336 34210 2 1 2 2 1 NA 0.576 1337 34220 1 2 1 2 1 46 3.705 1338 34223 2 2 1 2 1 39 0.099 1339 34245 1 2 2 1 1 18 3.277 1340 34253 2 2 2 2 1 30 0.688
On dispose d’une variable dichotomique indiquant si un assuré – suite à un accident de la route – a été représenté par un avocat (1 si oui, 2 si non). On connaît le sexe de l’assuré (1 pour les hommes et 2 pour les femmes), le statut marital (1 s’il est marié, 2 s’il est célibataire, 3 pour un veuf, et 4 pour un assuré divorcé). On sait aussi si l’assuré portait ou non une ceinture de sécurité lorsque l’accident s’est produit (1 si oui, 2 si non et 3 si l’information n’est pas connue). Enfin, une information pour savoir si le conducteur du véhicule était ou non assuré (1 si oui, 2 si non et 3 si l’information n’est pas connue). On va recoder un peu les données afin de les rendre plus claires à lire.
Les transparents du cours sont en ligne ici,
Sur les arbres de régression, je mettrais en ligne un billet, afin d’illustrer la méthode. En attendant des compléments théoriques peuvent se trouver en ligne http://genome.jouy.inra.fr/…, http://ensmp.fr/…, ou http://ujf-grenoble.fr/… (pour information, nous ne verrons que la méthode CART). Je peux renvoyer au livre (et au blog) de Stéphane Tuffery, ou (en anglais) au livre de Richard Berk, dont un résumé se trouve en ligne sur http://crim.upenn.edu/….
La semaine suivante, nous aborderons la régression de Poisson, les méthodes de biais minimal, et introduirons les modèles linéaires généralisés. Je renvoie au chapitre sur la tarification a priori du Denuit & Charpentier (2005), aux chapitres 12 et 13 de Frees (2010) ou encore les chapitres 5 et 6 du De Jong & Heller (2008). Pour les plus curieux qui veulent comprendre les liens entre les modèles linéaires généralisés et la tarification par crédibilité, je renvoie à l’article de Klinker (2010)
Mais qui connaît la loi logistique ?
Certains élèves ont pris peur en voyant apparaitre la loi logisitque dans le cours d’économétrie des variables qualitatives, car autant j’avais passé pas mal de temps lors des remises à niveau sur la loi normale, la loi de Student, la loi de Poisson…. jamais je n’avais présenté la loi logistiquue.. Mais qui connaît la loi logisitique ?
- La loi logistique (en tant que loi)
En maths, on a introduit très tôt les fonctions de la forme
où K et r sont des réels positifs et a un réel quelconque, qui sont ce que l’on appelle les solutions (en temps continu) du modèle de Verhulst.
C’est Pierre François Verhulst qui a donné le nom de “courbe logistique“, en expliquant, en 1845 « nous donnerons le terme de logistique à cette courbe ». Comme toujours, il faut revenir aux racines étymologique, en notant que « logistique » a même racine que logarithme, et que logistikos signifie “calcul” en grec (on peut retrouver plus d’information ici , là ou encore là).
Ces fonctions ont été introduites par Pierre François Verhulst, qui était un élève d’Adolphe Quetelet, et qui dévait étudier un modèle d’évolution de population qui ne soit pas exponentielle. Cette courbe a été utilisée dans l’étude des populations à partir des années 20 par Raymond Pearl et Lowell Jacob Reed, mais qui ne créditent Verhulst de la paternité de la découverte qu’en 1922. En revanche, on ne retrouve le terme de “logistique” qu’en 1924 dans une correspondance entre George Yule et Reed. On notera que c’est Joseph Berkson qui défendra l’idée d’ajuster certaines courbes par une fonction logistique (modèle logit) plutôt que par la fonction de répartition de la loi de Gauss (modèle probit, introduit par Bliss en 1934, mais c’est une autre histoire que je raconterais sûrement un jour….).
En dynamique des populations, le modèle de Verhulst est un modèle de croissance proposé vers 1840. Pour reprendre l’explication de wikipedia sur le sujet (inspiré de plusieurs ouvrages sur le sujet dont de Handbook on logisitic distribution) Verhulst a proposé ce modèle en réponse au modèle de Malthus qui proposait un taux d’accroissement constant sans frein conduisant à une croissance exponentielle de la population. Le modèle de Verhulst imagine que le taux de natalité et le taux de mortalité sont des fonctions affines respectivement décroissante et croissante de la taille de la population. Autrement dit, plus la taille de la population augmente, plus son taux de natalité diminue et son taux de mortalité augmente.
On appelle y la taille de la population, m(y) le taux de mortalité et n(y) le taux de natalité. La taille de la population suit l’équation différentielle
Si m et n sont des fonction affines respectivement croissante et décroissante (ce qui une hypothèse un peu forte, mais qui a longtemps été utilisée comme on l’avait vu en cours d’actuariat en M1) alors n – m est une fonction affine décroissante. L’équation peut s’écrire
avec a et b deux réels positifs, si la croissance est positive quand y est proche de 0. En posant K=a/b, l’équation devient
avec a > 0 et K > 0 . Aussi, histoire de formaliser, on recherche des fonctions strictement positives définies sur 
et vérifiant le système 
et
ce qui conduit à la solution dite logistique
- La cote au lieu de la probabilité
Bref, voilà pour la loi…. mais comment voit-on apparaître cette loi en économétrie ? Pour rappel, on cherche à modéliser la probabilité de survenant d’un évènement. Mais il existe d’autre outils que la probabilité en tant que telle. En particulier, la cote d’un évènement – ou “odds in favor of an event” – est la quantité
où p est la probabilité de réalisation de cet évènement.
Cette grandeur est bien connue par les turfistes…. Bon, dans ce cadre (les paris), cette probabilité n’a rien d’une probabilité historique, mais il s’agit plutôt d’une probabilité risque neutre comme on dit en finance: il s’agit de la probabilité extraite des prix. Considérons un exemple simple (où on ne mise que sur le vainqueur)
- Le parieur A mise 10 € sur le cheval numéro 1.
- Le parieur B mise 5 € sur le cheval numéro 2.
- Le parieur C mise 20 € sur le cheval numéro 3.
- Le parieur D mise 15 € sur le cheval numéro 4.
On suppose que les sommes misées reflètent les croyances des joueurs (qui ont tous la même aversion pour le risques). Autrement dit, les joueurs pensent que le cheval 1 a deux fois plus de chances de gagner que le cheval 2. Bref, La cote d’un cheval est obtenue en divisant la somme des enjeux joués sur les autres chevaux par la somme des enjeux joués sur le cheval choisi (ou l’inverse, ça dépend de
la manière dont on souhaite le présenter), soit ici
- Pour le cheval 1, une cote de [5+20+15]/10=4 (on dira 4 contre 1)
- Pour
le cheval 2, une cote de [10+20+15]/5=9 (on dira 9 contre 1) - Pour le cheval 3, une cote de [10+5+15]/20=3/2 (on dira 1.5 contre 1)
- Pour le cheval 4, une cote de [10+5+20]/15=7/3
Ca existe également en football (comme en témoigne le ticket ci-contre où figurent les cotes). Bref, cette cote est une manière très naturelle de parler de probabilité. L’énorme avantage est qu’une probabilité est définie sur [0,1] alors qu’une cote peut prendre n’importe quelle valeur positive.
- La régression logistique
On souhaite donc modéliser la cote à l’aide de variables explicatives. Mais comme le score peut être négative, il est assez “naturel” de prendre le logarithme (comme pour la régression Poissonnienne par exemple). Aussi, on peut comprendre l’écriture de la forme suivante
ou encore
Or la fonction qui apparaît à gauche est la fonction quantile de la loi logistique. Quod erat demonstrandi…
- Un cas particulier de GLM
La transformation logistique apparaît très naturelle quand on utilise les GLM (modèles linéaires généralisés). En effet, il faut que la loi s’écrive sous la forme exponentielle
(comme j’en avais parlé ici ou là, par exemple). Pour reprendre les notations du cours d’actuariat 2, il faut écrire la loi binomiale
sous la forme^
soit ici (on n’a pas de paramètre de surdispersion)
En notant que la loi binomiale se réécrit
ou encore
On en déduit que le paramètre naturel est
et que
(moyennant un peu de réécriture). Bref, on note une relation naturelle (on parlera de fonction de lien canonique) de la forme entre le paramètre usuel de la loi binomiale et le paramètre dit naturel de la loi exponentielle,
ou encore
La fonction de répartition logistique apparaît là encore presque par enchantemant…. enfin, presque…