{"word_count": 1275, "character_count": 10074}
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 lancereg = glm(Good~ABV+IBU+Style,data=base,family=binomial)summary(reg)Coefficients:              Estimate Std. Error z value Pr(&#62;|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.9934Si on calcule la probabilité à partir de ce modèle, on obtientpredict(reg,newdata=newbeer,type="response")        1 0.2682531On 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(&#62;|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.7428et comme IBU n'est toujours pas significatif, on utilisereg = glm(Good~ABV+IBU,data=base,family=binomial)summary(reg)Coefficients:            Estimate Std. Error z value Pr(&#62;|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 suivantepredict(reg,newdata=newbeer,type="response")        1 0.2707103très proche de l'ancienne, de l'ordre de 27%. On peut juste regarder s'il n'y aurait pas d'effet non-linéairelibrary(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 populationmean(base$Good)[1] 0.2727273que l'on aurait obtenue avecreg = glm(Good~1,data=base,family=binomial)summary(reg)Coefficients:            Estimate Std. Error z value Pr(&#62;|z|)   (Intercept)  -0.9808     0.3385  -2.898  0.00376 **predict(reg,newdata=newbeer,type="response")        1 0.2727273Si 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.9679383Si on reprend le modèle avec trop de variables (non-significatives), on obtiendraitreg = glm(Good~ABV+IBU+Style,data=base,family=binomial)predict(reg,newdata=newbeer,type="response")        1 0.9639461dont on peut également demander la précision, en calculant l'écart-type de cette prévisionpredict(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éesLa seconde question était un peu plus technique, avec la base suivante, contenant des données agrégéesdownload.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 rejetsBA = 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       taux1  Female    A   89   19 0.824074072  Female    B   17    8 0.680000003  Female    C  202  391 0.340640814  Female    D  131  244 0.349333335  Female    E   94  299 0.239185756  Female    F   24  317 0.070381237    Male    A  512  313 0.620606068    Male    B  353  207 0.630357149    Male    C  120  205 0.3692307710   Male    D  138  279 0.3309352511   Male    E   53  138 0.2774869112   Male    F   22  351 0.05898123On retrouve, sur la première ligne, le 82.4% mentionné. Pour construite le modèle, on utilise la syntaxe suivantereggroup= glm(cbind(NAdm,NRej)~Gender+Dept,data=B,family=binomial)summary(reggroup)Coefficients:             Estimate Std. Error z value Pr(&#62;|z|)    (Intercept)   0.58205    0.06899   8.436   &#60;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   &#60;2e-16 ***DeptD        -1.29461    0.10582 -12.234   &#60;2e-16 ***DeptE        -1.73931    0.12611 -13.792   &#60;2e-16 ***DeptF        -3.30648    0.16998 -19.452   &#60;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")&#62; B   Gender Dept NAdm NRej       taux  predgroup1  Female    A   89   19 0.82407407 0.664167422  Female    B   17    8 0.68000000 0.654419633  Female    C  202  391 0.34064081 0.358776944  Female    D  131  244 0.34933333 0.351446965  Female    E   94  299 0.23918575 0.257809646  Female    F   24  317 0.07038123 0.067574507    Male    A  512  313 0.62060606 0.641539308    Male    B  353  207 0.63035714 0.631499129    Male    C  120  205 0.36923077 0.3361393110   Male    D  138  279 0.33093525 0.3290345111   Male    E   53  138 0.27748691 0.2391665412   Male    F   22  351 0.05898123 0.06154717Et 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 = NULLfor(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 baseregindiv = glm(Y~Gender+Dept,data=base,family=binomial)B$predindiv=predict(regindiv,type="response",newdata=B)B   Gender Dept NAdm NRej       taux  predgroup  predindiv1  Female    A   89   19 0.82407407 0.66416742 0.664167422  Female    B   17    8 0.68000000 0.65441963 0.654419633  Female    C  202  391 0.34064081 0.35877694 0.358776944  Female    D  131  244 0.34933333 0.35144696 0.351446965  Female    E   94  299 0.23918575 0.25780964 0.257809646  Female    F   24  317 0.07038123 0.06757450 0.067574507    Male    A  512  313 0.62060606 0.64153930 0.641539308    Male    B  353  207 0.63035714 0.63149912 0.631499129    Male    C  120  205 0.36923077 0.33613931 0.3361393110   Male    D  138  279 0.33093525 0.32903451 0.3290345111   Male    E   53  138 0.27748691 0.23916654 0.2391665412   Male    F   22  351 0.05898123 0.06154717 0.06154717Les 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(&#62;|z|)    (Intercept)   0.58205    0.06899   8.436   &#60;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   &#60;2e-16 ***DeptD        -1.29461    0.10582 -12.234   &#60;2e-16 ***DeptE        -1.73931    0.12611 -13.792   &#60;2e-16 ***DeptF        -3.30648    0.16998 -19.452   &#60;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 freedomResidual deviance: 5187.5  on 4519  degrees of freedomAIC: 5201.5Number of Fisher Scoring iterations: 5Bref, 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 &#38; O’Connell (1975), en lien avec le paradoxe de Simpson. Et je recommande largement la lecture !