Vendredi dernier, j’ai mis en ligne un petit quiz pour prédire qui a survécu au naufrage du Titanic, en 1912. J’avais mis en ligne une sous-base avec quelques passagers
loc_fichier = "http://freakonometrics.free.fr/titanic.RData" download.file(loc_fichier, "titanic.RData") load("titanic.RData") str(base) |
et plus précisément
Survived
: Passenger survival indicator (1 if survived)Pclass
: Passenger classSex
: Sex of the passengerAge
: Age of the passengerSibSp
: Number of siblings/spouses aboardParch
: Number of parents/children aboardEmbarked
: Port of EmbarkationName
: Name of passenger
La dernière variable ne sert pas à grand chose, donc on la vire,
base = base[,1:7] |
On va maintenant pouvoir répondre aux questions.
- quelle proportion de passagers ont survécu ?
La réponse simple est
mean(base$Survived) [1] 0.3838384 |
que l’on peut retrouver dans le tableau de contingence suivant
table(base$Survived)/nrow(base) 0 1 0.6161616 0.3838384 |
soit 38.38\% de survivants ici. Cela dit, on peut aussi l’obtenir en faisant une régression logistique sur aucune variable explicative (autrement dit en régressant seulement sur la constante). En effet, la régression donne ici
glm(Survived~1, family=binomial,data=base) Coefficients: (Intercept) -0.4733 Degrees of Freedom: 890 Total (i.e. Null); 890 Residual Null Deviance: 1187 Residual Deviance: 1187 AIC: 1189 |
ce qui donne la valeur de \widehat{\beta}_0, et comme la probabilité de survie est \mathbb{P}[Y=1]=\frac{\exp(\beta_0)}{1+\exp(\beta_0)}on en déduit une prévision de cette probabilité en considérant\frac{\exp(\widehat{\beta}_0)}{1+\exp(\widehat{\beta}_0)}soit ici
exp(-0.4733)/(1+exp(-0.4733)) [1] 0.3838355 |
(l’erreur de précision vient du fait que j’utilise ici une version très approchée de \widehat{\beta}_0)On peut aussi utiliser la fonction predict
sous R (et prendre n’importe quelle prévision car elles sont ici toutes égales)
predict(glm(Survived~1, family=binomial, data=base), type="response")[1] 1 0.3838384 |
Si on note \widehat{\mu}_i=\frac{\exp(\widehat{\beta}_0)}{1+\exp(\widehat{\beta}_0)}, on retrouve ici le fait que pour une régression logistique, \frac{1}{n}\sum_{i-1}^n\widehat{\mu}_i=\frac{1}{n}\sum_{i-1}^n y_iCette propriété (qui était déjà vraie pour la méthode des moindres carrés) vient de la condition du premier ordre, lorsqu’on cherche à maximiser la (log)-vraisemblance. En effet, on l’avait vu au retour de la semaine de relâche, la log-vraisemblance est ici (je me place généralement dans le cas d’une régression sur des variables \boldsymbol{x},\log\mathcal{L}=\sum_{i=1}^ny_i \log p(\boldsymbol{x}_i)+(1-y_i) \log(1-p(\boldsymbol{x}_i))soit\log\mathcal{L}=\sum_{i=1}^ny_i[\boldsymbol{x}_i^\top\boldsymbol{\beta}]-\log(1+\exp[\boldsymbol{x}_i^\top\boldsymbol{\beta}])et la condition du premier ordre donne ici\left.\frac{\partial\log\mathcal{L}}{\partial \beta_j}\right\vert_{\boldsymbol{\beta}=\widehat{\boldsymbol{\beta}}}=0=\sum_{i=1}^n [y_i-\widehat{\mu}_i] x_{i,j}où \widehat{\mu}_i=\frac{\exp(\boldsymbol{x}_i^\top\widehat{\boldsymbol{\beta}})}{1+\exp(\boldsymbol{x}_i^\top\widehat{\boldsymbol{\beta}})}donc dans le cas particulier de la première condition du premier ordre, on obtient \left.\frac{\partial\log\mathcal{L}}{\partial \beta_0}\right\vert_{\boldsymbol{\beta}=\widehat{\boldsymbol{\beta}}}=0=\sum_{i=1}^n [y_i-\widehat{\mu}_i]et donc\frac{1}{n}\sum_{i-1}^n\widehat{\mu}_i=\frac{1}{n}\sum_{i-1}^n y_iD’ailleurs on peut noter que cette propriété reste vraie avec une régression probit,
reg=glm(Survived~1, family=binomial(link="probit"),data=base) > predict(reg,type="response")[1] 1 0.3838384 |
-
- quelle proportion de passagers de première classe ont survécu ?
Cette fois, on commence un peu à conditionner: on ne regarde que les personnes en première classe,
mean(base$Survived[base$Pclass=="1"]) [1] 0.6296296 |
Environ 63\% ont survécu. Mais là encore, on peut faire une régression logistique
glm(Survived~Pclass, family=binomial,data=base) Coefficients: (Intercept) Pclass2 Pclass3 0.5306 -0.6394 -1.6704 Degrees of Freedom: 890 Total (i.e. Null); 888 Residual Null Deviance: 1187 Residual Deviance: 1083 AIC: 1089 |
Comme la classe 1 est la classe de référence, on fait comme auparavant, en considérant\frac{\exp(\widehat{\beta}_0)}{1+\exp(\widehat{\beta}_0)}soit ici
exp(0.5306)/(1+exp(0.5306)) [1] 0.629623 |
(avec toujours la petite imprécision numérique). Une prévision plus propre donne ici
reg=glm(Survived~Pclass, family=binomial,data=base) predict(reg,newdata=data.frame(Pclass="1"),type="response") 1 0.6296296 |
Et là encore, on peut tenter une régression probit, et on obtient la même chose,
reg=glm(Survived~Pclass, family=binomial(link="probit"),data=base) predict(reg,newdata=data.frame(Pclass="1"),type="response") 1 0.6296296 |
-
- quelle est la statistique de test pour un test d’indépendance du chi-deux entre le fait de survivre (ou pas) et la classe ?
La commande pour le test du chi-deux est la suivante
chisq.test(table(base$Survived,base$Pclass)) Pearson's Chi-squared test data: table(base$Survived, base$Pclass) X-squared = 102.89, df = 2, p-value < 2.2e-16 |
Juste un petit rappel (on a revu le test du chi-deux au dernier cours), on a un tableau de contingence, avec n_{i,j} de nombre de personnes de type i (pour la première variable) et j (pour la seconde). Si les variables étaient indépendantes, on aura n_{i,j}^\perp=n\frac{n_{i\cdot}}{n}\frac{n_{\cdot j}}{n}où\sum_{i=1}^n n_{i,j}=n_{\cdot j}~~~\sum_{j=1}^n n_{i,j}=n_{i\cdot}La statistique de test est alors \displaystyle{\chi ^{2}=\sum _{i=1}^{I}\sum _{j=1}^{J}\frac{(n_{i,j}-n^\perp_{i,j})^{2}}{n^\perp_{i,j}}}On peut noter que les contributions au test sont ici\varepsilon_{i,j}=\frac{n_{i,j}-n^\perp_{i,j}}{\sqrt{n^\perp_{i,j}}}
chisq.test(table(base$Survived,base$Pclass))$residuals 1 2 3 0 -4.601993 -1.537771 3.993703 1 5.830678 1.948340 -5.059981 |
ce qui nous donne beaucoup d’information : on observe deux fortes valeurs positives, correspondant à une forte association (positive) entre “survivre” et “première classe”, et “ne pas survivre” et “troisième classe”, et deux fortes valeurs négatives, correspondant à une forte répulsion entre “survivre” et “troisième classe”, et “ne pas survivre” et “première classe”. On peut visualiser ces valeurs sur le dessin ci-dessous
library("graphics") mosaicplot(table(base$Pclass,base$Survived), shade = TRUE, las=2) |
ou celui là,
library("vcd") assoc(table(base$Survived,base$Pclass), shade = TRUE, las=3) |
- ensuite il faut faire une régression logistique, et prévoir la probabilité de survie pour deux passagers (fictifs) du paquebot
On nous donne la base avec deux passagers
newbase = data.frame( Pclass = as.factor(c(1,3)), Sex = as.factor(c("female","male")), Age = c(17,20), SibSp = c(1,0), Parch = c(2,0), Embarked = as.factor(c("S","S")), Name = as.factor(c("Winslet, Miss. Kate","DiCaprio, Mr. Leonardo"))) |
Faisons une simple régression, avec toutes les variables,
reg=glm(Survived~., family=binomial,data=base) summary(reg) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 16.830381 607.655774 0.028 0.97790 Pclass2 -1.268362 0.298428 -4.250 2.14e-05 *** Pclass3 -2.493756 0.296219 -8.419 < 2e-16 *** Sexmale -2.641145 0.222801 -11.854 < 2e-16 *** Age -0.043725 0.008294 -5.272 1.35e-07 *** SibSp -0.355755 0.128529 -2.768 0.00564 ** Parch -0.044628 0.120705 -0.370 0.71159 EmbarkedC -12.260112 607.655693 -0.020 0.98390 EmbarkedQ -13.104581 607.655894 -0.022 0.98279 EmbarkedS -12.687791 607.655674 -0.021 0.98334 --- 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: 632.67 on 704 degrees of freedom (177 observations deleted due to missingness) AIC: 652.67 Number of Fisher Scoring iterations: 13 |
Deux variables ne sont pas significatives. On va les enlever, une à une, et on arrive au modèle suivant
reg=glm(Survived~Pclass+Sex+Age+SibSp, family=binomial,data=base) summary(reg) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 4.334201 0.450700 9.617 < 2e-16 *** Pclass2 -1.414360 0.284727 -4.967 6.78e-07 *** Pclass3 -2.652618 0.285832 -9.280 < 2e-16 *** Sexmale -2.627679 0.214771 -12.235 < 2e-16 *** Age -0.044760 0.008225 -5.442 5.27e-08 *** SibSp -0.380190 0.121516 -3.129 0.00176 ** --- 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: 636.56 on 708 degrees of freedom (177 observations deleted due to missingness) AIC: 648.56 Number of Fisher Scoring iterations: 5 |
Maintenant, comme souvent, quand on a des variables continues, comme l’âge, il n’y a aucune raison qu’il intervienne ici de manière linéaire. On peut par exemple faire une régression polynomiale
reg3=glm(Survived~Pclass+Sex+poly(Age,3)+SibSp, family=binomial,data=base[!is.na(base$Age),]) summary(reg3) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.0213 0.2903 10.408 < 2e-16 *** Pclass2 -1.3603 0.2842 -4.786 1.70e-06 *** Pclass3 -2.5569 0.2853 -8.962 < 2e-16 *** Sexmale -2.6582 0.2176 -12.216 < 2e-16 *** poly(Age, 3)1 -17.7668 3.2583 -5.453 4.96e-08 *** poly(Age, 3)2 6.0044 3.0021 2.000 0.045491 * poly(Age, 3)3 -5.9181 3.0992 -1.910 0.056188 . SibSp -0.5041 0.1317 -3.828 0.000129 *** --- 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: 627.55 on 706 degrees of freedom AIC: 643.55 Number of Fisher Scoring iterations: 5 |
mais interpréter les paramètres devient alors complexe (comme expliqué dan Function basis and regression). On note que le degré 3 est ici significatif, donc on va faire la régression à la main
reg3=glm(Survived~Pclass+Sex+Age+I(Age^2)+I(Age^3)+SibSp, family=binomial,data=base[!is.na(base$Age),]) > summary(reg3) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.616e+00 6.565e-01 8.554 < 2e-16 *** Pclass2 -1.360e+00 2.842e-01 -4.786 1.7e-06 *** Pclass3 -2.557e+00 2.853e-01 -8.962 < 2e-16 *** Sexmale -2.658e+00 2.176e-01 -12.216 < 2e-16 *** Age -1.905e-01 5.528e-02 -3.446 0.000569 *** I(Age^2) 4.290e-03 1.854e-03 2.314 0.020669 * I(Age^3) -3.520e-05 1.843e-05 -1.910 0.056188 . SibSp -5.041e-01 1.317e-01 -3.828 0.000129 *** --- 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: 627.55 on 706 degrees of freedom AIC: 643.55 Number of Fisher Scoring iterations: 5 |
On notera (comme je le disais dans Function basis and regression) que la p-value associée au degré le plus élevé est identique. Bref, la transformation de l’âge en une fonction de l’âge, non linéaire, a du sens. On peut même visualiser cette fonction
beta = coefficients(reg3) xage = 0:80 yage = beta[5]*xage+beta[6]*xage^2+beta[7]*xage^3 plot(xage,yage,xlab="Age",ylab="",type="l") |
En fait, on peut (comme pour les modèles linéaires “standards”) on peut aussi utiliser les splines. Et les gam sont un outils parfait en terme de visualisation
library(gam) regam=gam(Survived~Pclass+Sex+bs(Age)+SibSp, family=binomial,data=base) summary(regam) (Dispersion Parameter for binomial family taken to be 1) Null Deviance: 964.516 on 713 degrees of freedom Residual Deviance: 627.5525 on 706 degrees of freedom AIC: 643.5525 177 observations deleted due to missingness Number of Local Scoring Iterations: 4 Anova for Parametric Effects Df Sum Sq Mean Sq F value Pr(>F) Pclass 2 26.72 13.361 11.3500 1.407e-05 *** Sex 1 131.57 131.573 111.7678 < 2.2e-16 *** bs(Age) 3 22.76 7.588 6.4455 0.0002620 *** SibSp 1 14.66 14.659 12.4525 0.0004445 *** Residuals 706 831.10 1.177 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 |
On peut voir la transformation de la variable âge,
plot(reg) |
et on retrouve une transformation proche de notre polynôme de degré 3. En fait, on a même mieux ici parce qu’on peut rajouter des bandes de confiance, ce qui permet de vérifier que la fonction n’est pas vraiment linéaire
plot(reg,se=TRUE) |
Bref, on a maintenant trois modèles. Et (étrangement) les deux derniers donnent des prévisions sur nos deux personnes très proches,
predict(reg,newdata=newbase,type="response") 1 2 0.9605736 0.1368988 predict(reg3,newdata=newbase,type="response") 1 2 0.9497834 0.1218426 predict(regam,newdata=newbase,type="response") 1 2 0.9497834 0.1218426 |
Autrement dit, Leonardo DiCaprio avait de l’ordre de 12% de chances de survie (comme tenu de son âge, du fait qu’il avait un billet de troisième classe, et pas de famille sur le bateau).
OpenEdition suggests that you cite this post as follows:
Arthur Charpentier (March 25, 2020). Il était prévisible que Léonard DiCaprio ne survive pas au naufrage du Titanic. Freakonometrics. Retrieved February 9, 2025 from https://doi.org/10.58079/ovez
One thought on “Il était prévisible que Léonard DiCaprio ne survive pas au naufrage du Titanic”