Tag Archives: quiz

STT5100, quiz (régression de Poisson #1)

Pour la fin de la semaine, j’avais donné un petit quiz, basé sur la base suivante, qui donne le nombre de cyclistes à une intersection, à New York City,

download.file("http://freakonometrics.free.fr/NYCVelo.RData","velo.RData")
load("velo.RData")
str(base)
'data.frame':	214 obs. of  7 variables:
 $ Date    : chr  "1-Apr-17" "2-Apr-17" "3-Apr-17" "4-Apr-17" ...
 $ HIGH_T  : num  46 62.1 63 51.1 63 48.9 48 55.9 66 73.9 ...
 $ LOW_T   : num  37 41 50 46 46 41 43 39.9 45 55 ...
 $ PRECIP  : num  0 0 0.03 1.18 0 0.73 0.01 0 0 0 ...
 $ BB_COUNT: int  606 2021 2470 723 2807 461 1222 1674 2375 3324 ...
 $ DAY     : chr  "Sam" "Dim" "Lun" "Mar" ...
 $ DIFF_T  : num  9 21.1 13 5.1 17 7.9 5 16 21 18.9 ......

En utiisant une régression de Poisson, il fallait prédire combien de cyclistes passeront un dimanche, s’il fait une température maximale de 85F, minimale de 70F, et s’il ne pleut pas. Et il fallait voir ce que donnerait la prévision pour un lundi.

newbase = data.frame(DAY=as.factor(c("Lun","Dim")),
 HIGH_T=c(85,85),LOW_T=c(70,70),
 PRECIP=c(0,0))

Faisons un modèle avec toutes les variables explicatives.

reg0 = glm(BB_COUNT~HIGH_T+LOW_T+PRECIP+DAY,data=base,family=poisson)

Rajoutons aussi une indicatrice pour indiquer les jours où il ne pleut pas,

reg = glm(BB_COUNT~HIGH_T+LOW_T+PRECIP+I(PRECIP==0)+DAY,data=base,family=poisson)
summary(reg) 
 
Coefficients:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)         6.8844970  0.0110463 623.241   <2e-16 ***
HIGH_T              0.0210950  0.0003133  67.328   <2e-16 ***
LOW_T              -0.0114006  0.0003351 -34.024   <2e-16 ***
PRECIP             -0.6570450  0.0071899 -91.384   <2e-16 ***
I(PRECIP == 0)TRUE  0.1303908  0.0033283  39.176   <2e-16 ***
DAYJeu              0.1683475  0.0049690  33.880   <2e-16 ***
DAYLun              0.1144129  0.0050480  22.665   <2e-16 ***
DAYMar              0.1655886  0.0049936  33.160   <2e-16 ***
DAYMer              0.1688035  0.0049190  34.317   <2e-16 ***
DAYSam              0.0466447  0.0051838   8.998   <2e-16 ***
DAYVen              0.1003536  0.0050919  19.709   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for poisson family taken to be 1)
 
    Null deviance: 70021  on 213  degrees of freedom
Residual deviance: 26493  on 203  degrees of freedom
AIC: 28580
 
Number of Fisher Scoring iterations: 4

Tout semble vraiment significatif. Et même avoir du sens (si on regarde les signes, par exemple). Et si on a peur de rater un effet non-linéaire, on peut mettre des splines sur toutes les variables continues

library(gam)
reggam = gam(BB_COUNT~bs(HIGH_T)+bs(LOW_T)+bs(PRECIP)+I(PRECIP==0)+DAY,data=base,family=poisson)
plot(reggam, se=TRUE)

pour la température maximale, ou pour la température minimale

et la courbe suivante pour la précipitation, avec un lissage linéaire entre l’observation maximale (3) et celle juste avant (vers 1.8)

On peut aussi régresser sur la température minimale, et l’écart entre le maximum et le minimum (dans un modèle linéaire les modèles sont équivalents, mais avec de transformation nonlinéaire, la différence pourrait donner un modèle plus simple)

library(gam)
reggam2 = gam(BB_COUNT~bs(DIFF_T)+bs(LOW_T)+bs(PRECIP)+I(PRECIP==0)+DAY,data=base,family=poisson)
plot(reggam2, se=TRUE)

On peut maintenant comparer les quatre modèles, et leurs prévision. Par exemple, pour le modèle linéaire (avec la variable indicatrice indiquant qu’il ne pleut pas),

P = predict(reg,newdata=newbase,type="response",se.fit=TRUE)

pour le lundi, on obtient l’intervalle de confiance à 95% pour \widehat{\lambda}

P$fit[1]+c(-2,2)*P$se.fit[1]
[1] 3349.842 3401.395

et pour le dimanche l’intervalle de confiance à 95% est

P$fit[2]+c(-2,2)*P$se.fit[2]
[1] 2987.497 3033.861

On peut visualiser les quatre intervalles de confiance (avec, de haut en bas, le second modèle gam, le premier, le modèle linéaire avec l’indicatrice d’absence de pluie, puis le premier modèle linéaire)

alors que pour le dimanche, on a

Autrement dit, en changeant de modèle, on change beaucoup les intervalles de confiance sur la prévision (et les intervalles sont parfois complètement disjoints). Ce qui n’est pas forcément une bonne nouvelle.

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 !