Il était prévisible que Léonard DiCaprio ne survive pas au naufrage du Titanic

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 class
  • Sex: Sex of the passenger
  • Age: Age of the passenger
  • SibSp: Number of siblings/spouses aboard
  • Parch: Number of parents/children aboard
  • Embarked: Port of Embarkation
  • Name: 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}\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}\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).


Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.