Mission 1 — Titanic: modéliser la survie (logistique)

source("../../utils/requirements.R")
library(titanic); library(dplyr); library(broom); library(ggplot2); library(janitor)
data("titanic_train")
df <- titanic_train %>%
  clean_names() %>%
  mutate(
    survived = factor(survived, levels=c(0,1), labels=c("No","Yes")),
    sex = factor(sex),
    pclass = factor(pclass),
    embarked = factor(embarked)
  ) %>%
  select(survived, pclass, sex, age, sib_sp, parch, fare, embarked)
set.seed(123)

Objectifs

EDA rapide — comprendre les prédicteurs clés et le sens attendu des effets
Modèle de régression logistique glm(..., family = binomial)
Interprétation des coefficients via odds ratios (rapports de cotes)
Évaluation de base : prédictions, matrice de confusion (seuil 0.5), sensibilité/spécificité
Préparer le terrain pour M2 (ROC, AUC, choix de seuil)

1) EDA éclair

summary(df)

 survived  pclass      sex           age            sib_sp     
 No :549   1:216   female:314   Min.   : 0.42   Min.   :0.000  
 Yes:342   2:184   male  :577   1st Qu.:20.12   1st Qu.:0.000  
           3:491                Median :28.00   Median :0.000  
                                Mean   :29.70   Mean   :0.523  
                                3rd Qu.:38.00   3rd Qu.:1.000  
                                Max.   :80.00   Max.   :8.000  
                                NA's   :177                    
     parch             fare        embarked
 Min.   :0.0000   Min.   :  0.00    :  2   
 1st Qu.:0.0000   1st Qu.:  7.91   C:168   
 Median :0.0000   Median : 14.45   Q: 77   
 Mean   :0.3816   Mean   : 32.20   S:644   
 3rd Qu.:0.0000   3rd Qu.: 31.00           
 Max.   :6.0000   Max.   :512.33

ggplot(df, aes(sex, fill=survived)) + geom_bar(position="fill") + labs(y="Proportion", title="Survie par sexe")

ggplot(df, aes(pclass, fill=survived)) + geom_bar(position="fill") + labs(y="Proportion", title="Survie par classe")

ggplot(df, aes(age, fill=survived)) + geom_histogram(alpha=.6, bins=30, position="identity") + labs(title="Âge: distribution par statut de survie")

Note

À faire (équipe)
- Quelles variables semblent fortement associées à la survie ?
- Y a-t-il des valeurs manquantes importantes (ex. age) ? Comment les géreriez-vous pour un premier modèle simple (supprimer lignes manquantes vs imputation grossière) ?

Astuce : pour aller vite, on peut filtrer les NA dans un premier temps, puis discuter des alternatives plus rigoureuses plus tard.

df_complete <- df %>% tidyr::drop_na()

2) Modèle logistique de base

mod0 <- glm(survived ~ sex + pclass + age + fare, data=df_complete, family=binomial())
summary(mod0)


Call:
glm(formula = survived ~ sex + pclass + age + fare, family = binomial(), 
    data = df_complete)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.7225052  0.4645113   8.014 1.11e-15 ***
sexmale     -2.5185052  0.2082017 -12.096  < 2e-16 ***
pclass2     -1.2765903  0.3126370  -4.083 4.44e-05 ***
pclass3     -2.5415762  0.3277677  -7.754 8.89e-15 ***
age         -0.0367302  0.0077325  -4.750 2.03e-06 ***
fare         0.0005226  0.0022579   0.231    0.817    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 964.52  on 713  degrees of freedom
Residual deviance: 647.23  on 708  degrees of freedom
AIC: 659.23

Number of Fisher Scoring iterations: 5

Tip

Interprétation rapide
- Le signe de \(\\beta\) indique si la variable augmente ou diminue l’odds de survie.
- Pour une interprétation multiplicative, calculez \(\\exp(\\beta)\) : c’est le odds ratio.

tidy(mod0, exponentiate = TRUE, conf.int = TRUE) %>%
  mutate(across(c(estimate, conf.low, conf.high), round, 3)) %>%
  knitr::kable(caption="Odds ratios (exp(coefficients)) avec IC à 95%")

Odds ratios (exp(coefficients)) avec IC à 95%
term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	41.368	0.4645113	8.0138092	0.0000000	16.939	105.090
sexmale	0.081	0.2082017	-12.0964697	0.0000000	0.053	0.120
pclass2	0.279	0.3126370	-4.0832994	0.0000444	0.150	0.512
pclass3	0.079	0.3277677	-7.7541995	0.0000000	0.041	0.148
age	0.964	0.0077325	-4.7501079	0.0000020	0.949	0.978
fare	1.001	0.0022579	0.2314361	0.8169761	0.996	1.005

Note

Réflexion
- Interprétez l’odds ratio de sex et pclass. Est-ce cohérent avec le contexte historique ?
- age est-il cliniquement important même si son effet paraît plus faible ?

3) Prédictions et matrice de confusion (seuil 0.5)

df_complete <- df_complete %>% mutate(prob = predict(mod0, type="response"),
                                      pred = ifelse(prob >= 0.5, "Yes", "No") %>% factor(levels=c("No","Yes")))
tab <- table(Truth = df_complete$survived, Pred = df_complete$pred)
acc <- sum(diag(tab))/sum(tab)
sens <- tab["Yes","Yes"]/sum(tab[,"Yes"])        # Sensibilité (rappel des 'Yes')
spec <- tab["No","No"]/sum(tab[,"No"])           # Spécificité
knitr::kable(as.data.frame.matrix(tab), caption=sprintf("Matrice de confusion (seuil 0.5) — Acc=%.3f, Sens=%.3f, Spec=%.3f", acc, sens, spec))

Matrice de confusion (seuil 0.5) — Acc=0.793, Sens=0.757, Spec=0.815
	No	Yes
No	357	67
Yes	81	209

Warning

Attention — piège classique
Le seuil 0.5 n’est pas toujours pertinent. Sa pertinence dépend : (i) de la prévalence des classes, (ii) des coûts d’erreurs (faux positifs/négatifs).

Note

Question
- Si rater un survivant (faux négatif) coûte plus cher que classer à tort comme survivant (faux positif), voudriez-vous un seuil plus bas ou plus haut ? Expliquez.

4) (optionnel) - Effets non linéaires et interactions (aperçu)

mod1 <- glm(survived ~ sex * pclass + splines::ns(age, df = 3) + fare, data=df_complete, family=binomial())
broom::glance(mod0); broom::glance(mod1)

# A tibble: 1 × 8
  null.deviance df.null logLik   AIC   BIC deviance df.residual  nobs
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
1          965.     713  -324.  659.  687.     647.         708   714

# A tibble: 1 × 8
  null.deviance df.null logLik   AIC   BIC deviance df.residual  nobs
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
1          965.     713  -304.  628.  673.     608.         704   714

Tip

À discuter
- Pourquoi une spline sur age ? (relation non linéaire plausible)
- L’interaction sex:pclass a-t-elle un sens ?

Travail d’équipe

Rédigez deux interprétations d’odds ratios (une pour une variable binaire, une pour age continu).
Donnez un argument pour/contre l’inclusion de l’interaction sex:pclass.
Proposez une modification du modèle (transformation, interaction, variable) et justifiez-la.

Points de discussion (retour groupe)

Différence probabilité vs odds vs log-odds.
Pourquoi l’odds ratio est multiplicatif et dépend du niveau de référence.
Limites : séparation quasi-parfaite, classes rares, gestion des NA.

Questions bonus

Montrez que \(\\text{logit}(p)=\\log\\frac{p}{1-p}\) est une fonction croissante de \(p\).
Calculez la variation d’odds pour +10 ans d’âge (utilisez \(\\exp(10\\beta_\\text{age})\)).
Donnez un exemple où un coefficient est significatif statistiquement mais peu utile pour la décision au seuil 0.5.