On continue sur les GLM avec la déviance et la présentation des résidus de déviance, et de Pearson. Les slides sont en ligne (slides 15) et la vidéo aussi (slides 15)
Category Archives: STT5100
Slides 14 – Lien et variance
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 !