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 !