Corrigé PC8
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 342ou d’exprimer ce résultat en pourcentage
prop.table(table(train_data$Survived))##
## 0 1
## 0.6161616 0.3838384Sans 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 109prop.table(table(train_data$Sex,train_data$Survived),margin=1)##
## 0 1
## female 0.2579618 0.7420382
## male 0.8110919 0.1889081l’â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 119Afin 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.00000000Quelles 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
- \({\bf Y}\) le vecteur colonne des \(Y_i\),
- \({\bf X}\) la matrice \(n\times(p+1)\) des covariables,
- \({\bf \varphi}\) le vecteur des \(\varphi(\theta_t^T{\bf x}_i)\),
- \({\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é
- Calcul de \({\bf \varphi}\) donc les éléments valent \({\bf \varphi}_i=\varphi(\theta^T{\bf x}_i)\).
Calcul de la matrice diagonale \({\bf W}\) telle que \({\bf W}_{i,i}=h(\theta^T{\bf x}_i)\).
- Tant que \({\left\|{\bf X}^T({\bf Y}-{\bf \varphi})\right\|}\geq \rho\) Faire
- \({\bf z}={\bf X}\theta+{\bf W}^{-1}({\bf Y}-{\bf \varphi})\)
- \(\theta=({\bf X}^T{\bf W}{\bf X})^{-1}{\bf X}^T{\bf W}{\bf z}\)
- Calcul de \({\bf \varphi}\) donc les éléments valent \({\bf \varphi}_i=\varphi(\theta^T{\bf x}_i)\).
- Calcul de la matrice diagonale \({\bf W}\) telle que \({\bf W}_{i,i}=h(\theta^T{\bf x}_i)\).
Fin
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: 55.
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\} \] où \(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: 5Every 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: 5we’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