logo

La régression logistique avec R, sa vie, son oeuvre

La régression logistique n’est pas, en soi, un modèle si compliqué. Mais les résultats qu’on peut en tirer sont nombreux, et la seule fonction glm du package {stats} ne les présente pas toujours de façon simple. Petit tour des possibilités offertes par des packages complémentaires pour tout comprendre des résultats d’un modèle.

Les données : passagers du Titanic

Le package {carData} contient un data.frame appelé TitanicSurvival. Il recense les 1309 passagers du navire avec leur âge, leur sexe et la classe à bord de laquelle ils voyageaient (1e, 2e ou 3e classe). Pour chaque passager, une indication est donnée de leur survie.

library(carData)   # données
library(tibble)    # visualisation des données

data("TitanicSurvival")
titanic <- TitanicSurvival # copie du data.frame sous un nom plus court

glimpse(titanic)
## Observations: 1,309
## Variables: 4
## $ survived       <fct> yes, yes, no, no, no, yes, yes, no, yes, no, no, yes...
## $ sex            <fct> female, male, female, male, female, male, female, ma...
## $ age            <dbl> 29.0000, 0.9167, 2.0000, 30.0000, 25.0000, 48.0000, ...
## $ passengerClass <fct> 1st, 1st, 1st, 1st, 1st, 1st, 1st, 1st, 1st, 1st, 1s...

Un petit peu de statistique descriptive nous informe sur le contenu de la table.

library(skimr)
pander(skim(titanic), split.tables=Inf)
skim_type skim_variable n_missing complete_rate factor.ordered factor.n_unique factor.top_counts numeric.mean numeric.sd numeric.p0 numeric.p25 numeric.p50 numeric.p75 numeric.p100
factor survived 0 1 FALSE 2 no: 809, yes: 500 NA NA NA NA NA NA NA
factor sex 0 1 FALSE 2 mal: 843, fem: 466 NA NA NA NA NA NA NA
factor passengerClass 0 1 FALSE 3 3rd: 709, 1st: 323, 2nd: 277 NA NA NA NA NA NA NA
numeric age 263 0.7991 NA NA NA 29.88 14.41 0.1667 21 28 39 80

On voit que les naufragés sont majoritaires (809/1309, soit environ 62%), que les femmes sont minoritaires (466 pour 1309 passagers) et que 263 passagers sont d’âge inconnu.

Que faire des données quantitatives ?

Une régression logistique peut utiliser une variable quantitative telle quelle, à condition : * de n’utiliser que les observations renseignées (les passagers sans âge ne seront pas pris en compte) * que la relation entre cette variable quantitative et le weight of evidence (ou woe, ou logit, ou log(probabilité(évènement)/(1-probabilité(évènement)))) soit linéaire.

Dans notre cas, abandonner 263 observations sur 1309 (soit 20% des passagers) serait une grosse perte d’information. On va donc découper la variable age en tranches, et créer une tranche spécifique pour les âges inconnus.

Commençons par regarder, pour les âges renseignés, la forme de la relation avec le weight of evidence.

library(dplyr)   # préparation des données
library(ggplot2) # graphiques
titanic %>% 
  filter(! is.na(age)) %>% 
  group_by(age) %>% 
  summarise(prob_survie = mean(survived == "yes"),
            n = n()) %>% 
  mutate(logit = log((prob_survie+0.5/n)/(1-prob_survie+0.5/n))) %>% 
  ggplot() +
  aes(x=age, y=logit, size=n) +
  geom_point() +
  geom_smooth(method="loess", se=FALSE, aes(weight=n), color="red") +
  geom_smooth(method="lm"   , se=FALSE, aes(weight=n), color="gray")
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

 

La relation entre les bulles (chacune représente un âge, la taille de la bulle étant proportionnelle au nombre de passagers de cet âge) ne semble pas spécialement linéaire. Si on opte pour un découpage en tranches, plutôt que d’ajuster le nuage des bulles par une droite oblique on l’ajustera par une courbe en escalier, dont chaque palier correspond à une tranche. On fait donc des tranches peu larges quand la relation évolue rapidement et des tranches plus larges quand les bulles sont à peu près toutes à la même hauteur dans le graphique.

Un premier découpage serait :

library(forcats)  # gestion des factors
titanic <- titanic %>% 
              mutate(age_cl = cut(age, breaks=c(-Inf,5,10,15,20,40,60,70,Inf)),
                     age_cl = fct_explicit_na(age_cl))

On obtient ce nouveau graphique.

titanic %>% 
  group_by(age_cl) %>% 
  summarise(prob_survie = mean(survived == "yes"),
            n = n()) %>% 
  mutate(logit = log((prob_survie+0.5)/(1-prob_survie+0.5))) %>% 
  ggplot() +
  aes(x=age_cl, y=logit) +
  geom_bar(stat="identity")

 

Il importe sur ce graphique que des tranches adjacentes ne correspondent pas à des barres de même hauteur. Cela semble correct ici.

Le nombre actuel de tranches peut sembler élevé : n’oublions pas qu’un modèle peut suggérer de fusionner des tranches a posteriori, mais pas d’éclater une tranche existante en plusieurs tranches plus petites.

Un premier modèle

La fonction glm permet de construire un premier modèle dans lequel on inclut toutes les variables potentiellement explicatives : l’âge (mis en tranches ci-dessus), le sexe ou la classe. Le premier argument de glm est, comme dans beaucoup de fonctions de modélisation de R, une formule de type Y ~ X1 + X2 où les variables Y, X1 et X2 doivent être dans le même data.frame

reg1 <- glm(survived ~ age_cl + sex + passengerClass, titanic, family = binomial)
pander(summary(reg1))
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.855 0.3965 9.721 2.444e-22
age_cl(5,10] -1.007 0.5584 -1.803 0.07139
age_cl(10,15] -0.7562 0.5742 -1.317 0.1879
age_cl(15,20] -1.544 0.403 -3.832 0.000127
age_cl(20,40] -1.477 0.3505 -4.213 2.518e-05
age_cl(40,60] -1.991 0.392 -5.08 3.784e-07
age_cl(60,70] -3.032 0.6642 -4.564 5.01e-06
age_cl(70, Inf] -2.027 1.11 -1.827 0.06776
age_cl(Missing) -1.676 0.3714 -4.514 6.372e-06
sexmale -2.5 0.1499 -16.68 1.847e-62
passengerClass2nd -1.152 0.2112 -5.458 4.827e-08
passengerClass3rd -2.036 0.1939 -10.5 8.204e-26

(Dispersion parameter for binomial family taken to be 1 )

Null deviance: 1741 on 1308 degrees of freedom
Residual deviance: 1218 on 1297 degrees of freedom

On peut récupérer les principaux résultats (coefficients et leur significativité) avec la fonction summary comme ci-dessus ou avec la fonction stargazer du package éponyme. Ce dernier sera particulièrement utile si on souhaite comparer plusieurs modèles : leurs caractéristiques s’afficheront côte à côte dans un même tableau.

library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
stargazer(reg1, type="html")
Dependent variable:
survived
age_cl(5,10] -1.007*
(0.558)
age_cl(10,15] -0.756
(0.574)
age_cl(15,20] -1.544***
(0.403)
age_cl(20,40] -1.477***
(0.350)
age_cl(40,60] -1.991***
(0.392)
age_cl(60,70] -3.032***
(0.664)
age_cl(70, Inf] -2.027*
(1.110)
age_cl(Missing) -1.676***
(0.371)
sexmale -2.500***
(0.150)
passengerClass2nd -1.152***
(0.211)
passengerClass3rd -2.036***
(0.194)
Constant 3.855***
(0.397)
Observations 1,309
Log Likelihood -608.936
Akaike Inf. Crit. 1,241.871
Note: p<0.1; p<0.05; p<0.01

Un autre package s’avère utile pour organiser les sorties du modèle au niveau de ses coefficients : il s’agit de {broom}.

library(broom)
kable(tidy(reg1))
term estimate std.error statistic p.value
(Intercept) 3.854830 0.3965304 9.721399 0.0000000
age_cl(5,10] -1.006734 0.5583757 -1.802969 0.0713931
age_cl(10,15] -0.756158 0.5742408 -1.316796 0.1879070
age_cl(15,20] -1.544280 0.4029821 -3.832130 0.0001270
age_cl(20,40] -1.476657 0.3504871 -4.213156 0.0000252
age_cl(40,60] -1.991385 0.3920417 -5.079525 0.0000004
age_cl(60,70] -3.031519 0.6641721 -4.564358 0.0000050
age_cl(70, Inf] -2.027379 1.1099279 -1.826586 0.0677621
age_cl(Missing) -1.676190 0.3713592 -4.513663 0.0000064
sexmale -2.500169 0.1498944 -16.679531 0.0000000
passengerClass2nd -1.152482 0.2111712 -5.457573 0.0000000
passengerClass3rd -2.036394 0.1938525 -10.504860 0.0000000

Dans tous les cas on remarque que pour chaque variable explicative qualitative (types texte et facteur) un coefficient n’est pas affiché, celui de la modalité de référence. Il est égal à 0 par construction. Par défaut R choisit la 1e valeur par ordre alphabétique comme référence ; la fonction relevel permet de changer ce choix.

titanic$age_cl <- relevel(titanic$age_cl, "(20,40]")
reg2 <- glm(survived ~ age_cl + sex + passengerClass, titanic, family = binomial)
kable(tidy(reg2))
term estimate std.error statistic p.value
(Intercept) 2.3781735 0.2080062 11.4331882 0.0000000
age_cl(-Inf,5] 1.4766569 0.3504871 4.2131555 0.0000252
age_cl(5,10] 0.4699227 0.4648991 1.0108057 0.3121094
age_cl(10,15] 0.7204989 0.4829013 1.4920209 0.1356937
age_cl(15,20] -0.0676231 0.2546823 -0.2655193 0.7906095
age_cl(40,60] -0.5147285 0.2204532 -2.3348655 0.0195504
age_cl(60,70] -1.5548623 0.5775239 -2.6922911 0.0070963
age_cl(70, Inf] -0.5507216 1.0611700 -0.5189758 0.6037776
age_cl(Missing) -0.1995336 0.2018029 -0.9887547 0.3227832
sexmale -2.5001688 0.1498944 -16.6795309 0.0000000
passengerClass2nd -1.1524823 0.2111712 -5.4575727 0.0000000
passengerClass3rd -2.0363937 0.1938525 -10.5048596 0.0000000

Avant de regarder en détail le sens de ces coefficients et des tests associés, intéressons-nous à une éventuelle simplification de ce modèle : peut-on faire aussi bien avec moins d’informations ?

Simplification du modèle

Simplification des variables

Un outil privilégié pour épurer un modèle est de faire des tests sur l’utilité de chaque variable explicative compte tenu de la présence des autres ; on sait ainsi quelle variable s’avère superflue, soit parce qu’elle ne fournit pas d’information pertinente, soit parce qu’elle s’avère redondante avec d’autres. On élimine ensuite la variable avec la p-value la plus élevée à un tel test, communément appelé test “de type III”. Une automatisation de ce processus conduit à une sélection Backward mais dans ce document nous nous limiterons à la stratégie manuelle. La fonction Anova du package {car} permet de réaliser ce test ; on repart de l’objet produit par la fonction glm.

library(car)
kable(Anova(reg2, type="III"))
LR Chisq Df Pr(>Chisq)
age_cl 39.35064 8 4.2e-06
sex 335.50113 1 0.0e+00
passengerClass 122.83504 2 0.0e+00

Ici toutes nos variables sont significatives, ce qui nous amène à toutes les conserver simultanément dans le modèle.

Simplification des modalités d’une variable

Une fois ce travail fait, on peut encore simplifier, cette fois à l’intérieur de chacune des variables explicatives : est-il nécessaire de conserver le détail de toutes les classes de passagers, de toutes les tranches d’âge dans le modèle ? Ou peut-on fusionner certaines de ces modalités en une valeur unique ? On va encore répondre à cette question par des tests, cette fois-ci avec l’aide du package {emmeans}. Celui-ci permet de combiner les coefficients d’une régression pour produire des prédictions “moyennes”. Dans ces constructions mathématiques, une seule des variables du modèle varie, les autres étant “moyennées” :

  • littéralement positionnées à leur valeur moyenne des données pour les variables X quantitatives du modèle
  • combinées comme si toutes leurs valeurs étaient représentées à parts égales dans les données pour les X qualitatives, par exemple sexe aura ses coefficients combinés avec chacun un poids de 0,5 comme si les données étaient composées pour moitié d’hommes et de femmes.

Le réalisme des valeurs produites par de telles constructions n’est pas toujours avéré, mais elles permettent au moins de produire des ordres de grandeur.

library(emmeans)
kable(emmeans(reg2, "passengerClass", type="response"))
passengerClass prob SE df asymp.LCL asymp.UCL
1st 0.7509343 0.0376548 Inf 0.6701813 0.8173055
2nd 0.4877822 0.0536589 Inf 0.3846611 0.5919538
3rd 0.2823585 0.0365134 Inf 0.2165359 0.3590210

Dans le tableau ci-dessus, on fait varier uniquement la classe des passagers. La colonne PROB indique un niveau de probabilité de survie prédit par le modèle : 75% chez les passagers de première classe, 49% chez ceux de deuxième classe et 28% chez ceux de troisième classe.

library(emmeans)
kable(prop.table(table(titanic$passengerClass,
                       titanic$survived),
                 1))
no yes
1st 0.3808050 0.6191950
2nd 0.5703971 0.4296029
3rd 0.7447109 0.2552891

On peut les comparer avec les pourcentages-lignes du tableau de fréquence obtenu sur les données brutes : 62% de survivants chez les passagers de 1e classe, 43% pour les 2e classe, 26% pour les 3e classe. Dans notre cas, les prédictions moyennes du modèle surévaluent les chances de survie des passagers. Le biais n’est pas systématiquement dans ce sens ; ici il est principalement dû aux disparités d’âge et de sexe au sein des classes, et au fait qu’il n’y a pas de parité hommes-femmes chez les passagers, ce qui rend irréalistes les projections de la fonction emmeans. Néanmoins on tire de cette fonction un moyen simple de comparer toutes les classes de passagers entre elles pour savoir s’il est utile au modèle de toutes les distinguer.

library(emmeans)
kable(pairs(emmeans(reg2, "passengerClass")))
contrast estimate SE df z.ratio p.value
1st – 2nd 1.1524823 0.2111712 Inf 5.457573 1.0e-07
1st – 3rd 2.0363937 0.1938525 Inf 10.504860 0.0e+00
2nd – 3rd 0.8839114 0.1881167 Inf 4.698738 7.8e-06

Toutes les comparaisons sont significatives : il n’y a pas lieu de simplifier. Faisons maintenant les comparaisons entre tranches d’âge.

library(emmeans)
kable(pairs(emmeans(reg2, "age_cl")))
contrast estimate SE df z.ratio p.value
(20,40] – (-Inf,5] -1.4766569 0.3504871 Inf -4.2131555 0.0008450
(20,40] – (5,10] -0.4699227 0.4648991 Inf -1.0108057 0.9849954
(20,40] – (10,15] -0.7204989 0.4829013 Inf -1.4920209 0.8597040
(20,40] – (15,20] 0.0676231 0.2546823 Inf 0.2655193 0.9999993
(20,40] – (40,60] 0.5147285 0.2204532 Inf 2.3348655 0.3209575
(20,40] – (60,70] 1.5548623 0.5775239 Inf 2.6922911 0.1503740
(20,40] – (70, Inf] 0.5507216 1.0611700 Inf 0.5189758 0.9998696
(20,40] – (Missing) 0.1995336 0.2018029 Inf 0.9887547 0.9869926
(-Inf,5] – (5,10] 1.0067341 0.5583757 Inf 1.8029690 0.6803412
(-Inf,5] – (10,15] 0.7561580 0.5742408 Inf 1.3167961 0.9266785
(-Inf,5] – (15,20] 1.5442799 0.4029821 Inf 3.8321302 0.0040416
(-Inf,5] – (40,60] 1.9913854 0.3920417 Inf 5.0795249 0.0000134
(-Inf,5] – (60,70] 3.0315192 0.6641721 Inf 4.5643580 0.0001733
(-Inf,5] – (70, Inf] 2.0273785 1.1099279 Inf 1.8265857 0.6644085
(-Inf,5] – (Missing) 1.6761905 0.3713592 Inf 4.5136629 0.0002196
(5,10] – (10,15] -0.2505761 0.6493946 Inf -0.3858611 0.9999866
(5,10] – (15,20] 0.5375458 0.5042753 Inf 1.0659768 0.9789559
(5,10] – (40,60] 0.9846512 0.4964923 Inf 1.9832153 0.5551651
(5,10] – (60,70] 2.0247850 0.7300606 Inf 2.7734481 0.1232186
(5,10] – (70, Inf] 1.0206443 1.1511905 Inf 0.8865990 0.9937014
(5,10] – (Missing) 0.6694563 0.4788847 Inf 1.3979487 0.8990258
(10,15] – (15,20] 0.7881219 0.5214683 Inf 1.5113515 0.8506649
(10,15] – (40,60] 1.2352274 0.5117643 Inf 2.4136646 0.2762289
(10,15] – (60,70] 2.2753612 0.7399818 Inf 3.0748881 0.0541304
(10,15] – (70, Inf] 1.2712205 1.1574324 Inf 1.0983108 0.9746458
(10,15] – (Missing) 0.9200325 0.4966329 Inf 1.8525401 0.6466726
(15,20] – (40,60] 0.4471055 0.3050664 Inf 1.4656003 0.8715346
(15,20] – (60,70] 1.4872393 0.6154593 Inf 2.4164705 0.2747069
(15,20] – (70, Inf] 0.4830985 1.0825849 Inf 0.4462454 0.9999588
(15,20] – (Missing) 0.1319105 0.2819856 Inf 0.4677918 0.9999408
(40,60] – (60,70] 1.0401338 0.5869071 Inf 1.7722289 0.7007376
(40,60] – (70, Inf] 0.0359931 1.0663269 Inf 0.0337543 1.0000000
(40,60] – (Missing) -0.3151949 0.2599726 Inf -1.2124162 0.9540635
(60,70] – (70, Inf] -1.0041407 1.1893316 Inf -0.8442900 0.9954916
(60,70] – (Missing) -1.3553287 0.5932384 Inf -2.2846275 0.3513664
(70, Inf] – (Missing) -0.3511880 1.0693059 Inf -0.3284261 0.9999962

Les comparaisons sont nombreuses et leur lecture, fastidieuse. Si on trie le tableau par p-value décroissante, on verra plus rapidement les fusions que le modèle encourage. Attention cependant, la fusion avec la p-value la plus élevée n’est pas forcément une bonne idée : elle peut correspondre concrètement à mélanger des profils incompatibles d’un point de vue métier. Ici, fusionner deux tranches d’âge non consécutives n’est pas souhaitable.

library(emmeans)
pairs(emmeans::emmeans(reg2,"age_cl")) %>% 
  as.data.frame() %>% 
  arrange(desc(p.value)) %>% 
  kable()
contrast estimate SE df z.ratio p.value
(40,60] – (70, Inf] 0.0359931 1.0663269 Inf 0.0337543 1.0000000
(20,40] – (15,20] 0.0676231 0.2546823 Inf 0.2655193 0.9999993
(70, Inf] – (Missing) -0.3511880 1.0693059 Inf -0.3284261 0.9999962
(5,10] – (10,15] -0.2505761 0.6493946 Inf -0.3858611 0.9999866
(15,20] – (70, Inf] 0.4830985 1.0825849 Inf 0.4462454 0.9999588
(15,20] – (Missing) 0.1319105 0.2819856 Inf 0.4677918 0.9999408
(20,40] – (70, Inf] 0.5507216 1.0611700 Inf 0.5189758 0.9998696
(60,70] – (70, Inf] -1.0041407 1.1893316 Inf -0.8442900 0.9954916
(5,10] – (70, Inf] 1.0206443 1.1511905 Inf 0.8865990 0.9937014
(20,40] – (Missing) 0.1995336 0.2018029 Inf 0.9887547 0.9869926
(20,40] – (5,10] -0.4699227 0.4648991 Inf -1.0108057 0.9849954
(5,10] – (15,20] 0.5375458 0.5042753 Inf 1.0659768 0.9789559
(10,15] – (70, Inf] 1.2712205 1.1574324 Inf 1.0983108 0.9746458
(40,60] – (Missing) -0.3151949 0.2599726 Inf -1.2124162 0.9540635
(-Inf,5] – (10,15] 0.7561580 0.5742408 Inf 1.3167961 0.9266785
(5,10] – (Missing) 0.6694563 0.4788847 Inf 1.3979487 0.8990258
(15,20] – (40,60] 0.4471055 0.3050664 Inf 1.4656003 0.8715346
(20,40] – (10,15] -0.7204989 0.4829013 Inf -1.4920209 0.8597040
(10,15] – (15,20] 0.7881219 0.5214683 Inf 1.5113515 0.8506649
(40,60] – (60,70] 1.0401338 0.5869071 Inf 1.7722289 0.7007376
(-Inf,5] – (5,10] 1.0067341 0.5583757 Inf 1.8029690 0.6803412
(-Inf,5] – (70, Inf] 2.0273785 1.1099279 Inf 1.8265857 0.6644085
(10,15] – (Missing) 0.9200325 0.4966329 Inf 1.8525401 0.6466726
(5,10] – (40,60] 0.9846512 0.4964923 Inf 1.9832153 0.5551651
(60,70] – (Missing) -1.3553287 0.5932384 Inf -2.2846275 0.3513664
(20,40] – (40,60] 0.5147285 0.2204532 Inf 2.3348655 0.3209575
(10,15] – (40,60] 1.2352274 0.5117643 Inf 2.4136646 0.2762289
(15,20] – (60,70] 1.4872393 0.6154593 Inf 2.4164705 0.2747069
(20,40] – (60,70] 1.5548623 0.5775239 Inf 2.6922911 0.1503740
(5,10] – (60,70] 2.0247850 0.7300606 Inf 2.7734481 0.1232186
(10,15] – (60,70] 2.2753612 0.7399818 Inf 3.0748881 0.0541304
(-Inf,5] – (15,20] 1.5442799 0.4029821 Inf 3.8321302 0.0040416
(20,40] – (-Inf,5] -1.4766569 0.3504871 Inf -4.2131555 0.0008450
(-Inf,5] – (Missing) 1.6761905 0.3713592 Inf 4.5136629 0.0002196
(-Inf,5] – (60,70] 3.0315192 0.6641721 Inf 4.5643580 0.0001733
(-Inf,5] – (40,60] 1.9913854 0.3920417 Inf 5.0795249 0.0000134

La première fusion concerne les tranches de 40 à 60 ans et les plus de 70 ans. Elle n’est pas à encourager car ces tranches ne sont pas adjacentes. La seconde fusion, en revanche, concerne les 15-20 ans et les 20-40 ans. On peut donc envisager une tranche unique de 15 à 40 ans.

titanic$age_cl <- fct_recode(titanic$age_cl,
                             "(15,40]"="(15,20]",
                             "(15,40]"="(20,40]")
reg2b <- glm(survived ~ age_cl + sex + passengerClass,
             titanic,
             family="binomial")
pairs(emmeans(reg2b,"age_cl")) %>% 
  as.data.frame() %>% 
  arrange(desc(p.value)) %>% 
  kable()
contrast estimate SE df z.ratio p.value
(40,60] – (70, Inf] 0.0367207 1.0664081 Inf 0.0344340 1.0000000
(70, Inf] – (Missing) -0.3549343 1.0693071 Inf -0.3319292 0.9999789
(5,10] – (10,15] -0.2502722 0.6493688 Inf -0.3854085 0.9999416
(15,40] – (70, Inf] 0.5408137 1.0606127 Inf 0.5099069 0.9996178
(60,70] – (70, Inf] -1.0038500 1.1893715 Inf -0.8440172 0.9905450
(5,10] – (70, Inf] 1.0247290 1.1511655 Inf 0.8901665 0.9870160
(15,40] – (Missing) 0.1858794 0.1951899 Inf 0.9523003 0.9807564
(15,40] – (5,10] -0.4839153 0.4618941 Inf -1.0476758 0.9670385
(10,15] – (70, Inf] 1.2750012 1.1574435 Inf 1.1015667 0.9566909
(40,60] – (Missing) -0.3182136 0.2597285 Inf -1.2251778 0.9246878
(-Inf,5] – (10,15] 0.7557670 0.5742275 Inf 1.3161456 0.8930759
(5,10] – (Missing) 0.6697947 0.4788389 Inf 1.3987891 0.8582081
(15,40] – (10,15] -0.7341875 0.4801801 Inf -1.5289837 0.7919186
(40,60] – (60,70] 1.0405708 0.5868308 Inf 1.7732041 0.6384021
(-Inf,5] – (5,10] 1.0060392 0.5583158 Inf 1.8019179 0.6187794
(-Inf,5] – (70, Inf] 2.0307682 1.1099472 Inf 1.8296079 0.5997048
(10,15] – (Missing) 0.9200669 0.4966361 Inf 1.8525977 0.5837882
(5,10] – (40,60] 0.9880082 0.4963088 Inf 1.9907127 0.4881226
(60,70] – (Missing) -1.3587843 0.5930383 Inf -2.2912254 0.2982793
(15,40] – (40,60] 0.5040930 0.2167912 Inf 2.3252460 0.2796266
(10,15] – (40,60] 1.2382804 0.5116609 Inf 2.4201194 0.2313830
(15,40] – (60,70] 1.5446637 0.5761783 Inf 2.6808778 0.1283360
(5,10] – (60,70] 2.0285790 0.7298583 Inf 2.7794147 0.1001215
(10,15] – (60,70] 2.2788512 0.7398365 Inf 3.0802091 0.0432004
(15,40] – (-Inf,5] -1.4899544 0.3469218 Inf -4.2947852 0.0004639
(-Inf,5] – (Missing) 1.6758338 0.3713210 Inf 4.5131676 0.0001722
(-Inf,5] – (60,70] 3.0346182 0.6640160 Inf 4.5700983 0.0001319
(-Inf,5] – (40,60] 1.9940474 0.3919127 Inf 5.0879888 0.0000100

Une fois la fusion réalisée avec le package {forcats}, on réitère les comparaisons de tranches, triées par p-value décroissante. La première fusion “acceptable” est la troisième, qui crée une classe unique de 5-15 ans à partir des 5-10 et 10-15.

titanic$age_cl <- fct_recode(titanic$age_cl,
                             "(5,15]"="(5,10]",
                             "(5,15]"="(10,15]")
reg2c <- glm(survived ~ age_cl + sex + passengerClass,
             titanic,
             family="binomial")
pairs(emmeans(reg2c,"age_cl")) %>% 
  as.data.frame() %>% 
  arrange(desc(p.value)) %>% 
  kable()
contrast estimate SE df z.ratio p.value
(40,60] – (70, Inf] 0.0369200 1.0665027 Inf 0.0346178 1.0000000
(70, Inf] – (Missing) -0.3557970 1.0693989 Inf -0.3327075 0.9998922
(15,40] – (70, Inf] 0.5415975 1.0607053 Inf 0.5106013 0.9987193
(60,70] – (70, Inf] -1.0039214 1.1894807 Inf -0.8439998 0.9802746
(15,40] – (Missing) 0.1858006 0.1952032 Inf 0.9518313 0.9639260
(5,15] – (70, Inf] 1.1455754 1.1073872 Inf 1.0344850 0.9461272
(40,60] – (Missing) -0.3188770 0.2597370 Inf -1.2276920 0.8835812
(40,60] – (60,70] 1.0408414 0.5868891 Inf 1.7734890 0.5661025
(15,40] – (5,15] -0.6039779 0.3400650 Inf -1.7760657 0.5643439
(-Inf,5] – (70, Inf] 2.0318876 1.1100395 Inf 1.8304641 0.5272403
(-Inf,5] – (5,15] 0.8863122 0.4630702 Inf 1.9139910 0.4708687
(5,15] – (Missing) 0.7897784 0.3628780 Inf 2.1764297 0.3083231
(60,70] – (Missing) -1.3597184 0.5930892 Inf -2.2926035 0.2474133
(15,40] – (40,60] 0.5046776 0.2168040 Inf 2.3278054 0.2305530
(15,40] – (60,70] 1.5455189 0.5762329 Inf 2.6821083 0.1027293
(5,15] – (40,60] 1.1086554 0.3843600 Inf 2.8844196 0.0599418
(5,15] – (60,70] 2.1494968 0.6586040 Inf 3.2637168 0.0189360
(15,40] – (-Inf,5] -1.4902900 0.3469426 Inf -4.2954943 0.0003498
(-Inf,5] – (Missing) 1.6760906 0.3713439 Inf 4.5135806 0.0001296
(-Inf,5] – (60,70] 3.0358090 0.6640695 Inf 4.5715227 0.0000988
(-Inf,5] – (40,60] 1.9949676 0.3919308 Inf 5.0901023 0.0000074

Le nouveau modèle encourage la fusion des 60-70 ans avec les 70 ans et plus. Voici les codes pour les fusions suivantes, qui ne seront pas commentées, elles sont toutes du même tonneau.

# 60-70 avec 70-Inf
titanic$age_cl <- fct_recode(titanic$age_cl,
                             "(60, Inf]"="(60,70]",
                             "(60, Inf]"="(70, Inf]")
reg2d <- glm(survived ~ age_cl + sex + passengerClass,
             titanic,
             family="binomial")
pairs(emmeans(reg2d,"age_cl")) %>% 
  as.data.frame() %>% 
  arrange(desc(p.value)) %>% 
  kable()
contrast estimate SE df z.ratio p.value
(15,40] – (Missing) 0.1861563 0.1951600 Inf 0.9538651 0.9322304
(40,60] – (Missing) -0.3189034 0.2597018 Inf -1.2279600 0.8232753
(40,60] – (60, Inf] 0.8453906 0.5240239 Inf 1.6132675 0.5896952
(15,40] – (5,15] -0.6039219 0.3399531 Inf -1.7764859 0.4810436
(-Inf,5] – (5,15] 0.8861312 0.4629356 Inf 1.9141563 0.3933051
(5,15] – (Missing) 0.7900782 0.3627665 Inf 2.1779248 0.2480386
(60, Inf] – (Missing) -1.1642941 0.5308100 Inf -2.1934292 0.2406756
(15,40] – (40,60] 0.5050598 0.2167652 Inf 2.3299854 0.1819492
(15,40] – (60, Inf] 1.3504504 0.5119808 Inf 2.6376972 0.0883530
(5,15] – (40,60] 1.1089816 0.3842558 Inf 2.8860503 0.0450935
(5,15] – (60, Inf] 1.9543723 0.6031152 Inf 3.2404628 0.0151378
(15,40] – (-Inf,5] -1.4900530 0.3468599 Inf -4.2958358 0.0002516
(-Inf,5] – (Missing) 1.6762094 0.3712539 Inf 4.5149949 0.0000926
(-Inf,5] – (60, Inf] 2.8405035 0.6087754 Inf 4.6659300 0.0000452
(-Inf,5] – (40,60] 1.9951128 0.3918493 Inf 5.0915314 0.0000053
# 40-60 avec 60-Inf
titanic$age_cl <- fct_recode(titanic$age_cl,
                             "(40, Inf]"="(40,60]",
                             "(40, Inf]"="(60, Inf]")
reg2e <- glm(survived ~ age_cl + sex + passengerClass,
             titanic,
             family="binomial")
pairs(emmeans(reg2e,"age_cl")) %>% 
  as.data.frame() %>% 
  arrange(desc(p.value)) %>% 
  kable()
contrast estimate SE df z.ratio p.value
(15,40] – (Missing) 0.1850310 0.1952045 Inf 0.9478829 0.8780711
(40, Inf] – (Missing) -0.4244489 0.2524640 Inf -1.6812258 0.4455949
(15,40] – (5,15] -0.6025311 0.3402284 Inf -1.7709609 0.3906355
(-Inf,5] – (5,15] 0.8859840 0.4632678 Inf 1.9124663 0.3105428
(5,15] – (Missing) 0.7875621 0.3630164 Inf 2.1694943 0.1913138
(15,40] – (40, Inf] 0.6094799 0.2080931 Inf 2.9288806 0.0280987
(5,15] – (40, Inf] 1.2120111 0.3798372 Inf 3.1908702 0.0123598
(15,40] – (-Inf,5] -1.4885151 0.3470500 Inf -4.2890514 0.0001744
(-Inf,5] – (Missing) 1.6735461 0.3714202 Inf 4.5058022 0.0000648
(-Inf,5] – (40, Inf] 2.0979950 0.3875480 Inf 5.4135098 0.0000006
# 5-15 avec 15-40
titanic$age_cl <- fct_recode(titanic$age_cl,
                             "(5,40]"="(5,15]",
                             "(5,40]"="(15,40]")
reg2f <- glm(survived ~ age_cl + sex + passengerClass,
             titanic,
             family="binomial")
pairs(emmeans(reg2f,"age_cl")) %>% 
  as.data.frame() %>% 
  arrange(desc(p.value)) %>% 
  kable()
contrast estimate SE df z.ratio p.value
(5,40] – (Missing) 0.2413447 0.1924441 Inf 1.254103 0.5922897
(40, Inf] – (Missing) -0.4077027 0.2522518 Inf -1.616253 0.3693494
(5,40] – (40, Inf] 0.6490474 0.2069428 Inf 3.136362 0.0092727
(5,40] – (-Inf,5] -1.4347233 0.3456244 Inf -4.151106 0.0001938
(-Inf,5] – (Missing) 1.6760680 0.3715880 Inf 4.510554 0.0000383
(-Inf,5] – (40, Inf] 2.0837707 0.3875032 Inf 5.377428 0.0000005

On en arrive à un point où plus rien n’est simplifiable, et où il nous reste 4 tranches : âge inconnu, moins de 5 ans, 5 à 40 ans et plus de 40 ans.

Commentaire du modèle

Utilisons maintenant un autre package, {blorr}, pour avoir les sorties détaillées de ce modèle. Elles sont organisées suivant la présentation de la proc LOGISTIC de SAS et permet aux utilisateurs de ce logiciel de retrouver leurs repères. (La mise en forme HTML pour intégration à une page Web comme dans {rmarkdown} n’est pour le moment pas disponible, ni avec {knitr} ni avec {pander}.)

library(blorr)
blr_regress(reg2f, odd_conf_limit=TRUE)
## - Creating model overview. 
## - Creating response profile. 
## - Extracting maximum likelihood estimates. 
## - Computing odds ratio estimates. 
## - Estimating concordant and discordant pairs.
##                              Model Overview                              
## ------------------------------------------------------------------------
## Data Set    Resp Var    Obs.    Df. Model    Df. Residual    Convergence 
## ------------------------------------------------------------------------
##   data      survived    1309      1308           1302           TRUE     
## ------------------------------------------------------------------------
## 
##                     Response Summary                     
## --------------------------------------------------------
## Outcome        Frequency        Outcome        Frequency 
## --------------------------------------------------------
##    0              809              1              500    
## --------------------------------------------------------
## 
##                       Maximum Likelihood Estimates                       
## ------------------------------------------------------------------------
##     Parameter        DF    Estimate    Std. Error    z value     Pr(>|z|) 
## ------------------------------------------------------------------------
##    (Intercept)       1      2.4014        0.2027     11.8463      0.0000 
##  age_cl(-Inf,5]      1      1.4347        0.3456      4.1511      0.0000 
##  age_cl(40, Inf]     1     -0.6490        0.2069     -3.1364      0.0017 
##  age_cl(Missing)     1     -0.2413        0.1924     -1.2541      0.2098 
##      sexmale         1     -2.5127        0.1493    -16.8278      0.0000 
## passengerClass2nd    1     -1.1353        0.2097     -5.4135      0.0000 
## passengerClass3rd    1     -2.0031        0.1917    -10.4486      0.0000 
## ------------------------------------------------------------------------
## 
##                         Odds Ratio Estimates                         
## --------------------------------------------------------------------
##      Effects              Estimate             95% Wald Conf. Limit 
## --------------------------------------------------------------------
##  age_cl(-Inf,5]            4.1985            2.1514           8.3734 
##  age_cl(40, Inf]           0.5225            0.3467           0.7810 
##  age_cl(Missing)           0.7856            0.5370           1.1428 
##      sexmale               0.0811            0.0602           0.1081 
## passengerClass2nd          0.3213            0.2121           0.4830 
## passengerClass3rd          0.1349            0.0922           0.1955 
## --------------------------------------------------------------------
## 
##  Association of Predicted Probabilities and Observed Responses  
## ---------------------------------------------------------------
## % Concordant          0.8049          Somers' D        0.7121   
## % Discordant          0.1353          Gamma            0.6696   
## % Tied                0.0597          Tau-a            0.3164   
## Pairs                 404500          c                0.8348   
## ---------------------------------------------------------------

On retrouve dans ces sorties les coefficients de la régression, avec la même présentation que dans la fonction tidy de {broom}.

kable(tidy(reg2f))
term estimate std.error statistic p.value
(Intercept) 2.4013701 0.2027105 11.846301 0.0000000
age_cl(-Inf,5] 1.4347233 0.3456244 4.151106 0.0000331
age_cl(40, Inf] -0.6490474 0.2069428 -3.136362 0.0017106
age_cl(Missing) -0.2413447 0.1924441 -1.254103 0.2098046
sexmale -2.5126650 0.1493161 -16.827829 0.0000000
passengerClass2nd -1.1353278 0.2097209 -5.413517 0.0000001
passengerClass3rd -2.0030860 0.1917088 -10.448589 0.0000000

Deux tableaux sont inédits par rapport aux sorties précédentes : les odds-ratios et le tableau de paires (“Association of Predicted Probabilities etc.”).

Odds-ratios

Les odds-ratios (OR) sont les exponentielles des coefficients de régression, qu’on peut interpréter approximativement comme un facteur multiplicatif des chances de survenance de l’évènement. Si on prend l’OR pour SEX=male, on se compare à la référence SEX=female. Cet OR vaut 0,08 – soit environ 1/12. On peut donc dire que les chances de survie d’un homme sont celles d’une femme multipliée par 0,08, autrement dit, divisées par 12. Pour le tourner encore autrement, les chances de survie d’un homme sont 12 fois plus faibles que celles d’une femme. Si l’OR est supérieur à 1, les chances de survenance de l’évènement sont en faveur de la catégorie considérée par rapport à la référence. Sur l’âge, on compare les enfants de moins de 5 ans à la référence des personnes de 5 à 40 ans. Les petits enfants ont 4,2 fois plus de chances de survie que leurs aînés. Ces OR sont assortis d’intervalles de confience ; quand ils incluent la valeur 1, l’OR n’est pas significatif donc pas commentable. Cela correspond à un coefficient dont la p-value est supérieure à 0,05.

On peut aussi les obtenir avec la fonction tidy de {broom} en ajoutant l’option exponentiate=TRUE.

# on change l'âge de référence pour 5-40, la classe pour 2e et le sexe, homme
titanic$age_cl         <- relevel(titanic$age_cl          , "(5,40]")
titanic$sex            <- relevel(titanic$sex             , "male")
titanic$passengerClass <- relevel(titanic$passengerClass  , "2nd")
reg2g <- glm(survived ~ age_cl + sex + passengerClass,
             titanic,
             family="binomial")
broom::tidy(reg2g, exponentiate=TRUE, conf.int=TRUE)
## # A tibble: 7 x 7
##   term              estimate std.error statistic  p.value conf.low conf.high
##   <chr>                <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)          0.287     0.164     -7.59 3.28e-14    0.207     0.395
## 2 age_cl(-Inf,5]       4.20      0.346      4.15 3.31e- 5    2.15      8.37 
## 3 age_cl(40, Inf]      0.523     0.207     -3.14 1.71e- 3    0.347     0.781
## 4 age_cl(Missing)      0.786     0.192     -1.25 2.10e- 1    0.537     1.14 
## 5 sexfemale           12.3       0.149     16.8  1.53e-63    9.25     16.6  
## 6 passengerClass1st    3.11      0.210      5.41 6.18e- 8    2.07      4.71 
## 7 passengerClass3rd    0.420     0.186     -4.67 3.00e- 6    0.291     0.604

Avec le résultat de tidy, il est plus simple de produire un graphique.

library(ggplot2)
reg2g %>% 
  tidy(exponentiate=TRUE, conf.int=TRUE) %>% 
  ggplot() +
    aes(x=estimate, y=term, xmin=conf.low, xmax=conf.high) +
    geom_errorbarh(height=0.2) +
    geom_point(shape=15, size=3)

 

Tableau des paires

L’autre résultat inédit de {blorr} est le tableau des paires. Pour rappel le voici.

## - Creating model overview. 
## - Creating response profile. 
## - Extracting maximum likelihood estimates. 
## - Estimating concordant and discordant pairs.
pairs concordance discordance tied somers_d gamma tau c
404500 0.8049172 0.135335 0.0597478 0.7121305 0.6695822 0.3163771 0.8347911

Il est construit en constituant depuis les données toutes les paires possibles, chacune composée d’un passager survivant et d’un passager naufragé. Les prédictions du modèle sont calculées sur chaque passager ; on notera pS la probabilité de survie du passager réellement survivant et pN celle du passager naufragé. On compte

  1. comme concordante une paire où pS > pN, ce qui est la logique espérée : qu’un passager ayant réellement survécu soit aux yeux du modèle celui qui avait le plus de chances de survivre.
  2. comme discordante une paire où pS < pN.
  3. comme liée (en anglais tied) une paire où pS=pN.

Le pourcentage de paires concordantes est une bonne approximation (légèrement optimiste) du pourcentage d’individus bien classés en prédiction par le modèle. Mais la combinaison du nombre de paires des différents types produit également d’autres statistiques : D de Somers (également appelé “indice de Gini”), c de Pregibon (ou AUC), ou encore Tau-alpha de Kendall. Plus ces indicateurs sont proches de 1, meilleure est la qualité discriminante du modèle (= sa capacité à reconnaître un survivant comme tel, et un naufragé comme tel).

Prédictions

Cela nous amène à calculer, justement, des prédictions du modèle pour chaque passager. La fonction predict va calculer pour chacun la probabilité de l’évènement, ici la survie.

titanic$p_surv <- predict(reg2g,
                          titanic,
                          type="response")
head(titanic$p_surv)
## [1] 0.9169317 0.7897515 0.9788780 0.4722050 0.9169317 0.3185720

De ce nombre entre 0 et 1, on peut faire plusieurs usages…

  • soit le transformer en décision binaire (“survie”, “naufrage”) en utilisant un seuil ; il n’y a pas de valeur standard pour ce seuil, il doit être déterminé au mieux pour chaque modèle
  • soit l’utiliser comme un critère de classement, s’il s’agit d’identifier les 100 passagers avec les meilleures chances de survie (un classement d’un cynisme achevé, dont on préfère ne pas réfléchir à l’usage). C’est l’usage qu’on fait d’un modèle pour réaliser un ciblage : clients auxquels on envoie une publicité, clients à contacter en priorité pour éviter un départ à la concurrence, entreprises à aider pour éviter un redressement judiciaire, contribuables les plus suspects de fraude… Dans ce cas, pas de seuil, juste un classement décroissant selon les prédictions brutes

Modèle pour classer

Pour classer et obtenir une prédiction binaire, il nous faut donc choisir un seuil. Une possibilité est d’utiliser le taux de survivants (le % d’évènements) dans la base initiale.

mean(titanic$survived == "yes")
## [1] 0.381971

On peut en déduire une matrice de confusion.

confusion <- table(titanic$survived,
                   ifelse(titanic$p_surv < 0.382,
                          "prédit  \"non\"",
                          "prédit  \"oui\""))
kable(confusion)
prédit “non” prédit “oui”
no 622 187
yes 112 388
# taux de bien classés
sum(diag(confusion))/sum(confusion)*100
## [1] 77.15814

En représentant graphiquement la répartition des vrais survivants et vrais naufragés selon la prédiction du modèle, on peut aussi trouver un seuil optimal.

ggplot(titanic) +
  aes(x=p_surv, group=survived, colour=survived) +
  geom_line(stat="density", lwd=1.2) +
  theme_minimal() +
  scale_x_continuous(breaks=seq(0,1,by=0.1))

 

On voit un croisement des courbes vers… 0,38 – exactement le seuil proposé plus haut. C’est une propriété des régressions logistiques qui se vérifiera systématiquement. On peut essayer un seuil à 0,42 pour voir s’il y a une grosse différence dans les taux de bien classés.

confusion <- table(titanic$survived,
      ifelse(titanic$p_surv < 0.42,
             "prédit  \"non\"",
             "prédit  \"oui\""))
# kable(confusion)
confusion
##      
##       prédit  "non" prédit  "oui"
##   no            642           167
##   yes           120           380
# taux de bien classés
sum(diag(confusion))/sum(confusion)*100
## [1] 78.07487

On peut aussi utiliser une autre visualisation pour définir un seuil “optimal” (au sens où le modèle se tromper de manière équitable dans les deux sens : autant de faux positifs que de faux négatifs) : un diagramme de distributions cumulées. C’est sur ce type de comparaisons qu’est basé le test de Kolmogorov-Smirnov. Ce graphique est facilement produit par le package {blorr}.

reg2g %>%
    blr_gains_table() %>%
    blr_ks_chart()

 

Dernier graphique pour les recherches de seuils, la courbe ROC où chaque point correspond à un seuil possible. Le point de la courbe le plus proche du coin supérieur gauche – le point (0,1) – correspond au seuil optimal. De nombreux packages dessinent des courbes ROC, mais pour une régression logistique on peut rester avec {blorr}.

reg2g %>%
    blr_gains_table() %>%
    blr_roc_curve()

 

L’aire sous la courbe ROC est l’indicateur AUC (ou c dans les sorties de blr_regress) vu plus haut.

Modèle pour cibler

Pour l’autre utilisation prédictive d’un modèle, le ciblage, où la prédiction ne sert que de critère de tri sans que ses valeurs dans l’absolu soient exploitées, le graphique pertinent est la courbe de lift. Elle sera aussi fournie par {blorr}.

reg2g %>%
  blr_gains_table() %>% 
  plot()

 

On y voit en abscisses la population exprimée en pourcentage, triée par prédiction décroissante : les 10% des passagers avec le plus fort pourcentage de chances de survie, les 20% des passagers avec le plus fort pourcentage de chances de survie, etc. (c’est cumulatif), jusqu’à l’intégralité des passagers à 100%. L’axe des ordonnées indique le pourcentage de vrais survivants dans cette fraction de la population.

N’importe quel modèle va du point (0,0) au point (100,100) car en ne sélectionnant personne on ne sélectionne aucun individu réellement cible ; et en sélectionnant tout le monde on sélectionne forcément toutes les cibles. Entre ces deux points, beaucoup de trajectoires sont possibles, y compris des courbes qui démarrent sous la diagonale avant de repasser au-dessus ; ce type de modèle est à fuir absolument. Un modèle parfait peut également être facilement dessiné sur cette courbe, il est composé de deux segments de droite. En notant p la proportion de cibles dans la population, le premier segment est oblique et va de (0,0) à (p,100) ; le second segment est horizontal et va de (p,100) à (100,100).

parfait <- titanic %>% 
  summarise(p=mean(survived == "yes"))
reg2g %>%
  blr_gains_table() %>% 
  plot() +
  geom_segment(data=parfait,
               aes(x=0, y=0,
                   xend=p, yend=1),
               linetype="dashed") +
  geom_segment(data=parfait,
               aes(x=p, y=1,
                   xend=1, yend=1),
               linetype="dashed")

 

26 found this helpful