STT5100, quiz (régression logistique)

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 ...

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")

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

Si on calcule la probabilité à partir de ce modèle, on obtient

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

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 **

Cette fois, la variable restante est significative. On obtient la prévision suivante

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)

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

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

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

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

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
  • 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 ...

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

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 ***

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

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 ...

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

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

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 !


Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.