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"))

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"))

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))

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")))

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")))

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()

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()

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()

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()

# 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()

# 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()

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.

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.

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).

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)

# 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")