Category Archives: STT5100

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 !

Slides 1 – loi binomiale

Allez, on reprend les cours, après un peu plus de deux semaines d’interruption… Mais je reprends tranquillement, en revenant sur des rappels que j’avais fait au premier cours, et dont on aura finalement besoin très bientôt ! On revient sur la loi de Bernoulli, la loi binomiale et la loi multinomiale. Et on parlera inférence, et tests. Les slides sont en ligne (slides 1) et la vidéo aussi (slides 1)

Il y a une coquille dans les slides (7 et 8) mais je n’ai pas interrompu la vidéo et j’essayais de continuer comme si de rien n’était. Pour la proportion \widehat{p}, on a juste la loi asymptotique Gaussienne (ou une transformation de la loi Binomiale à distance finie, mais pas de loi de Student, désolé pour le copier-coller depuis des vieux slides qui parlaient d’estimation de moyenne \overline{x} pour des observations supposées Gaussiennes)

J’essaye de limiter le nombre d’erreurs, mais ces cours à distance me prennent beaucoup de temps, et je n’ai pas le courage pour refaire les vidéos, pour l’instant. Toutes mes excuses.

 

 

Reprise des cours (à venir)

Lundi, les cours de la session d’hiver reprennent… à distance. Je vais mettre en ligne une série de capsules vidéos pour finir le cours, sur les GLM. J’ai mis en ligne une première vidéo (slides 0) pour annoncer le plan. J’ai fait des slides rapidement (ça changera des cours que je faisais au tableau) et je fais des enregistrement unique, sans montage, histoire de mettre en ligne rapidement le cours en ligne. Le pdf des slides est aussi en ligne (slides 0). J’essayerais de mettre les vidéos en ligne au fur et à mesure. Je ne suis pas particulièrement fier de moi, mais quitte à perdre du temps à faire le guignol devant la caméra, autant que ça serve au plus grand nombre.

Ah oui, il y a probablement des coquilles, voire des erreurs dans les slides… les commentaires sont ouverts pour me faire part de toute suggestion !