Tutoriel du jour 2

Régression logistique avec les données Titanic

Objectif

Ce tutoriel construit et évalue une régression logistique de la survie enregistrée dans l’échantillon titanic_train. Il suit la même progression que la présentation et les trois missions:

1. comprendre la réponse binaire et la portée des données;
2. passer de la probabilité à la cote et à la log-cote;
3. ajuster le modèle et interpréter les rapports de cotes;
4. produire des probabilités sur un jeu test;
5. distinguer confusion, discrimination et calibration;
6. choisir un seuil selon un mandat explicite;
7. vérifier les conditions de validité et les limites.

Les calculs restent accessibles. Les concepts nécessaires à une utilisation défendable du modèle sont néanmoins présentés au complet.

Carte théorique

La régression logistique relie quatre objets qu’il faut garder distincts.

Objet	Question
Réponse binaire	Quel événement est observé pour chaque unité?
Probabilité conditionnelle	Quelle probabilité est associée à un profil donné?
Rapport de cotes	Comment les cotes changent-elles lorsque le profil change?
Règle de seuil	À partir de quelle probabilité attribue-t-on une classe?

Le modèle produit d’abord une probabilité. La classe prédite n’apparaît qu’après le choix d’un seuil. Changer le seuil ne réajuste pas le modèle: cela modifie la décision construite à partir des mêmes probabilités.

Préparer l’environnement

library(tidyverse)
library(tidymodels)
library(titanic)
library(janitor)
library(broom)
library(performance)
library(scales)
library(patchwork)

theme_set(theme_minimal(base_size = 14))
tidymodels_prefer()

Le paquet titanic fournit des observations individuelles issues du jeu utilisé pour la compétition Kaggle. Le tableau titanic_train contient l’indicateur de survie et plusieurs caractéristiques enregistrées pour 891 passagers (P. Hendricks, 2015).

1. Charger et documenter les données

data("titanic_train")

titanic_eiom <- titanic_train |>
  clean_names() |>
  transmute(
    passenger_id,
    survie = survived,
    verite = factor(
      if_else(survived == 1, "survecu", "non_survecu"),
      levels = c("non_survecu", "survecu")
    ),
    classe = factor(
      pclass,
      levels = c(3, 2, 1),
      labels = c("3e", "2e", "1re")
    ),
    sexe = factor(
      sex,
      levels = c("male", "female"),
      labels = c("homme", "femme")
    ),
    age,
    tarif = fare,
    famille_bord = sib_sp + parch,
    embarquement = na_if(embarked, "")
  )

glimpse(titanic_eiom)

Rows: 891
Columns: 9
$ passenger_id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
$ survie       <int> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, …
$ verite       <fct> non_survecu, survecu, survecu, survecu, non_survecu, non_…
$ classe       <fct> 3e, 1re, 3e, 1re, 3e, 3e, 1re, 3e, 3e, 2e, 3e, 1re, 3e, 3…
$ sexe         <fct> homme, femme, femme, femme, homme, homme, homme, homme, f…
$ age          <dbl> 22, 38, 26, 35, 35, NA, 54, 2, 27, 14, 4, 58, 20, 39, 14,…
$ tarif        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.8625…
$ famille_bord <int> 1, 1, 0, 1, 0, 0, 0, 4, 2, 1, 2, 0, 0, 6, 0, 0, 5, 0, 1, …
$ embarquement <chr> "S", "C", "S", "S", "S", "Q", "S", "S", "S", "C", "S", "S…

Une ligne représente un passager de cet échantillon. La population visée par une conclusion doit donc rester proche de ce fichier. Il ne s’agit pas d’un manifeste complet ni d’une base permettant de reconstituer toutes les circonstances du naufrage.

La variable numérique survie sert à l’ajustement avec glm(). Le facteur verite sert aux mesures de classification. Dans les deux cas, l’événement positif est la survie enregistrée.

2. Comprendre la réponse binaire

titanic_eiom |>
  count(verite) |>
  mutate(proportion = n / sum(n))

       verite   n proportion
1 non_survecu 549  0.6161616
2     survecu 342  0.3838384

Dans cet échantillon, 342 passagers sur 891 sont enregistrés comme survivants, soit environ 38,4 %. Une règle qui prédit toujours non_survecu aurait donc une exactitude d’environ 61,6 %. Cette référence rappelle qu’une exactitude doit être comparée à une règle simple.

Une réponse binaire vaut 0 ou 1 pour chaque passager. Une probabilité, elle, est une quantité comprise entre 0 et 1 qui décrit l’incertitude conditionnelle avant d’observer la réponse.

Conceptuellement, le modèle suppose:

\[ Y_i \mid X_i \sim \operatorname{Bernoulli}(p_i), \]

où \(Y_i\) est la réponse du passager \(i\), \(X_i\) son profil et \(p_i\) sa probabilité conditionnelle de survie enregistrée.

3. Explorer sans confondre association et cause

survie_sexe <- titanic_eiom |>
  count(sexe, verite) |>
  group_by(sexe) |>
  mutate(proportion = n / sum(n)) |>
  ggplot(aes(sexe, proportion, fill = verite)) +
  geom_col(position = "fill") +
  scale_y_continuous(labels = label_percent()) +
  labs(
    x = NULL,
    y = "Proportion",
    fill = "Survie enregistrée"
  )

survie_classe <- titanic_eiom |>
  count(classe, verite) |>
  group_by(classe) |>
  mutate(proportion = n / sum(n)) |>
  ggplot(aes(classe, proportion, fill = verite)) +
  geom_col(position = "fill") +
  scale_y_continuous(labels = label_percent()) +
  labs(
    x = "Classe",
    y = "Proportion",
    fill = "Survie enregistrée"
  )

survie_sexe + survie_classe

Les proportions observées diffèrent fortement selon le sexe et la classe enregistrés. Ces graphiques montrent des associations marginales. Ils ne contrôlent pas l’âge, le tarif ou d’autres caractéristiques et ne démontrent pas un mécanisme causal.

Une formulation prudente est:

Dans cet échantillon, la survie enregistrée est plus fréquente dans certains groupes.

Une formulation non soutenue serait:

Changer artificiellement une caractéristique aurait produit la différence observée.

4. Examiner les valeurs manquantes

titanic_eiom |>
  summarise(
    lignes = n(),
    age_manquant = sum(is.na(age)),
    tarif_manquant = sum(is.na(tarif)),
    embarquement_manquant = sum(is.na(embarquement))
  )

  lignes age_manquant tarif_manquant embarquement_manquant
1    891          177              0                     2

L’âge manque pour 177 passagers. Pour garder le code du jour 2 lisible, le modèle commun utilise les 714 lignes complètes pour l’âge et les autres prédicteurs retenus.

titanic_analyse <- titanic_eiom |>
  drop_na(age, tarif, sexe, classe)

nrow(titanic_analyse)

[1] 714

L’analyse en cas complets n’est pas neutre. Elle est acceptable ici comme simplification pédagogique, mais une analyse substantielle devrait étudier pourquoi l’âge manque et comparer des stratégies d’imputation. Une disparition de lignes peut modifier la population effectivement analysée.

5. Pourquoi ne pas ajuster une droite ordinaire?

Un modèle linéaire appliqué directement à une réponse 0 ou 1 pose trois problèmes:

1. ses valeurs ajustées peuvent sortir de l’intervalle \([0,1]\);
2. la variance d’une réponse binaire dépend de sa probabilité;
3. une différence fixe sur l’échelle linéaire ne représente pas bien la saturation près de 0 et de 1.

La régression logistique conserve un prédicteur linéaire, mais elle ne le relie pas directement à la probabilité. Elle passe par une fonction de lien.

6. Passer par trois échelles

Pour une probabilité \(p\):

\[ \text{cote} = \frac{p}{1-p}, \qquad \text{logit}(p) = \log\left(\frac{p}{1-p}\right). \]

La cote compare la probabilité de l’événement à celle de son complément. Une cote de 3 signifie trois possibilités favorables pour une possibilité défavorable. Elle correspond à une probabilité de \(3/(1+3)=0{,}75\).

tibble(
  probabilite = c(0.10, 0.25, 0.50, 0.75, 0.90)
) |>
  mutate(
    cote = probabilite / (1 - probabilite),
    log_cote = log(cote)
  )

# A tibble: 5 × 3
  probabilite  cote log_cote
        <dbl> <dbl>    <dbl>
1        0.1  0.111    -2.20
2        0.25 0.333    -1.10
3        0.5  1         0   
4        0.75 3         1.10
5        0.9  9         2.20

La probabilité est bornée entre 0 et 1. La cote est positive. La log-cote peut prendre n’importe quelle valeur réelle. C’est sur cette dernière échelle que le modèle reste linéaire.

7. Écrire le modèle logistique

Le modèle commun est:

\[ \log\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1\,\text{sexe}_i + \beta_2\,\text{classe}_i + \beta_3\,\text{âge}_i + \beta_4\,\text{tarif}_i. \]

La fonction logistique ramène ensuite le prédicteur linéaire \(\eta_i\) à une probabilité:

\[ p_i = \frac{\exp(\eta_i)}{1+\exp(\eta_i)}. \]

Le mot « linéaire » décrit la relation entre les prédicteurs et la log-cote, pas nécessairement la relation entre un prédicteur et la probabilité.

8. Séparer entraînement et test

set.seed(2026)

partage <- initial_split(
  titanic_analyse,
  prop = 0.75,
  strata = verite
)

entrainement <- training(partage)
test <- testing(partage)

tibble(
  ensemble = c("Entraînement", "Test"),
  n = c(nrow(entrainement), nrow(test)),
  proportion_survecu = c(
    mean(entrainement$survie),
    mean(test$survie)
  )
)

# A tibble: 2 × 3
  ensemble         n proportion_survecu
  <chr>        <int>              <dbl>
1 Entraînement   535              0.406
2 Test           179              0.408

Le modèle apprend ses coefficients sur 535 lignes. Les 179 lignes de test servent à évaluer les probabilités sur des observations qui n’ont pas contribué à l’ajustement.

La stratification stabilise la proportion de l’événement dans les deux ensembles. Le germe fixe rend le partage reproductible. Le test doit rester fermé pendant le choix de la formule, sinon il cesse d’être une évaluation honnête (M. Kuhn et J. Silge, 2022).

9. Estimer par vraisemblance

En régression linéaire ordinaire, les moindres carrés minimisent une somme de résidus au carré. En régression logistique, les coefficients sont généralement estimés par maximum de vraisemblance.

L’idée est simple:

choisir les coefficients qui rendent les réponses 0 et 1 observées aussi plausibles que possible sous le modèle.

Une observation survivante contribue davantage à la vraisemblance si sa probabilité prédite est élevée. Une observation non survivante contribue davantage si cette probabilité est faible. Le calcul numérique est effectué par glm(); aucune dérivation n’est nécessaire pour interpréter le résultat.

10. Ajuster avec glm()

modele_logistique <- glm(
  survie ~ sexe + classe + age + tarif,
  data = entrainement,
  family = binomial()
)

summary(modele_logistique)

Call:
glm(formula = survie ~ sexe + classe + age + tarif, family = binomial(), 
    data = entrainement)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.268137   0.286928  -4.420 9.88e-06 ***
sexefemme    2.438833   0.238654  10.219  < 2e-16 ***
classe2e     1.337214   0.280871   4.761 1.93e-06 ***
classe1re    2.496966   0.385640   6.475 9.49e-11 ***
age         -0.037415   0.008986  -4.164 3.13e-05 ***
tarif        0.001997   0.002712   0.736    0.462    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 722.49  on 534  degrees of freedom
Residual deviance: 488.24  on 529  degrees of freedom
AIC: 500.24

Number of Fisher Scoring iterations: 5

L’argument family = binomial() indique une réponse binaire et utilise le lien logit par défaut. Les catégories de référence sont homme pour le sexe et 3e pour la classe.

La sortie summary() présente les coefficients sur l’échelle des log-cotes. Cette échelle est essentielle au modèle, mais elle est rarement la plus intuitive pour communiquer une association.

11. Lire les rapports de cotes

L’exponentielle d’un coefficient produit un rapport de cotes:

\[ \operatorname{OR}_j = \exp(\beta_j). \]

rapports_cotes <- tidy(modele_logistique) |>
  mutate(
    estimation_log = estimate,
    erreur_type_log = std.error,
    conf.low = exp(estimate - qnorm(0.975) * std.error),
    conf.high = exp(estimate + qnorm(0.975) * std.error),
    estimate = exp(estimate)
  ) |>
  select(
    term, estimation_log, erreur_type_log,
    estimate, conf.low, conf.high, p.value
  )

rapports_cotes

# A tibble: 6 × 7
  term       estimation_log erreur_type_log estimate conf.low conf.high  p.value
  <chr>               <dbl>           <dbl>    <dbl>    <dbl>     <dbl>    <dbl>
1 (Intercep…       -1.27            0.287      0.281    0.160     0.494 9.88e- 6
2 sexefemme         2.44            0.239     11.5      7.18     18.3   1.63e-24
3 classe2e          1.34            0.281      3.81     2.20      6.60  1.93e- 6
4 classe1re         2.50            0.386     12.1      5.70     25.9   9.49e-11
5 age              -0.0374          0.00899    0.963    0.946     0.980 3.13e- 5
6 tarif             0.00200         0.00271    1.00     0.997     1.01  4.62e- 1

Les intervalles affichés sont des intervalles de Wald, calculés sur l’échelle des log-cotes puis exponentiés. Ce choix les aligne avec les tests z de Wald de summary(). Des intervalles fondés sur la vraisemblance profilée peuvent être légèrement différents. La colonne erreur_type_log reste explicitement sur l’échelle des log-cotes et ne doit pas être combinée directement avec le rapport de cotes exponentié.

Pour sexefemme, le rapport compare femmes et hommes, à classe, âge et tarif identiques. Pour classe1re, il compare 1re et 3e classe, toutes les autres variables du modèle étant maintenues constantes.

Vérifier ce que « identiques » permet réellement

Fare est documenté comme le tarif du passager dans le paquet. Dans ce fichier, plusieurs personnes portant le même numéro de billet possèdent souvent exactement la même valeur. Cette variable ne doit donc pas être interprétée automatiquement comme une mesure de richesse personnelle ou comme un coût individuel indépendant du groupe de voyage.

titanic_train |>
  clean_names() |>
  count(ticket, fare) |>
  filter(n > 1) |>
  summarise(groupes_billet_tarif_repetes = n())

  groupes_billet_tarif_repetes
1                          133

La classe et le tarif sont liés. On peut vérifier à la fois la redondance globale et la sensibilité du rapport de cotes de la classe à la présence du tarif:

modele_sans_tarif <- update(
  modele_logistique,
  . ~ . - tarif
)

comparaison_classe <- bind_rows(
  tidy(modele_logistique) |>
    mutate(modele = "avec tarif"),
  tidy(modele_sans_tarif) |>
    mutate(modele = "sans tarif")
) |>
  filter(term %in% c("classe2e", "classe1re")) |>
  transmute(
    modele,
    comparaison = recode(
      term,
      classe2e = "2e contre 3e",
      classe1re = "1re contre 3e"
    ),
    rapport_cotes = exp(estimate)
  )

comparaison_classe

# A tibble: 4 × 3
  modele     comparaison   rapport_cotes
  <chr>      <chr>                 <dbl>
1 avec tarif 2e contre 3e           3.81
2 avec tarif 1re contre 3e         12.1 
3 sans tarif 2e contre 3e           3.90
4 sans tarif 1re contre 3e         14.1 

check_collinearity(modele_logistique)

# Check for Multicollinearity

Low Correlation

   Term  VIF   VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
   sexe 1.07 [1.02, 1.27]     1.03      0.94     [0.79, 0.98]
 classe 1.98 [1.76, 2.26]     1.19      0.51     [0.44, 0.57]
    age 1.42 [1.30, 1.61]     1.19      0.70     [0.62, 0.77]
  tarif 1.39 [1.26, 1.57]     1.18      0.72     [0.64, 0.79]

Ici, les VIF restent faibles à modérés, mais cela ne garantit pas que toutes les combinaisons classe-tarif soient abondantes. Le rapport de cotes de la 1re classe est conditionnel à un tarif donné et passe d’environ 14,1 sans tarif à 12,1 avec tarif. Ce résultat doit être communiqué comme une comparaison de modèle dans la zone couverte par les données, pas comme un effet causal de la classe.

Un rapport de cotes:

• supérieur à 1 indique des cotes plus élevées;
• inférieur à 1 indique des cotes plus faibles;
• égal à 1 indique l’absence de différence conditionnelle sur cette échelle.

Il ne s’agit ni d’un rapport de probabilités, ni d’une différence de points de pourcentage, ni d’une preuve causale.

12. Adapter l’unité d’une variable continue

Le rapport de cotes associé à l’âge porte sur une année. Pour une différence de dix ans:

or_age_10 <- exp(10 * coef(modele_logistique)[["age"]])
or_age_10

[1] 0.6878716

La formulation correcte est:

À sexe, classe et tarif identiques, dix années supplémentaires sont associées à des cotes de survie enregistrée multipliées par la valeur calculée.

Cette interprétation suppose une relation linéaire entre l’âge et la log-cote. Elle ne garantit pas que la même pente soit adéquate pour les très jeunes enfants et les adultes plus âgés.

13. Distinguer coefficient, intervalle et valeur p

Le coefficient décrit une direction et une taille sur l’échelle logit. Son exponentielle donne le rapport de cotes. L’intervalle de confiance décrit l’incertitude d’estimation sous les hypothèses du modèle. La valeur p examine la compatibilité avec l’absence d’association conditionnelle.

Une petite valeur p ne garantit pas:

1. une grande différence de probabilité;
2. une bonne prédiction hors échantillon;
3. une bonne calibration;
4. une importance historique ou pratique;
5. une interprétation causale.

Inversement, une incertitude importante peut venir d’un échantillon restreint, de groupes rares, de redondance entre prédicteurs ou d’une forme mal spécifiée.

14. Revenir aux probabilités

Deux profils permettent de voir pourquoi un rapport de cotes n’est pas une différence de probabilités.

profils <- tibble(
  profil = c("Homme, 3e classe", "Femme, 3e classe"),
  sexe = factor(
    c("homme", "femme"),
    levels = levels(entrainement$sexe)
  ),
  classe = factor(
    c("3e", "3e"),
    levels = levels(entrainement$classe)
  ),
  age = c(30, 30),
  tarif = c(15, 15)
)

profils |>
  mutate(
    probabilite_survie = predict(
      modele_logistique,
      newdata = profils,
      type = "response"
    )
  ) |>
  select(profil, probabilite_survie)

# A tibble: 2 × 2
  profil           probabilite_survie
  <chr>                         <dbl>
1 Homme, 3e classe             0.0862
2 Femme, 3e classe             0.520 

La variation de probabilité dépend de tous les termes qui déterminent le profil de départ. Un même rapport de cotes peut correspondre à une petite différence près de 0, à une différence plus grande au centre de l’échelle, puis à une petite différence près de 1.

15. Produire les probabilités du test

predictions_test <- test |>
  mutate(
    probabilite_survie = predict(
      modele_logistique,
      newdata = test,
      type = "response"
    )
  )

predictions_test |>
  select(
    passenger_id,
    verite,
    sexe,
    classe,
    age,
    probabilite_survie
  ) |>
  slice_head(n = 10)

   passenger_id      verite  sexe classe   age probabilite_survie
1             2     survecu femme    1re 38.00          0.9159234
2            13 non_survecu homme     3e 20.00          0.1191641
3            25 non_survecu femme     3e  8.00          0.7137060
4            36 non_survecu homme    1re 42.00          0.4405844
5            50 non_survecu femme     3e 18.00          0.6301266
6            64 non_survecu homme     3e  4.00          0.2039003
7            71 non_survecu homme     2e 32.00          0.2483842
8            79     survecu homme     2e  0.83          0.5239625
9            81 non_survecu homme     3e 22.00          0.1117187
10           85     survecu femme     2e 17.00          0.8690760

Une probabilité élevée n’est pas une certitude. Elle signifie que le profil reçoit une valeur plus proche de 1 sous le modèle ajusté. La qualité de ces probabilités doit être vérifiée sur le test.

16. Transformer une probabilité en classe

Un seuil \(s\) définit une règle:

\[ \widehat{Y} = \begin{cases} 1, & \widehat{p} \ge s,\\ 0, & \widehat{p} < s. \end{cases} \]

predictions_05 <- predictions_test |>
  mutate(
    classe_predite = factor(
      if_else(
        probabilite_survie >= 0.50,
        "survecu",
        "non_survecu"
      ),
      levels = levels(verite)
    )
  )

conf_mat(
  predictions_05,
  truth = verite,
  estimate = classe_predite
)

             Truth
Prediction    non_survecu survecu
  non_survecu          85      15
  survecu              21      58

Le seuil de 0,50 est facile à expliquer, mais il n’est pas une propriété naturelle de la régression logistique. Il donne implicitement le même poids aux deux catégories d’erreur et ignore les contraintes du mandat.

17. Lire les mesures de classification

metric_set(accuracy, sens, spec, ppv, npv)(
  predictions_05,
  truth = verite,
  estimate = classe_predite,
  event_level = "second"
)

# A tibble: 5 × 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy binary         0.799
2 sens     binary         0.795
3 spec     binary         0.802
4 ppv      binary         0.734
5 npv      binary         0.85 

Les mesures répondent à des questions différentes:

Mesure	Question
Exactitude	Quelle proportion de toutes les classes est correcte?
Sensibilité	Parmi les survivants enregistrés, quelle proportion est détectée?
Spécificité	Parmi les non-survivants enregistrés, quelle proportion est écartée?
Précision	Parmi les prédictions positives, quelle proportion correspond à un survivant enregistré?
Valeur prédictive négative	Parmi les prédictions négatives, quelle proportion correspond à un non-survivant enregistré?

La sensibilité et la précision ont des dénominateurs différents. Les confondre change la question.

18. Comparer les seuils

evaluer_seuil <- function(donnees, seuil) {
  donnees |>
    summarise(
      seuil = seuil,
      proportion_positive = mean(probabilite_survie >= seuil),
      sensibilite = sum(
        probabilite_survie >= seuil & survie == 1
      ) / sum(survie == 1),
      specificite = sum(
        probabilite_survie < seuil & survie == 0
      ) / sum(survie == 0),
      precision = sum(
        probabilite_survie >= seuil & survie == 1
      ) / sum(probabilite_survie >= seuil)
    )
}

table_seuils <- map_dfr(
  seq(0.10, 0.90, by = 0.05),
  ~ evaluer_seuil(predictions_test, .x)
)

table_seuils |>
  filter(seuil %in% c(0.30, 0.50, 0.70))

  seuil proportion_positive sensibilite specificite precision
1   0.5           0.4413408   0.7945205   0.8018868 0.7341772

table_seuils |>
  select(seuil, sensibilite, specificite, precision) |>
  pivot_longer(
    -seuil,
    names_to = "mesure",
    values_to = "valeur"
  ) |>
  ggplot(aes(seuil, valeur, color = mesure)) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = label_percent()) +
  labs(x = "Seuil", y = "Valeur", color = NULL)

Lorsque le seuil augmente, les prédictions positives deviennent moins nombreuses. La sensibilité diminue généralement, tandis que la spécificité et parfois la précision augmentent. Le seuil doit donc être relié à une erreur plus grave, à une capacité ou à un usage analytique.

Dans le contexte Titanic, ce choix sert uniquement à comprendre la mécanique de classification. Il ne reconstitue pas une décision de sauvetage.

19. Mesurer la discrimination

roc_auc(
  predictions_test,
  truth = verite,
  probabilite_survie,
  event_level = "second"
)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.863

pr_auc(
  predictions_test,
  truth = verite,
  probabilite_survie,
  event_level = "second"
)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 pr_auc  binary         0.823

roc_curve(
  predictions_test,
  truth = verite,
  probabilite_survie,
  event_level = "second"
) |>
  autoplot() +
  coord_equal()

L’aire ROC mesure la probabilité qu’un survivant enregistré choisi au hasard reçoive un score plus élevé qu’un non-survivant enregistré choisi au hasard (T. Fawcett, 2006). Elle évalue le classement sur tous les seuils.

Une aire ROC élevée ne garantit ni des probabilités exactes, ni un seuil utile, ni une validité causale. La courbe précision-rappel complète l’analyse lorsque l’événement est rare ou que les faux positifs sont particulièrement importants.

20. Examiner la calibration

La calibration compare probabilités prédites et fréquences observées. Si un ensemble de profils reçoit une probabilité moyenne de 0,70, une bonne calibration signifie qu’environ 70 % de ces observations connaissent l’événement, dans des données comparables.

calibration <- predictions_test |>
  mutate(groupe = ntile(probabilite_survie, 5)) |>
  summarise(
    probabilite_moyenne = mean(probabilite_survie),
    frequence_observee = mean(survie),
    n = n(),
    .by = groupe
  )

ggplot(
  calibration,
  aes(probabilite_moyenne, frequence_observee)
) +
  geom_abline(
    slope = 1,
    intercept = 0,
    linetype = 2,
    color = "#52606d"
  ) +
  geom_line(color = "#0f7c80") +
  geom_point(aes(size = n), color = "#b5222e") +
  coord_equal(xlim = c(0, 1), ylim = c(0, 1)) +
  labs(
    x = "Probabilité moyenne prédite",
    y = "Fréquence observée",
    size = "n"
  )

La représentation groupée est pédagogique et instable avec seulement 179 lignes de test. Elle permet de repérer des écarts grossiers, pas de déclarer définitivement le modèle calibré.

Le score de Brier est l’erreur quadratique moyenne des probabilités. Pour lui donner un point de comparaison, on peut utiliser une règle constante égale à la proportion de survivants enregistrés dans l’entraînement:

probabilite_reference <- mean(entrainement$survie)

comparaison_brier <- predictions_test |>
  summarise(
    brier_modele = mean(
      (survie - probabilite_survie)^2
    ),
    brier_reference = mean(
      (survie - probabilite_reference)^2
    )
  )

comparaison_brier

  brier_modele brier_reference
1    0.1391384        0.241508

Un score de 0 correspondrait à des probabilités parfaites. Dans ce partage reproductible, le modèle obtient environ 0,139 contre 0,242 pour la règle constante. Cette amélioration ne suffit pas à démontrer une bonne calibration dans chaque partie de l’échelle. Le score combine discrimination et calibration; il doit être comparé à une référence et interprété avec la prévalence (B. Van Calster et al., 2019).

21. Conditions de validité

Une régression logistique simple demande plus qu’un appel réussi à glm().

Condition	Pourquoi elle importe	Vérification possible
Réponse et unité correctement définies	Une mauvaise cible rend le modèle incohérent	documentation et protocole
Observations suffisamment indépendantes	La dépendance réduit l’information effective	plan de collecte et structure des groupes
Forme adéquate sur l’échelle logit	Une pente linéaire peut masquer une courbure	résidus, lissages et extensions
Absence de séparation complète	Des groupes parfaitement séparés rendent certains coefficients instables	avertissements, cellules vides et coefficients extrêmes
Prédicteurs sans redondance extrême	La colinéarité élargit l’incertitude	tableaux, VIF et stabilité
Nombre suffisant d’événements	Les estimations et mesures doivent être stables	effectifs par groupe et intervalles
Données pertinentes pour l’usage	Une bonne performance interne ne garantit pas le transfert	validation externe et comparaison des populations

La normalité des prédicteurs n’est pas une hypothèse de la régression logistique. La normalité des résidus n’est pas exigée comme en régression linéaire gaussienne.

22. Diagnostics et extensions

Des résidus de déviance importants peuvent signaler des observations mal représentées. Le levier et la distance de Cook aident à repérer les cas qui influencent fortement les coefficients. Aucun seuil automatique ne justifie une suppression.

diagnostic <- augment(modele_logistique) |>
  transmute(
    ligne = row_number(),
    valeur_ajustee = .fitted,
    residu_deviance = .resid,
    levier = .hat,
    distance_cook = .cooksd
  )

diagnostic |>
  slice_max(distance_cook, n = 8)

# A tibble: 8 × 5
  ligne valeur_ajustee residu_deviance  levier distance_cook
  <int>          <dbl>           <dbl>   <dbl>         <dbl>
1   482          0.905           0.824 0.298          0.0409
2   494          0.942           0.811 0.293          0.0381
3    10          1.04           -1.64  0.0626         0.0337
4    51          0.825          -1.54  0.0586         0.0252
5   177          3.03           -2.48  0.00616        0.0216
6   470         -1.70            1.93  0.0224         0.0215
7   457         -2.23            2.16  0.0126         0.0201
8    75          1.85           -2.00  0.0120         0.0131

Deux extensions naturelles sont:

modele_interaction <- glm(
  survie ~ sexe * classe + age + tarif,
  data = entrainement,
  family = binomial()
)

modele_age_courbe <- glm(
  survie ~ sexe + classe + splines::ns(age, df = 3) + tarif,
  data = entrainement,
  family = binomial()
)

L’interaction autorise l’association du sexe à varier selon la classe. La spline autorise une relation non linéaire entre l’âge et la log-cote. Ces extensions doivent être justifiées, comparées hors échantillon et communiquées sans masquer leur complexité.

23. Ce que le modèle ne permet pas d’affirmer

Ce prototype ne permet pas de conclure:

1. qu’un prédicteur a causé la survie ou le décès;
2. que les coefficients décrivent tous les passagers du Titanic;
3. que la suppression des âges manquants est sans conséquence;
4. qu’un seuil de 0,50 est historiquement ou moralement optimal;
5. que la performance resterait identique dans une autre population;
6. que les variables disponibles résument toutes les circonstances.

Le jeu est utile parce que les variables sont faciles à expliquer. Cette accessibilité ne doit pas faire oublier la portée limitée de l’échantillon et la gravité du contexte historique.

24. Sept questions avant de défendre le modèle

1. Quelle est l’unité statistique et quel événement est positif?
2. Sur quelle population les données permettent-elles de conclure?
3. Comment les valeurs manquantes ont-elles été traitées?
4. Que signifie chaque rapport de cotes et quelle est sa catégorie de référence?
5. La performance a-t-elle été évaluée sur des données réservées?
6. Le modèle discrimine-t-il et calibre-t-il suffisamment pour l’usage visé?
7. Quel seuil, quelle erreur et quelle limite accompagnent la conclusion?

Lectures complémentaires

Lecture	Pour approfondir
James et ses collègues, chapitre 4	Régression logistique, classification et interprétation.
Documentation officielle de glm()	Syntaxe, famille, prédictions et objet retourné.
Familles de modèles dans R	Lien logit, réponse binomiale et codage de l’événement.
Kuhn et Silge, séparation des données	Jeu d’entraînement, jeu test et budget de données.
Fawcett, courbes ROC	Sensibilité, spécificité et aire ROC.
Van Calster et ses collègues, calibration	Calibration, discrimination et utilité prédictive.
Documentation CRAN de titanic	Structure, variables et provenance du jeu.

Conclusion

Une régression logistique ne se résume pas à transformer une réponse en 0 et 1. Elle relie une probabilité conditionnelle à des prédicteurs, impose une échelle d’interprétation, demande une validation hors échantillon et ne produit une classe qu’après un choix de seuil.

Une conclusion défendable précise donc la population, les catégories de référence, l’incertitude, la discrimination, la calibration, le seuil et les limites.

Les références

T. Fawcett. (2006). An Introduction to ROC Analysis. Pattern Recognition Letters. https://doi.org/10.1016/j.patrec.2005.10.010.

P. Hendricks. (2015). titanic: Titanic Passenger Survival Data Set. https://CRAN.R-project.org/package=titanic.

M. Kuhn et J. Silge. (2022). Tidy Modeling with R. https://www.tmwr.org/.

B. Van Calster, D. J. McLernon, M. van Smeden, L. Wynants, et E. W. Steyerberg. (2019). Calibration: The Achilles Heel of Predictive Analytics. BMC Medicine. https://doi.org/10.1186/s12916-019-1466-7.