Near, far,…

Lors de son voyage inaugural, de Southampton à New York via Cherbourg et Cobh (Queenstown à l’époque), le paquebot Titanic a fait naufrage en heurtant avec sa coque avant un iceberg à tribord.

Entre 1 490 et 1 520 personnes disparaissent.

Le jeu de données ‘titanic’ recense quelques informations sur 891 passagers montés à bord du paquebot.

train_data <- read.csv('titanic.csv',header=T,na.strings=c(""))
str(train_data)
## 'data.frame':    891 obs. of  12 variables:
##  $ PassengerId: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Survived   : int  0 1 1 1 0 0 0 0 1 1 ...
##  $ Pclass     : int  3 1 3 1 3 3 1 3 3 2 ...
##  $ Name       : Factor w/ 891 levels "Abbing, Mr. Anthony",..: 109 191 358 277 16 559 520 629 417 581 ...
##  $ Sex        : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...
##  $ Age        : num  22 38 26 35 35 NA 54 2 27 14 ...
##  $ SibSp      : int  1 1 0 1 0 0 0 3 0 1 ...
##  $ Parch      : int  0 0 0 0 0 0 0 1 2 0 ...
##  $ Ticket     : Factor w/ 681 levels "110152","110413",..: 524 597 670 50 473 276 86 396 345 133 ...
##  $ Fare       : num  7.25 71.28 7.92 53.1 8.05 ...
##  $ Cabin      : Factor w/ 147 levels "A10","A14","A16",..: NA 82 NA 56 NA NA 130 NA NA NA ...
##  $ Embarked   : Factor w/ 3 levels "C","Q","S": 3 1 3 3 3 2 3 3 3 1 ...

On veut expliquer “survie” : \(0\) si mort \(1\) si survivant, grâce à des covariables : “classe” \(1,2,3\), âge des passagers, etc…)

Il est aisé de compter les survivants

train_data$Age[is.na(train_data$Age)] = mean(train_data$Age, na.rm = TRUE)
table(train_data$Survived)
## 
##   0   1 
## 549 342

ou d’exprimer ce résultat en pourcentage

prop.table(table(train_data$Survived))
## 
##         0         1 
## 0.6161616 0.3838384

Sans surprise, le sexe des passagers permet de discriminer entre les survivants (les femmes…)

table(train_data$Sex,train_data$Survived)
##         
##            0   1
##   female  81 233
##   male   468 109
prop.table(table(train_data$Sex,train_data$Survived),margin=1)
##         
##                  0         1
##   female 0.2579618 0.7420382
##   male   0.8110919 0.1889081

l’âge aussi (… et les enfants d’abord!!)

hist(table(train_data$Age,train_data$Survived),nclass = 20)

wherever you are

De même, on voit une forte corrélation entre le taux de survie et la classe des passagers

table(train_data$Pclass,train_data$Survived)
##    
##       0   1
##   1  80 136
##   2  97  87
##   3 372 119

Afin d’analyser maintenant plus en détail ce jeu de données, on en retire les données non numériques (le numéro des passager, leur nom, leur numéro de cabine, etc…)

nonvars = c("PassengerId","Name","Ticket","Embarked","Cabin")
train_data = train_data[,!(names(train_data) %in% nonvars)]
str(train_data)
## 'data.frame':    891 obs. of  7 variables:
##  $ Survived: int  0 1 1 1 0 0 0 0 1 1 ...
##  $ Pclass  : int  3 1 3 1 3 3 1 3 3 2 ...
##  $ Sex     : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...
##  $ Age     : num  22 38 26 35 35 ...
##  $ SibSp   : int  1 1 0 1 0 0 0 3 0 1 ...
##  $ Parch   : int  0 0 0 0 0 0 0 1 2 0 ...
##  $ Fare    : num  7.25 71.28 7.92 53.1 8.05 ...

Les corrélations à ‘survie’ donnent un premier indice.

train_data$Sex = as.numeric(train_data$Sex); cor(train_data)
##             Survived      Pclass         Sex         Age       SibSp
## Survived  1.00000000 -0.33848104 -0.54335138 -0.06980852 -0.03532250
## Pclass   -0.33848104  1.00000000  0.13190049 -0.33133877  0.08308136
## Sex      -0.54335138  0.13190049  1.00000000  0.08415344 -0.11463081
## Age      -0.06980852 -0.33133877  0.08415344  1.00000000 -0.23262459
## SibSp    -0.03532250  0.08308136 -0.11463081 -0.23262459  1.00000000
## Parch     0.08162941  0.01844267 -0.24548896 -0.17919092  0.41483770
## Fare      0.25730652 -0.54949962 -0.18233283  0.09156609  0.15965104
##                Parch        Fare
## Survived  0.08162941  0.25730652
## Pclass    0.01844267 -0.54949962
## Sex      -0.24548896 -0.18233283
## Age      -0.17919092  0.09156609
## SibSp     0.41483770  0.15965104
## Parch     1.00000000  0.21622494
## Fare      0.21622494  1.00000000

Quelles données semblent le plus corrélées à la survie?

Régression logistique.

Il s’agit d’expliquer \(Y_i\) ‘survie’ à valeurs \(0\) ou \(1\) à partir des différentes covariables réunies dans les vecteurs \({\bf x}_i\).

Les \(Y_i\) sont des Bernoulli indépendantes, de paramètres respectifs \(\varphi(\theta^T{\bf x}_i)\).

La fonction logit \(\varphi(x)=e^{x}/(1+e^x)\) est une bijection strictement croissante de \({\mathbb{R}}\) sur \([0,1]\).

Once more, you open the door…

On a l’idée que les réponses \(0\) et \(1\) sont bien séparées par un hyperplan.

La direction de cet hyperplan est orthogonale à un vecteur \(\theta\) qu’il s’agit donc d’estimer.

Exercice 1

1.

Par indépendance des v.a. \(\{Y_k, k \geq 1 \}\), \[{\text{L}}(\theta) = \sum_{k=1}^n \ln \left( \varphi(\theta^T \mathbf{x}_k)^{Y_k} \left(1- \varphi(\theta^T \mathbf{x}_k) \right)^{1-Y_k}\right) = \sum_{k=1}^n Y_k \ln \left( \frac{\varphi(\theta^T \mathbf{x}_k)}{1-\varphi(\theta^T\mathbf{x}_k)} \right) + \sum_{k=1}^n \ln \{1-\varphi(\theta^T\mathbf{x}_k)\} \] et donc

\[{\text{L}}(\theta) = \theta^T \left( \sum_{k=1}^n Y_k \mathbf{x}_k \right) - \sum_{k=1}^n \ln\left(1 + \exp(\theta^T \mathbf{x}_k) \right)\enspace.\]

2.

Calculons le gradient: \[ \nabla {\text{L}}(\theta) = \left( \sum_{k=1}^n Y_k \mathbf{x}_k \right) - \sum_{k=1}^n \nabla \ln\left(1 + \exp(\theta^T \mathbf{x}_k) \right) = \left( \sum_{k=1}^n Y_k \mathbf{x}_k \right) - \sum_{k=1}^n \frac{\mathbf{x}_k \exp(\theta^T \mathbf{x}_k)}{1+\exp(\theta^T \mathbf{x}_k)}\enspace, \] et donc \[ \nabla {\text{L}}(\theta) = \sum_{k=1}^n \left(Y_k - \varphi(\theta^T\mathbf{x}_k) \right) \mathbf{x}_k\ . \]

I believe that my heart will go on

Calculons maintenant la hessienne \[ {\operatorname{H}_{{\text{L}}}(\theta)} = - \sum_{k=1}^n \mathbf{x}_k \left( \nabla \varphi(\theta^T \mathbf{x}_k) \right)^T = - \sum_{k=1}^n \varphi(\theta^T\mathbf{x}_k) \mathbf{x}_k \ \left( \nabla \ln \varphi(\theta^T\mathbf{x}_k) \right)^T \enspace,\] et donc

\[{\operatorname{H}_{{\text{L}}}(\theta)} = - \sum_{k=1}^n \varphi(\theta^T \mathbf{x}_k) \mathbf{x}_k \ \frac{\mathbf{x}_k^T \exp(-\theta^T \mathbf{x}_k)}{1+\exp(-\theta^T \mathbf{x}_k)} = - \sum_{k=1}^n \varphi(\theta^T \mathbf{x}_k) \{1-\varphi(\theta^T \mathbf{x}_k)\} \ \mathbf{x}_k \ \mathbf{x}_k^T\enspace. \]

Pour tout \(\lambda \in {\mathbb{R}}^p\), on a \[ \lambda^T {\operatorname{H}_{\text{L}}(\theta)} \lambda = - \sum_{k=1}^n \varphi(\theta^T \mathbf{x}_k) \{1-\varphi(\theta^T \mathbf{x}_k)\} \ \left( \lambda^T \mathbf{x}_k \right)^2 \leq 0. \] Si \(\lambda \neq 0\), alors \(\mathbf{X} \lambda \neq 0\) puisque \(\mathbf{X}\) est de rang plein en colonnes (\(\mathbf{X}\) est la matrice \(n \times p\) dont les lignes sont les vecteurs \(\mathbf{x}_1, \ldots, \mathbf{x}_n\)).

Donc il existe \(k\) tel que \(\left( \lambda^T \mathbf{x}_k \right)^2 > 0\), ainsi \(\lambda^T {\operatorname{H}_{{\text{L}}}(\theta)} \lambda< 0\).

\({\text{L}}\) est strictement concave.

3.

La fonction \(\theta \mapsto {\text{L}}(\theta)\) est strictement concave et différentiable. Donc tout point critique est un maximum global. Donc s’il existe, l’EMV est caractérisé par \[ \nabla {\text{L}}(\widehat \theta_n) = 0 \qquad \text{i.e.} \ \ \sum_{k=1}^n X_k Y_k = \sum_{k=1}^n X_k p_k(\widehat \theta_n). \]

4.

Le logiciel R utilise un algorithme appelé IWLS : dans l’aide \({\tt ??glm}\) de R, vous trouverez en particulier

\({\tt An\ object\ of\ class\ "glm" is\ a\ list\ containing\ at\ least\\ the\ following\ components\ :\\ ...\\ iter\ :\ the\ number\ of\ iterations\ of\ IWLS\ used.}\)

IWLS est un algorithme de Newton-Raphson. Rappelons que cet algorithme itératif de recherche de solution de \(\nabla{\text{L}}(\theta)=0\), partant d’une valeur \(\theta_t\), propose la mise à jour \[ \theta_{t+1}=\theta_t-({\operatorname{H}_{{\text{L}}}(\theta_t)})^{-1}\nabla{\text{L}}(\theta_t)\ . \]

Dans le cas de la régression logistique, cette mise à jour peut s’écrire sous forme matricielle, en posant

  1. \({\bf Y}\) le vecteur colonne des \(Y_i\),
  2. \({\bf X}\) la matrice \(n\times(p+1)\) des covariables,
  3. \({\bf \varphi}\) le vecteur des \(\varphi(\theta_t^T{\bf x}_i)\),
  4. \({\bf W}\) la matrice diagonale \(n\times n\) d’éléments diagonaux \(h(\theta_t^T{\bf x}_i)\),

On a \[ \nabla{\text{L}}(\theta_t)={\bf X}^T({\bf Y}-{\bf \varphi}),\quad {\operatorname{H}_{{\text{L}}}(\theta_t)}=-{\bf X}^T{\bf W}{\bf X}\ . \]

Ainsi, la mise à jour de NR se réécrit \[ \theta_{t+1}=\theta_t+({\bf X}^T{\bf W}{\bf X})^{-1}{\bf X}^T({\bf Y}-{\bf \varphi})\ . \]

En posant \({\bf z}={\bf X}\theta_t+{\bf W}^{-1}({\bf Y}-{\bf \varphi})\), on a donc \[ \theta_{t+1}=({\bf X}^T{\bf W}{\bf X})^{-1}{\bf X}^T{\bf W}{\bf z}\ . \]

\(\theta_{t+1}\) est donc solution du problème suivant \[ \theta_{t+1}=\operatorname*{argmin}_{\theta\in{\mathbb{R}}^{p+1}}({\bf z}-{\bf X}\theta)^T{\bf W}({\bf z}-{\bf X}\theta)\ . \]

Comme la mise à jour est solution d’un problème quadratique (moindres carrés à ‘poids’ \({\bf W}\)), on donne à cet algorithme le nom de IWLS pour (iteratively weighted least-squares).

Input : \(\theta=0\), \({\bf X}\), \({\bf Y}\), \(\rho>0\).

Output : EMV approché

  1. Calcul de \({\bf \varphi}\) donc les éléments valent \({\bf \varphi}_i=\varphi(\theta^T{\bf x}_i)\).
  2. Calcul de la matrice diagonale \({\bf W}\) telle que \({\bf W}_{i,i}=h(\theta^T{\bf x}_i)\).

  3. Tant que \({\left\|{\bf X}^T({\bf Y}-{\bf \varphi})\right\|}\geq \rho\) Faire
  4. \({\bf z}={\bf X}\theta+{\bf W}^{-1}({\bf Y}-{\bf \varphi})\)
  5. \(\theta=({\bf X}^T{\bf W}{\bf X})^{-1}{\bf X}^T{\bf W}{\bf z}\)
  6. Calcul de \({\bf \varphi}\) donc les éléments valent \({\bf \varphi}_i=\varphi(\theta^T{\bf x}_i)\).
  7. Calcul de la matrice diagonale \({\bf W}\) telle que \({\bf W}_{i,i}=h(\theta^T{\bf x}_i)\).
  8. Fin

  9. Retourner \(\theta\).

TitanicLog1 = glm(Survived~., data = train_data, family = binomial)
summary(TitanicLog1)
## 
## Call:
## glm(formula = Survived ~ ., family = binomial, data = train_data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7129  -0.6032  -0.4273   0.6191   2.4186  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  7.723375   0.645814  11.959  < 2e-16 ***
## Pclass      -1.084297   0.139119  -7.794 6.49e-15 ***
## Sex         -2.762930   0.199011 -13.883  < 2e-16 ***
## Age         -0.039702   0.007797  -5.092 3.55e-07 ***
## SibSp       -0.350725   0.109552  -3.201  0.00137 ** 
## Parch       -0.111963   0.117400  -0.954  0.34024    
## Fare         0.002852   0.002361   1.208  0.22718    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.66  on 890  degrees of freedom
## Residual deviance:  788.73  on 884  degrees of freedom
## AIC: 802.73
## 
## Number of Fisher Scoring iterations: 5

5.

Posons \(F(t)=\nabla {\text{L}}(t\widehat \theta_n+(1-t)\theta)\). Puisque \(\nabla {\text{L}}(\widehat \theta_n) = 0\) et \(\nabla {\text{L}}\) est dérivable sur \({\mathbb{R}}^p\), \(F(1)-F(0)=\int_0^1F'(t)dt\), c’est à dire \[ - \nabla {\text{L}}(\theta) = {\left[\int_0^1{\operatorname{H}_{{\text{L}}}(t\widehat\theta_{n}+(1-t)\theta)}dt\right]} (\widehat\theta_n-\theta) \enspace.\]

On en déduit \[ - \frac{\nabla {\text{L}}(\theta)}{\sqrt{n}} = \left(\frac{{\operatorname{H}_{{\text{L}}}(\theta)}}n + R_n \right) \sqrt{n} (\widehat\theta_n-\theta)\enspace, \] avec \[ R_n=\int_0^1\frac{{\operatorname{H}_{{\text{L}}}(t\widehat\theta_{n}+(1-t)\theta)}-{\operatorname{H}_{{\text{L}}}(\theta)}}ndt \]

Pour conclure la preuve, on montre que

\[ \forall n\ge 1,\qquad n^{-1}{\operatorname{H}_{{\text{L}}}(\cdot)} \text{ est $C$-Lipshitzienne} \ . \]

Soient \(\widetilde{\theta}\) et \(\theta\) dans \({\mathbb{R}}^p\), on a \[ n^{-1} \left\| {\operatorname{H}_{{\text{L}}}(\widetilde \theta)} -{\operatorname{H}_{{\text{L}}}(\theta)} \right\| = \frac{1}{n} \left\| \sum_{k=1}^n \left\{ h\left( \widetilde \theta^T {\mathbf{x}}_k \right) - h\left( \theta^T {\mathbf{x}}_k \right) \right\} {\mathbf{x}}_k {\mathbf{x}}_k^T \right\| \ . \]

Or, \(h\) est \(1\)-Lipschitzienne : \[ h(t) = \frac{\exp(t)}{(1+ \exp(t))^2},\quad \text{donc}\quad h'(t) = \frac{\exp(t)}{1 +\exp(t)} \frac{1-\exp(t)}{(1+\exp(t))^2}\ . \] Donc \(h'\) est majorée par \(1\). \[ n^{-1} \left\| {\operatorname{H}_{{\text{L}}}(\widetilde \theta)} -{\operatorname{H}_{{\text{L}}}(\theta)} \right\| \leq \|\widetilde \theta - \theta \| \frac{1}{n} \sum_{k=1}^n \|{\mathbf{x}}_k\| \|{\mathbf{x}}_k {\mathbf{x}}_k^T \|\ . \] Par équivalence des normes sur \(\mathcal{M}_p({\mathbb{R}})\), on a de plus, l’existence d’une constante \(C_p\) telle que \[ \|{\mathbf{x}}_k {\mathbf{x}}_k^T \|\le C_p\sqrt{\sum_{i,j=1}^n({\mathbf{x}}_k {\mathbf{x}}_k^T)^2_{i,j}}\le C_p{\left\|{\mathbf{x}}_k\right\|}^2\ . \] Ainsi \[ n^{-1} \left\| {\operatorname{H}_{{\text{L}}}(\widetilde \theta)} -{\operatorname{H}_{{\text{L}}}(\theta)} \right\| \leq C_p\|\widetilde \theta - \theta \|\sup_{n\ge 1} \frac{1}{n} \sum_{k=1}^n \|{\mathbf{x}}_k\|^3\ . \]

Puisque \(\sup_n n^{-1} \sum_{k=1}^n \|{\mathbf{x}}_k\|^3 < \infty\), les fonctions \(n^{-1}{\operatorname{H}_{{\text{L}}}(\cdot)}\) sont bien toutes \(C\)-Lipshitziennes.

On a donc, \[ R_n\le C{\left\|\widehat \theta_n - \theta\right\|}\int_0^1tdt\le C \|\widehat \theta_n - \theta \|\ . \] On conclut avec l’hypothèse \(\widehat \theta_n {\stackrel{{\mathbb{P}}_\theta-\operatorname{prob}}{\longrightarrow}} \theta\).

6.

Vu le calcul du gradient fait précédemment, on a \[ \frac {\nabla {\text{L}}(\theta)}{\sqrt n} = \sqrt{n} \left( \frac{1}{n} \sum_{k=1}^n \left(Y_k - \varphi(\theta^T {\mathbf{x}}_k) \right) {\mathbf{x}}_k \right)\enspace. \] Pour établir la normalité asymptotique, on applique le théorème de Lindeberg-Feller (cf. Théorème III-2.41 du polycopié) avec \[ Z_{n,k} = \frac{1}{\sqrt{n}} Y_k {\mathbf{x}}_k, \qquad \mu_{n,k} = \frac{\varphi(\theta^T {\mathbf{x}}_k)}{\sqrt{n}} {\mathbf{x}}_k\ ,\] \[ \Sigma_{n,k} = \frac{\varphi(\theta^T {\mathbf{x}}_k) (1-\varphi(\theta^T {\mathbf{x}}_k))}{n} {\mathbf{x}}_k {\mathbf{x}}_k^T\ . \]

Rappel : Théorème de Lindeberg-Feller

Supposons que, \(Z_{n,1}, \ldots, Z_{n,n}\) sont indépendantes, d’espérances \(\mu_{n,1}, \ldots, \mu_{n,n}\) et de matrices de covariance \(\Sigma_{n,1}, \ldots, \Sigma_{n,n}\). Si \[ \forall\epsilon>0, \quad \lim_n \sum_{k=1}^n {\mathbb{E}}\left[ \| Z_{n,k}\|^2 {\mathbf{1}}_{\|Z_{n,k}\|>\epsilon} \right] = 0 \ ,\] \[ \lim_n \sum_{k=1}^n \Sigma_{n,k} = \Sigma \ , \] alors \[ \sum_{k=1}^n \left( Z_{n,k} - \mu_{n,k} \right) {\stackrel{}{\Longrightarrow}} {\operatorname{N}}(0,\Sigma)\ . \]

On a, par hypothèses, \[ \sum_{k=1}^n \Sigma_{n,k} = \frac{1}{n}\sum_{k=1}^n \varphi(\theta^T {\mathbf{x}}_k) (1-\varphi(\theta^T {\mathbf{x}}_k)) \ {\mathbf{x}}_k {\mathbf{x}}_k^T \to Q_\theta\ . \] La condition 2 du théorème est donc bien satisfaite.

\[ \sum_{k=1}^n {\mathbb{E}}_\theta\left[ \| Z_{n,k}\|^2 {\mathbf{1}}_{\|Z_{n,k}\|>\epsilon} \right] = \frac{1}{n} \sum_{k=1}^n \|{\mathbf{x}}_k\|^2 {\mathbb{E}}_\theta\left[ Y_k^2 {\mathbf{1}}_{Y_k \|{\mathbf{x}}_k\|>\epsilon \sqrt{n}} \right] \enspace.\] Comme \({\mathbf{1}}_{Y_k \|{\mathbf{x}}_k\|>\epsilon\sqrt{n}} \leq |Y_k| \|{\mathbf{x}}_k\|/\epsilon\sqrt{n}\),

\[ \sum_{k=1}^n {\mathbb{E}}_\theta\left[ \| Z_{n,k}\|^2 {\mathbf{1}}_{\|Z_{n,k}\|>\epsilon} \right] \leq\frac{1}{\epsilon\sqrt{n}}\frac1n \sum_{k=1}^n \|{\mathbf{x}}_k\|^3 {\mathbb{E}}_\theta\left[ Y_k^3\right]\ .\] Donc, comme \(C:=\sup_n n^{-1} \sum_{k=1}^n \|{\mathbf{x}}_k\|^3 < \infty\) et \({\mathbb{E}}[Y_k^3]\le 1\), \[ \sum_{k=1}^n {\mathbb{E}}_\theta\left[ \| Z_{n,k}\|^2 {\mathbf{1}}_{\|Z_{n,k}\|>\epsilon} \right]\leq \frac{C}{\epsilon \sqrt{n}}\ . \] Ceci établit la condition 1. Par Lindeberg Feller, on conclut \[ \frac {\nabla {\text{L}}(\theta)}{\sqrt n} {\stackrel{{\mathbb{P}}_\theta}{\Longrightarrow}} {\operatorname{N}}(0, Q_\theta). \]

7.

On déduit des questions 5 et 6 que \[ \left(\frac {{\operatorname{H}_{{\text{L}}}(\theta)}}{n}+ R_n\right)\, \sqrt n(\widehat\theta_n-\theta) {\stackrel{{\mathbb{P}}_\theta}{\Longrightarrow}} {\operatorname{N}}(0, Q_\theta)\ . \] Par ailleurs \[ \frac{{\operatorname{H}_{{\text{L}}}(\theta)}}{n}+ R_n = - \frac{1}{n} \sum_{k=1}^n h(\theta^T {\mathbf{x}}_k) {\mathbf{x}}_k {\mathbf{x}}_k^T + R_n {\stackrel{{\mathbb{P}}_\theta-\operatorname{prob}}{\longrightarrow}} - Q_\theta\ . \] Donc, d’après le lemme de Slutsky, comme \(Q_\theta\) est symétrique et inversible, \[ \sqrt n(\widehat\theta_n-\theta) {\stackrel{{\mathbb{P}}_\theta}{\Longrightarrow}} {\operatorname{N}}(0, Q_\theta^{-1})\ . \]

8.

D’après la question précédente, \[ \frac{1}{\sqrt{(Q_\theta^{-1})_{k,k}}}\sqrt n(\widehat\theta_{n,k}-\theta_k) {\stackrel{{\mathbb{P}}_\theta}{\Longrightarrow}} {\operatorname{N}}(0, 1)\ . \]

Comme \(n^{-1}{\operatorname{H}_{{\text{L}}}(\cdot)}\) est Lipshitzienne et \(\widehat \theta_n {\stackrel{{\mathbb{P}}_\theta-\operatorname{prob}}{\longrightarrow}} \theta\), on a \[ {\left\|\frac{{\operatorname{H}_{{\text{L}}}(\widehat\theta_n)}}{n} - \frac{{\operatorname{H}_{{\text{L}}}(\theta)}}{n}\right\|}\le C{\left\|\widehat\theta_n-\theta\right\|} {\stackrel{{\mathbb{P}}_\theta-\operatorname{prob}}{\longrightarrow}} 0\ . \]

Puis par définition de \(Q_\theta\), \(\lim_n - \frac{{\operatorname{H}_{{\text{L}}}(\theta)}}{n} = Q_\theta\). Par suite, \[ - \frac{{\operatorname{H}_{{\text{L}}}(\widehat\theta_n)}}{n} {\stackrel{{\mathbb{P}}_\theta-\operatorname{prob}}{\longrightarrow}} Q_\theta\ . \] En utilisant le théorème de l’application continue pour la convergence en probabilité, on en déduit \[ \beta_{n,k} {\stackrel{{\mathbb{P}}_\theta-\operatorname{prob}}{\longrightarrow}} (Q_\theta^{-1})_{k,k}\ . \] On conclut en utilisant le Lemme de Slutsky.

9.

Soit \(T={\mathbf{1}}_{(Y_1,\ldots,Y_n)\in \mathcal{R}}\) le test pur de zone de rejet \[ \mathcal{R}=\left\{(y_1,\ldots,y_n) : \sqrt{\frac{n}{\beta_{n,k}}}\ \left| \widehat\theta_{n,k} \right| > z_{1-\alpha/2} \right\} \]\(z_q\) est le quantile d’ordre \(q\) de la loi \({\operatorname{N}}(0,1)\). On a bien d’après la question que, pour tout \(\theta\in{\mathbb{R}}^p\) tel que \(\theta_k=0\), \[ {\mathbb{P}}_\theta{\left(\sqrt{\frac{n}{\beta_{n,k}}}\ \left| \widehat\theta_{n,k} \right| > z_{1-\alpha/2}\right)}\to {\mathbb{P}}{\left(|{\operatorname{N}}(0,1)|> z_{1-\alpha/2}\right)}=\alpha\ . \] Le test est bien de niveau asymptotique \(\alpha\).

10.

Etant donnée une réalisation de \(\sqrt{\frac{n}{\beta_{n,k}}}\ \left| \widehat\theta_{n,k} \right|\) la \(p\)-valeur asymptotique du test précédent est \[ {\mathbb{P}}{\left(|{\operatorname{N}}(0,1)|>\sqrt{\frac{n}{\beta_{n,k}}}\ \left| \widehat\theta_{n,k} \right|\ |\ Y_1,\ldots,Y_n\right)}\ . \]

TitanicLog1 = glm(Survived~., data = train_data, family = binomial)
summary(TitanicLog1)
## 
## Call:
## glm(formula = Survived ~ ., family = binomial, data = train_data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7129  -0.6032  -0.4273   0.6191   2.4186  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  7.723375   0.645814  11.959  < 2e-16 ***
## Pclass      -1.084297   0.139119  -7.794 6.49e-15 ***
## Sex         -2.762930   0.199011 -13.883  < 2e-16 ***
## Age         -0.039702   0.007797  -5.092 3.55e-07 ***
## SibSp       -0.350725   0.109552  -3.201  0.00137 ** 
## Parch       -0.111963   0.117400  -0.954  0.34024    
## Fare         0.002852   0.002361   1.208  0.22718    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.66  on 890  degrees of freedom
## Residual deviance:  788.73  on 884  degrees of freedom
## AIC: 802.73
## 
## Number of Fisher Scoring iterations: 5

Every night in my dreams

Expliquer les dernières colonnes.

Commentez le résultat en regard de l’analyse des corrélations faite au début, commentez en particulier les conclusions potentielles dans les 2 cas sur les variables ‘Age’, ‘SibSp’ et ‘Fare’.

Considérons un modèle restreint.

TitanicLog2 = glm(Survived ~ . - Parch, data = train_data, family = binomial)
summary(TitanicLog2)
## 
## Call:
## glm(formula = Survived ~ . - Parch, family = binomial, data = train_data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7458  -0.5948  -0.4170   0.6109   2.4501  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  7.668775   0.641018  11.963  < 2e-16 ***
## Pclass      -1.098189   0.137969  -7.960 1.72e-15 ***
## Sex         -2.726408   0.194561 -14.013  < 2e-16 ***
## Age         -0.039385   0.007773  -5.067 4.05e-07 ***
## SibSp       -0.378646   0.106212  -3.565 0.000364 ***
## Fare         0.002373   0.002250   1.054 0.291707    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.66  on 890  degrees of freedom
## Residual deviance:  789.65  on 885  degrees of freedom
## AIC: 801.65
## 
## Number of Fisher Scoring iterations: 5

we’ll stay forever this way

Quelle est la différence avec le modèle précédent?

Commenter en particulier : le choix de garder la variable ‘Fare’, l’évolution de la significativité de ‘SibSp’.

Quel modèle proposez-vous de tester ensuite?

TitanicLog3 = glm(Survived ~ . - Parch - Fare, data = train_data, family = binomial)
summary(TitanicLog3)
## 
## Call:
## glm(formula = Survived ~ . - Parch - Fare, family = binomial, 
##     data = train_data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6869  -0.6055  -0.4169   0.6111   2.4547  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  7.931782   0.594109  13.351  < 2e-16 ***
## Pclass      -1.172391   0.119725  -9.792  < 2e-16 ***
## Sex         -2.739806   0.194142 -14.112  < 2e-16 ***
## Age         -0.039793   0.007755  -5.131 2.88e-07 ***
## SibSp       -0.357788   0.104033  -3.439 0.000583 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.66  on 890  degrees of freedom
## Residual deviance:  790.84  on 886  degrees of freedom
## AIC: 800.84
## 
## Number of Fisher Scoring iterations: 5