{"word_count": 1470, "character_count": 13826}
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 passagersloc_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 passengerLa 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 estmean(base$Survived)[1] 0.3838384que l'on peut retrouver dans le tableau de contingence suivanttable(base$Survived)/nrow(base)        0         1 0.6161616 0.3838384soit [latex]38.38\%[/latex] 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 iciglm(Survived~1, family=binomial,data=base)Coefficients:(Intercept)      -0.4733  Degrees of Freedom: 890 Total (i.e. Null);  890 ResidualNull Deviance:    1187 Residual Deviance: 1187 AIC: 1189ce qui donne la valeur de [latex]\widehat{\beta}_0[/latex], et comme la probabilité de survie  est [latex display="true"]\mathbb{P}[Y=1]=\frac{\exp(\beta_0)}{1+\exp(\beta_0)}[/latex]on en déduit une prévision de cette probabilité en considérant[latex display="true"]\frac{\exp(\widehat{\beta}_0)}{1+\exp(\widehat{\beta}_0)}[/latex]soit iciexp(-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 [latex]\widehat{\beta}_0[/latex])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.3838384Si on note [latex display="true"]\widehat{\mu}_i=\frac{\exp(\widehat{\beta}_0)}{1+\exp(\widehat{\beta}_0)}[/latex], on retrouve ici le fait que pour une régression logistique, [latex display="true"]\frac{1}{n}\sum_{i-1}^n\widehat{\mu}_i=\frac{1}{n}\sum_{i-1}^n y_i[/latex]Cette 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 [latex]\boldsymbol{x}[/latex],[latex display="true"]\log\mathcal{L}=\sum_{i=1}^ny_i \log p(\boldsymbol{x}_i)+(1-y_i) \log(1-p(\boldsymbol{x}_i))[/latex]soit[latex display="true"]\log\mathcal{L}=\sum_{i=1}^ny_i[\boldsymbol{x}_i^\top\boldsymbol{\beta}]-\log(1+\exp[\boldsymbol{x}_i^\top\boldsymbol{\beta}])[/latex]et la condition du premier ordre donne ici[latex display="true"]\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}[/latex]où [latex display="true"]\widehat{\mu}_i=\frac{\exp(\boldsymbol{x}_i^\top\widehat{\boldsymbol{\beta}})}{1+\exp(\boldsymbol{x}_i^\top\widehat{\boldsymbol{\beta}})}[/latex]donc dans le cas particulier de la première condition du premier ordre, on obtient [latex display="true"]\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][/latex]et donc[latex display="true"]\frac{1}{n}\sum_{i-1}^n\widehat{\mu}_i=\frac{1}{n}\sum_{i-1}^n y_i[/latex]D'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)&#62; 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.6296296Environ [latex]63\%[/latex] ont survécu. Mais là encore, on peut faire une régression logistiqueglm(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 ResidualNull Deviance:    1187 Residual Deviance: 1083 AIC: 1089Comme la classe 1 est la classe de référence,  on fait comme auparavant, en considérant[latex display="true"]\frac{\exp(\widehat{\beta}_0)}{1+\exp(\widehat{\beta}_0)}[/latex]soit iciexp(0.5306)/(1+exp(0.5306))[1] 0.629623(avec toujours la petite imprécision numérique). Une prévision plus propre donne icireg=glm(Survived~Pclass, family=binomial,data=base)predict(reg,newdata=data.frame(Pclass="1"),type="response")        1 0.6296296Et 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 suivantechisq.test(table(base$Survived,base$Pclass))Pearson's Chi-squared testdata:  table(base$Survived, base$Pclass)X-squared = 102.89, df = 2, p-value &#60; 2.2e-16Juste un petit rappel (on a revu le test du chi-deux au dernier cours), on a un tableau de contingence, avec [latex]n_{i,j}[/latex] de nombre de personnes de type [latex]i[/latex] (pour la première variable) et [latex]j[/latex] (pour la seconde). Si les variables étaient indépendantes, on aura [latex display="true"]n_{i,j}^\perp=n\frac{n_{i\cdot}}{n}\frac{n_{\cdot j}}{n}[/latex]où[latex display="true"]\sum_{i=1}^n n_{i,j}=n_{\cdot j}~~~\sum_{j=1}^n n_{i,j}=n_{i\cdot}[/latex]La statistique de test est alors [latex display="true"]\displaystyle{\chi ^{2}=\sum _{i=1}^{I}\sum _{j=1}^{J}\frac{(n_{i,j}-n^\perp_{i,j})^{2}}{n^\perp_{i,j}}}[/latex]On peut noter que les contributions au test sont ici[latex display="true"]\varepsilon_{i,j}=\frac{n_{i,j}-n^\perp_{i,j}}{\sqrt{n^\perp_{i,j}}}[/latex]chisq.test(table(base$Survived,base$Pclass))$residuals               1         2         3  0 -4.601993 -1.537771  3.993703  1  5.830678  1.948340 -5.059981ce 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-dessouslibrary("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 paquebotOn nous donne la base avec deux passagersnewbase = 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(&#62;|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  &#60; 2e-16 ***Sexmale      -2.641145   0.222801 -11.854  &#60; 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 freedomResidual deviance: 632.67  on 704  degrees of freedom  (177 observations deleted due to missingness)AIC: 652.67Number of Fisher Scoring iterations: 13Deux variables ne sont pas significatives. On va les enlever, une à une, et on arrive au modèle suivantreg=glm(Survived~Pclass+Sex+Age+SibSp, family=binomial,data=base)summary(reg)Coefficients:             Estimate Std. Error z value Pr(&#62;|z|)    (Intercept)  4.334201   0.450700   9.617  &#60; 2e-16 ***Pclass2     -1.414360   0.284727  -4.967 6.78e-07 ***Pclass3     -2.652618   0.285832  -9.280  &#60; 2e-16 ***Sexmale     -2.627679   0.214771 -12.235  &#60; 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 freedomResidual deviance: 636.56  on 708  degrees of freedom  (177 observations deleted due to missingness)AIC: 648.56Number of Fisher Scoring iterations: 5Maintenant, 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 polynomialereg3=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(&#62;|z|)    (Intercept)     3.0213     0.2903  10.408  &#60; 2e-16 ***Pclass2        -1.3603     0.2842  -4.786 1.70e-06 ***Pclass3        -2.5569     0.2853  -8.962  &#60; 2e-16 ***Sexmale        -2.6582     0.2176 -12.216  &#60; 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 freedomResidual deviance: 627.55  on 706  degrees of freedomAIC: 643.55Number of Fisher Scoring iterations: 5mais 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 mainreg3=glm(Survived~Pclass+Sex+Age+I(Age^2)+I(Age^3)+SibSp, family=binomial,data=base[!is.na(base$Age),])&#62; summary(reg3)Coefficients:              Estimate Std. Error z value Pr(&#62;|z|)    (Intercept)  5.616e+00  6.565e-01   8.554  &#60; 2e-16 ***Pclass2     -1.360e+00  2.842e-01  -4.786  1.7e-06 ***Pclass3     -2.557e+00  2.853e-01  -8.962  &#60; 2e-16 ***Sexmale     -2.658e+00  2.176e-01 -12.216  &#60; 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 freedomResidual deviance: 627.55  on 706  degrees of freedomAIC: 643.55Number of Fisher Scoring iterations: 5On notera (comme je le disais dans Function basis and regression) que la [latex]p[/latex]-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 fonctionbeta = coefficients(reg3)xage = 0:80yage = beta[5]*xage+beta[6]*xage^2+beta[7]*xage^3plot(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 visualisationlibrary(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 freedomResidual Deviance: 627.5525 on 706 degrees of freedomAIC: 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(&#62;F)    Pclass      2  26.72  13.361  11.3500 1.407e-05 ***Sex         1 131.57 131.573 111.7678 &#60; 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 ‘ ’ 1On 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éaireplot(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.1218426Autrement 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).