Je mets en ligne des éléments de correction pour le quiz donné lundi dernier, en lien avec la régression logistique.
- Quiz 1 – régression logistique (sur données individuelles)
On disposait de la base suivante, avec 44 bières testées, et un jugement porté (bon, ou pas),
download.file("http://freakonometrics.free.fr/Beer.RData","Beer.RData")
load("Beer.RData")
str(base)
'data.frame': 44 obs. of 7 variables:
$ Brewery : Factor w/ 8 levels "Bauhaus","Bent Paddle",..: 1 1 1 1 2 2 2 2 2 4 ...
$ Beer : Factor w/ 44 levels "14* ESB","B-Side Pils",..: 44 38 43 36 19 42 5 1 7 11 ...
$ Description: Factor w/ 37 levels "American Blonde Ale",..: 22 15 37 20 32 28 16 13 3 26 ...
$ Style : Factor w/ 3 levels "Ale","IPA","Lager": 3 3 3 2 1 3 2 1 1 1 ...
$ ABV : num 5.4 5 5.4 6.7 4.8 5 6.2 5.6 6 5.4 ...
$ IBU : int 48 28 55 70 48 38 68 32 34 45 ...
$ Good : int 0 0 0 0 0 0 0 0 0 1 ... |
download.file("http://freakonometrics.free.fr/Beer.RData","Beer.RData")
load("Beer.RData")
str(base)
'data.frame': 44 obs. of 7 variables:
$ Brewery : Factor w/ 8 levels "Bauhaus","Bent Paddle",..: 1 1 1 1 2 2 2 2 2 4 ...
$ Beer : Factor w/ 44 levels "14* ESB","B-Side Pils",..: 44 38 43 36 19 42 5 1 7 11 ...
$ Description: Factor w/ 37 levels "American Blonde Ale",..: 22 15 37 20 32 28 16 13 3 26 ...
$ Style : Factor w/ 3 levels "Ale","IPA","Lager": 3 3 3 2 1 3 2 1 1 1 ...
$ ABV : num 5.4 5 5.4 6.7 4.8 5 6.2 5.6 6 5.4 ...
$ IBU : int 48 28 55 70 48 38 68 32 34 45 ...
$ Good : int 0 0 0 0 0 0 0 0 0 1 ...
Le but est d’estimer la probabilité de trouver une bière bonne, à partir de trois caractéristiques: le degré d’alcool (ABV-Alcool By Volume), la mesure d’amertume (IBU-International Bitterness Unit) et le style (trois styles sont ici présent, les lager, les ale, et les IPA – qui techniquement sont des ale, mais peu importe),
newbeer = data.frame(ABV=6,IBU=78,Style="IPA") |
newbeer = data.frame(ABV=6,IBU=78,Style="IPA")
Pour répondre, on se lance
reg = glm(Good~ABV+IBU+Style,data=base,family=binomial)
summary(reg)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -10.98183 4.71577 -2.329 0.0199 *
ABV 1.42984 0.79147 1.807 0.0708 .
IBU 0.04905 0.04005 1.225 0.2207
StyleIPA -2.42680 1.67669 -1.447 0.1478
StyleLager -16.69972 2009.50187 -0.008 0.9934 |
reg = glm(Good~ABV+IBU+Style,data=base,family=binomial)
summary(reg)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -10.98183 4.71577 -2.329 0.0199 *
ABV 1.42984 0.79147 1.807 0.0708 .
IBU 0.04905 0.04005 1.225 0.2207
StyleIPA -2.42680 1.67669 -1.447 0.1478
StyleLager -16.69972 2009.50187 -0.008 0.9934
Si on calcule la probabilité à partir de ce modèle, on obtient
predict(reg,newdata=newbeer,type="response")
1
0.2682531 |
predict(reg,newdata=newbeer,type="response")
1
0.2682531
On peut tenter de simplifier un peu le modèle, en notant que le style semblant peu significatif,
reg = glm(Good~ABV+IBU,data=base,family=binomial)
summary(reg)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -9.307699 3.678245 -2.530 0.0114 *
ABV 1.308989 0.721717 1.814 0.0697 .
IBU 0.008248 0.025139 0.328 0.7428 |
reg = glm(Good~ABV+IBU,data=base,family=binomial)
summary(reg)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -9.307699 3.678245 -2.530 0.0114 *
ABV 1.308989 0.721717 1.814 0.0697 .
IBU 0.008248 0.025139 0.328 0.7428
et comme IBU n’est toujours pas significatif, on utilise
reg = glm(Good~ABV+IBU,data=base,family=binomial)
summary(reg)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -9.7881 3.3959 -2.882 0.00395 **
ABV 1.4662 0.5481 2.675 0.00747 ** |
reg = glm(Good~ABV+IBU,data=base,family=binomial)
summary(reg)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -9.7881 3.3959 -2.882 0.00395 **
ABV 1.4662 0.5481 2.675 0.00747 **
Cette fois, la variable restante est significative. On obtient la prévision suivante
predict(reg,newdata=newbeer,type="response")
1
0.2707103 |
predict(reg,newdata=newbeer,type="response")
1
0.2707103
très proche de l’ancienne, de l’ordre de 27%. On peut juste regarder s’il n’y aurait pas d’effet non-linéaire
library(gam)
reg = gam(Good~bs(ABV),data=base,family=binomial)
plot(reg,se=TRUE) |
library(gam)
reg = gam(Good~bs(ABV),data=base,family=binomial)
plot(reg,se=TRUE)
et non, on peut garder l’idée d’un modèle linéaire.
Par contre, on est un peu frustrés parce que la prévision qu’on ferait correspond à la fréquence de bières jugées bonnes dans la population
mean(base$Good)
[1] 0.2727273 |
mean(base$Good)
[1] 0.2727273
que l’on aurait obtenue avec
reg = glm(Good~1,data=base,family=binomial)
summary(reg)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.9808 0.3385 -2.898 0.00376 **
predict(reg,newdata=newbeer,type="response")
1
0.2727273 |
reg = glm(Good~1,data=base,family=binomial)
summary(reg)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.9808 0.3385 -2.898 0.00376 **
predict(reg,newdata=newbeer,type="response")
1
0.2727273
Si le taux d’alcool passe de 6% à 9%, cette probabilité passe de 27% à 97%
reg = glm(Good~ABV,data=base,family=binomial)
newbeer = data.frame(ABV=9,IBU=78,Style="IPA")
predict(reg,newdata=newbeer,type="response")
1
0.9679383 |
reg = glm(Good~ABV,data=base,family=binomial)
newbeer = data.frame(ABV=9,IBU=78,Style="IPA")
predict(reg,newdata=newbeer,type="response")
1
0.9679383
Si on reprend le modèle avec trop de variables (non-significatives), on obtiendrait
reg = glm(Good~ABV+IBU+Style,data=base,family=binomial)
predict(reg,newdata=newbeer,type="response")
1
0.9639461 |
reg = glm(Good~ABV+IBU+Style,data=base,family=binomial)
predict(reg,newdata=newbeer,type="response")
1
0.9639461
dont on peut également demander la précision, en calculant l’écart-type de cette prévision
predict(reg,newdata=newbeer,type="response",se.fit=TRUE)
$fit
1
0.9639461
$se.fit
1
0.0680774 |
predict(reg,newdata=newbeer,type="response",se.fit=TRUE)
$fit
1
0.9639461
$se.fit
1
0.0680774
- Quiz 2 : régression logistique sur données agrégées
La seconde question était un peu plus technique, avec la base suivante, contenant des données agrégées
download.file("http://freakonometrics.free.fr/BerkeleyAdmissions.RData","BerkeleyAdmissions.RData")
load("BerkeleyAdmissions.RData")
str(database)
'data.frame': 24 obs. of 4 variables:
$ Admit : Factor w/ 2 levels "Admitted","Rejected": 1 2 1 2 1 2 1 2 1 2 ...
$ Gender: Factor w/ 2 levels "Male","Female": 1 1 2 2 1 1 2 2 1 1 ...
$ Dept : Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 2 2 2 2 3 3 ...
$ Freq : num 512 313 89 19 353 207 17 8 120 205 ... |
download.file("http://freakonometrics.free.fr/BerkeleyAdmissions.RData","BerkeleyAdmissions.RData")
load("BerkeleyAdmissions.RData")
str(database)
'data.frame': 24 obs. of 4 variables:
$ Admit : Factor w/ 2 levels "Admitted","Rejected": 1 2 1 2 1 2 1 2 1 2 ...
$ Gender: Factor w/ 2 levels "Male","Female": 1 1 2 2 1 1 2 2 1 1 ...
$ Dept : Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 2 2 2 2 3 3 ...
$ Freq : num 512 313 89 19 353 207 17 8 120 205 ...
Par exemple, pour une personne de sexe féminin, ayant postulé dans le département A, il y a eu 89 admissions, sur 108 dossiers reçus, soit 82.4% de chances d’être admise. La question était d’estimer la probabilité d’admission avec une régression logistique, sur les facteurs Gender et Dept (sans effet croisé) était utilisée ?
Pour commencer, j’avoue que ce format de données ne me plait pas trop : par genre et par département, on va rajouter deux variables : le nombre d’acceptation et le nombre de rejets
BA = database[database$Admit=="Admitted",]
BA=BA[,2:4]
names(BA)=c("Gender","Dept","NAdm")
BR = database[database$Admit=="Rejected",]
BR=BR[,2:4]
names(BR)=c("Gender","Dept","NRej")
B=merge(BA,BR)
B$taux = B$NAdm/(B$NAdm+B$NRej)
B
Gender Dept NAdm NRej taux
1 Female A 89 19 0.82407407
2 Female B 17 8 0.68000000
3 Female C 202 391 0.34064081
4 Female D 131 244 0.34933333
5 Female E 94 299 0.23918575
6 Female F 24 317 0.07038123
7 Male A 512 313 0.62060606
8 Male B 353 207 0.63035714
9 Male C 120 205 0.36923077
10 Male D 138 279 0.33093525
11 Male E 53 138 0.27748691
12 Male F 22 351 0.05898123 |
BA = database[database$Admit=="Admitted",]
BA=BA[,2:4]
names(BA)=c("Gender","Dept","NAdm")
BR = database[database$Admit=="Rejected",]
BR=BR[,2:4]
names(BR)=c("Gender","Dept","NRej")
B=merge(BA,BR)
B$taux = B$NAdm/(B$NAdm+B$NRej)
B
Gender Dept NAdm NRej taux
1 Female A 89 19 0.82407407
2 Female B 17 8 0.68000000
3 Female C 202 391 0.34064081
4 Female D 131 244 0.34933333
5 Female E 94 299 0.23918575
6 Female F 24 317 0.07038123
7 Male A 512 313 0.62060606
8 Male B 353 207 0.63035714
9 Male C 120 205 0.36923077
10 Male D 138 279 0.33093525
11 Male E 53 138 0.27748691
12 Male F 22 351 0.05898123
On retrouve, sur la première ligne, le 82.4% mentionné. Pour construite le modèle, on utilise la syntaxe suivante
reggroup= glm(cbind(NAdm,NRej)~Gender+Dept,data=B,family=binomial)
summary(reggroup)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.58205 0.06899 8.436 <2e-16 ***
GenderFemale 0.09987 0.08085 1.235 0.217
DeptB -0.04340 0.10984 -0.395 0.693
DeptC -1.26260 0.10663 -11.841 <2e-16 ***
DeptD -1.29461 0.10582 -12.234 <2e-16 ***
DeptE -1.73931 0.12611 -13.792 <2e-16 ***
DeptF -3.30648 0.16998 -19.452 <2e-16 *** |
reggroup= glm(cbind(NAdm,NRej)~Gender+Dept,data=B,family=binomial)
summary(reggroup)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.58205 0.06899 8.436 <2e-16 ***
GenderFemale 0.09987 0.08085 1.235 0.217
DeptB -0.04340 0.10984 -0.395 0.693
DeptC -1.26260 0.10663 -11.841 <2e-16 ***
DeptD -1.29461 0.10582 -12.234 <2e-16 ***
DeptE -1.73931 0.12611 -13.792 <2e-16 ***
DeptF -3.30648 0.16998 -19.452 <2e-16 ***
On peut alors faire notre prévision de la probabilité d’être acceptée, pour une femme postulant dans le département A. On peut rajouter les probabilités prédite dans notre base,
B$predgroup=predict(reggroup,type="response")
> B
Gender Dept NAdm NRej taux predgroup
1 Female A 89 19 0.82407407 0.66416742
2 Female B 17 8 0.68000000 0.65441963
3 Female C 202 391 0.34064081 0.35877694
4 Female D 131 244 0.34933333 0.35144696
5 Female E 94 299 0.23918575 0.25780964
6 Female F 24 317 0.07038123 0.06757450
7 Male A 512 313 0.62060606 0.64153930
8 Male B 353 207 0.63035714 0.63149912
9 Male C 120 205 0.36923077 0.33613931
10 Male D 138 279 0.33093525 0.32903451
11 Male E 53 138 0.27748691 0.23916654
12 Male F 22 351 0.05898123 0.06154717 |
B$predgroup=predict(reggroup,type="response")
> B
Gender Dept NAdm NRej taux predgroup
1 Female A 89 19 0.82407407 0.66416742
2 Female B 17 8 0.68000000 0.65441963
3 Female C 202 391 0.34064081 0.35877694
4 Female D 131 244 0.34933333 0.35144696
5 Female E 94 299 0.23918575 0.25780964
6 Female F 24 317 0.07038123 0.06757450
7 Male A 512 313 0.62060606 0.64153930
8 Male B 353 207 0.63035714 0.63149912
9 Male C 120 205 0.36923077 0.33613931
10 Male D 138 279 0.33093525 0.32903451
11 Male E 53 138 0.27748691 0.23916654
12 Male F 22 351 0.05898123 0.06154717
Et on obtient une probabilité d’être acceptatée de 66%. Mais notons qu’on peut procéder autrement (on reviendra sur cette idée plus tard dans le cours): on peut recréer la base originale, qui a permis d’obtenir le tableau de contingence.
base = NULL
for(i in 1:nrow(B)){
L = B[i,1:4]
vL = data.frame(Gender=L[1],Dept=L[2],Y=c(rep(1,L[3]),rep(0,L[4])))
base = rbind(base,vL)}
str(base)
'data.frame': 4526 obs. of 3 variables:
$ Gender: Factor w/ 2 levels "Male","Female": 2 2 2 2 2 2 2 2 2 2 ...
$ Dept : Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Y : num 1 1 1 1 1 1 1 1 1 1 ... |
base = NULL
for(i in 1:nrow(B)){
L = B[i,1:4]
vL = data.frame(Gender=L[1],Dept=L[2],Y=c(rep(1,L[3]),rep(0,L[4])))
base = rbind(base,vL)}
str(base)
'data.frame': 4526 obs. of 3 variables:
$ Gender: Factor w/ 2 levels "Male","Female": 2 2 2 2 2 2 2 2 2 2 ...
$ Dept : Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Y : num 1 1 1 1 1 1 1 1 1 1 ...
On peut alors faire la régression sur cette base, avec les 4526 individus, et rajouter les prévisions, dans la base
regindiv = glm(Y~Gender+Dept,data=base,family=binomial)
B$predindiv=predict(regindiv,type="response",newdata=B)
B
Gender Dept NAdm NRej taux predgroup predindiv
1 Female A 89 19 0.82407407 0.66416742 0.66416742
2 Female B 17 8 0.68000000 0.65441963 0.65441963
3 Female C 202 391 0.34064081 0.35877694 0.35877694
4 Female D 131 244 0.34933333 0.35144696 0.35144696
5 Female E 94 299 0.23918575 0.25780964 0.25780964
6 Female F 24 317 0.07038123 0.06757450 0.06757450
7 Male A 512 313 0.62060606 0.64153930 0.64153930
8 Male B 353 207 0.63035714 0.63149912 0.63149912
9 Male C 120 205 0.36923077 0.33613931 0.33613931
10 Male D 138 279 0.33093525 0.32903451 0.32903451
11 Male E 53 138 0.27748691 0.23916654 0.23916654
12 Male F 22 351 0.05898123 0.06154717 0.06154717 |
regindiv = glm(Y~Gender+Dept,data=base,family=binomial)
B$predindiv=predict(regindiv,type="response",newdata=B)
B
Gender Dept NAdm NRej taux predgroup predindiv
1 Female A 89 19 0.82407407 0.66416742 0.66416742
2 Female B 17 8 0.68000000 0.65441963 0.65441963
3 Female C 202 391 0.34064081 0.35877694 0.35877694
4 Female D 131 244 0.34933333 0.35144696 0.35144696
5 Female E 94 299 0.23918575 0.25780964 0.25780964
6 Female F 24 317 0.07038123 0.06757450 0.06757450
7 Male A 512 313 0.62060606 0.64153930 0.64153930
8 Male B 353 207 0.63035714 0.63149912 0.63149912
9 Male C 120 205 0.36923077 0.33613931 0.33613931
10 Male D 138 279 0.33093525 0.32903451 0.32903451
11 Male E 53 138 0.27748691 0.23916654 0.23916654
12 Male F 22 351 0.05898123 0.06154717 0.06154717
Les prévisions sont identiques, et les modèles aussi,
summary(regindiv)
Call:
glm(formula = Y ~ Gender + Dept, family = binomial, data = base)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4773 -0.9306 -0.3741 0.9588 2.3613
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.58205 0.06899 8.436 <2e-16 ***
GenderFemale 0.09987 0.08085 1.235 0.217
DeptB -0.04340 0.10984 -0.395 0.693
DeptC -1.26260 0.10663 -11.841 <2e-16 ***
DeptD -1.29461 0.10582 -12.234 <2e-16 ***
DeptE -1.73931 0.12611 -13.792 <2e-16 ***
DeptF -3.30648 0.16998 -19.452 <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: 6044.3 on 4525 degrees of freedom
Residual deviance: 5187.5 on 4519 degrees of freedom
AIC: 5201.5
Number of Fisher Scoring iterations: 5 |
summary(regindiv)
Call:
glm(formula = Y ~ Gender + Dept, family = binomial, data = base)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4773 -0.9306 -0.3741 0.9588 2.3613
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.58205 0.06899 8.436 <2e-16 ***
GenderFemale 0.09987 0.08085 1.235 0.217
DeptB -0.04340 0.10984 -0.395 0.693
DeptC -1.26260 0.10663 -11.841 <2e-16 ***
DeptD -1.29461 0.10582 -12.234 <2e-16 ***
DeptE -1.73931 0.12611 -13.792 <2e-16 ***
DeptF -3.30648 0.16998 -19.452 <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: 6044.3 on 4525 degrees of freedom
Residual deviance: 5187.5 on 4519 degrees of freedom
AIC: 5201.5
Number of Fisher Scoring iterations: 5
Bref, la probabilité prédite par le modèle linéaire, sans effet croisé, serait de l’ordre de 66%. Petit moment culturel: cette dernière base est tirée d’un article relativement connu, Bickel, Hammel & O’Connell (1975), en lien avec le paradoxe de Simpson. Et je recommande largement la lecture !
