Exercices théoriques

Cette page rassemble des exercices théoriques pour chaque chapitre des notes de cours. Les exercices sont classés en trois niveaux de difficulté : facile, moyen et difficile. La solution de chaque exercice est repliée par défaut ; cliquez sur « ✅ Solution » pour l’afficher.

Accès rapide :

• Chapitre 1 : Rappels de régression linéaire
• Chapitre 3 : Modèles de régression pour données binaires
• Chapitre 4 : Modèles linéaires généralisés (GLM)

Chapitre 1 : Rappels de régression linéaire

Niveau facile

Exercice 1. Soit le modèle de régression simple \(Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\), pour \(i=1,\ldots,n\).

1. Écrivez la matrice de design \(\boldsymbol{X}\), puis calculez explicitement \(\boldsymbol{X}^\top\boldsymbol{X}\) et \(\boldsymbol{X}^\top\boldsymbol{Y}\) en fonction de \(\bar x\), \(\bar Y\), \(\sum x_i^2\) et \(\sum x_i Y_i\).

2. À partir de \(\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{Y}\), retrouvez les formules usuelles \[\hat\beta_1 = \frac{\sum_i (x_i-\bar x)(Y_i - \bar Y)}{\sum_i (x_i - \bar x)^2} \qquad\text{et}\qquad \hat\beta_0 = \bar Y - \hat\beta_1 \bar x.\]

Solution

On a \(\boldsymbol{X}_i^\top = (1, x_i)\), donc

\[\boldsymbol{X}^\top\boldsymbol{X} = \begin{pmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2 \end{pmatrix}, \qquad \boldsymbol{X}^\top\boldsymbol{Y} = \begin{pmatrix} \sum Y_i \\ \sum x_i Y_i \end{pmatrix}.\]

Le déterminant de \(\boldsymbol{X}^\top\boldsymbol{X}\) est \(n\sum x_i^2 - (\sum x_i)^2 = n\sum_i(x_i-\bar x)^2\), de sorte que

\[(\boldsymbol{X}^\top\boldsymbol{X})^{-1} = \frac{1}{n\sum_i(x_i-\bar x)^2}\begin{pmatrix} \sum x_i^2 & -\sum x_i \\ -\sum x_i & n \end{pmatrix}.\]

En multipliant par \(\boldsymbol{X}^\top\boldsymbol{Y}\), la seconde composante donne

\[\hat\beta_1 = \frac{n\sum x_i Y_i - \sum x_i \sum Y_i}{n\sum_i(x_i-\bar x)^2} = \frac{\sum_i(x_i-\bar x)(Y_i-\bar Y)}{\sum_i(x_i-\bar x)^2},\]

en utilisant l’identité \(n\sum x_iY_i - \sum x_i\sum Y_i = n\sum_i(x_i-\bar x)(Y_i-\bar Y)\). La première composante donne, après simplification, \(\hat\beta_0 = \bar Y - \hat\beta_1\bar x\).

Exercice 2. Lorsque le modèle contient une ordonnée à l’origine, les équations normales s’écrivent \[\boldsymbol{X}^\top(\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}}) = \boldsymbol{0}.\] Montrez qu’elles impliquent que \(\sum_{i=1}^n e_i = 0\) et que \(\sum_{i=1}^n x_{ij} e_i = 0\) pour chaque variable explicative \(j\), où \(e_i = Y_i - \hat Y_i\). Que peut-on en déduire sur la moyenne des valeurs prédites \(\hat Y_i\) ?

Solution

En notant \(\boldsymbol{e}=\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}}\), les équations normales s’écrivent \(\boldsymbol{X}^\top\boldsymbol{e}=\boldsymbol{0}\). La première colonne de \(\boldsymbol{X}\) ne contenant que des \(1\), la première ligne de ce système est exactement \(\sum_{i=1}^n e_i = 0\). Chaque autre ligne, associée à la colonne des \(x_{ij}\), donne \(\sum_{i=1}^n x_{ij} e_i = 0\).

Puisque \(\sum_i e_i = \sum_i (Y_i - \hat Y_i) = 0\), on a \(\sum_i Y_i = \sum_i \hat Y_i\), donc la moyenne des valeurs prédites est égale à la moyenne observée : \(\bar{\hat Y} = \bar Y\).

Exercice 3. On ajoute une variable explicative supplémentaire à un modèle de régression linéaire multiple. Montrez algébriquement (à partir de la définition de \(SSE\) comme minimum d’une somme de carrés) que le \(R^2\) ne peut pas diminuer. Expliquez ensuite, sans calcul, pourquoi le \(R^2\) ajusté peut, lui, diminuer.

Solution

Le modèle à \(p\) variables est un cas particulier du modèle à \(p+1\) variables (il suffit de fixer à \(0\) le coefficient de la variable ajoutée). Comme les moindres carrés minimisent \(SSE\) sur l’ensemble de tous les coefficients possibles, et que cet ensemble est plus grand pour le modèle à \(p+1\) variables (il contient en particulier la solution du modèle à \(p\) variables), on a nécessairement \(SSE_{p+1} \leq SSE_p\). Puisque \(R^2 = 1-SSE/SST\) et que \(SST\) ne dépend que de \(Y\) (donc reste inchangé), \(R^2\) ne peut qu’augmenter ou rester constant.

Le \(R^2\) ajusté, \[R^2_{adj} = 1 - \frac{SSE/(n-(p+1))}{SST/(n-1)},\] divise \(SSE\) par les degrés de liberté résiduels \(n-(p+1)\), qui diminuent quand on ajoute une variable. Si la diminution de \(SSE\) est trop faible pour compenser la diminution des degrés de liberté, le ratio \(SSE/(n-(p+1))\) augmente, et \(R^2_{adj}\) diminue : ajouter une variable peu informative est donc pénalisé.

Niveau moyen

Exercice 1. Soit \(\boldsymbol{H} = \boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\) la matrice chapeau. Montrez que \(\boldsymbol{H}\) est symétrique et idempotente (\(\boldsymbol{H}^2=\boldsymbol{H}\)), puis que \(\boldsymbol{I}-\boldsymbol{H}\) l’est également. En déduire que \(Var(\boldsymbol{e}) = \sigma^2(\boldsymbol{I}-\boldsymbol{H})\), où \(\boldsymbol{e} = \boldsymbol{Y}-\hat{\boldsymbol{Y}}\), et que \(Var(e_i) = \sigma^2(1-h_i)\).

Solution

Symétrie : \(\boldsymbol{H}^\top = \boldsymbol{X}\{(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\}^\top\boldsymbol{X}^\top = \boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top = \boldsymbol{H}\), car \((\boldsymbol{X}^\top\boldsymbol{X})^{-1}\) est symétrique.

Idempotence : \(\boldsymbol{H}^2 = \boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top = \boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top = \boldsymbol{H}\).

Pour \(\boldsymbol{I}-\boldsymbol{H}\) : la symétrie est immédiate, et \((\boldsymbol{I}-\boldsymbol{H})^2 = \boldsymbol{I}-2\boldsymbol{H}+\boldsymbol{H}^2 = \boldsymbol{I}-2\boldsymbol{H}+\boldsymbol{H}=\boldsymbol{I}-\boldsymbol{H}\).

Puisque \(\boldsymbol{e}=(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{Y}\) et \(Var(\boldsymbol{Y})=\sigma^2\boldsymbol{I}\),

\[Var(\boldsymbol{e}) = (\boldsymbol{I}-\boldsymbol{H})\,\sigma^2\boldsymbol{I}\,(\boldsymbol{I}-\boldsymbol{H})^\top = \sigma^2(\boldsymbol{I}-\boldsymbol{H})^2 = \sigma^2(\boldsymbol{I}-\boldsymbol{H}).\]

L’élément diagonal \(i\) de cette matrice donne \(Var(e_i) = \sigma^2(1-h_i)\).

Exercice 2. En supposant \(\boldsymbol{\epsilon} \sim \mathcal{N}(\boldsymbol{0}, \sigma^2\boldsymbol{I})\), écrivez la vraisemblance du modèle linéaire et dérivez les estimateurs du maximum de vraisemblance \(\hat{\boldsymbol{\beta}}_{MV}\) et \(\hat\sigma^2_{MV}\). Montrez que \(\hat{\boldsymbol{\beta}}_{MV}\) coïncide avec l’estimateur des moindres carrés, mais que \(\hat\sigma^2_{MV}\) est biaisé, et précisez le lien entre \(\hat\sigma^2_{MV}\) et l’estimateur sans biais \(\hat\sigma^2\).

Solution

La log-vraisemblance est

\[l(\boldsymbol{\beta},\sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}(\boldsymbol{Y}-\boldsymbol{X\beta})^\top(\boldsymbol{Y}-\boldsymbol{X\beta}).\]

Pour \(\sigma^2\) fixé, maximiser \(l\) en \(\boldsymbol{\beta}\) revient à minimiser \(SS_E(\boldsymbol{\beta})=(\boldsymbol{Y}-\boldsymbol{X\beta})^\top(\boldsymbol{Y}-\boldsymbol{X\beta})\), donc \(\hat{\boldsymbol{\beta}}_{MV} = \hat{\boldsymbol{\beta}}\), l’estimateur des moindres carrés.

En dérivant par rapport à \(\sigma^2\) et en annulant :

\[\frac{\partial l}{\partial\sigma^2} = -\frac{n}{2\sigma^2} + \frac{SS_E}{2\sigma^4} = 0 \quad\Longrightarrow\quad \hat\sigma^2_{MV} = \frac{SS_E}{n}.\]

Comme \(\hat\sigma^2 = SS_E/(n-(p+1))\), on a \[\hat\sigma^2_{MV} = \frac{n-(p+1)}{n}\,\hat\sigma^2.\] Puisque \(E[\hat\sigma^2]=\sigma^2\) (exercice suivant), \(E[\hat\sigma^2_{MV}] = \frac{n-(p+1)}{n}\sigma^2 < \sigma^2\) : l’estimateur du maximum de vraisemblance sous-estime systématiquement \(\sigma^2\).

Exercice 3. Montrez que \(E[SS_E] = \sigma^2 (n-(p+1))\), ce qui justifie que \(\hat\sigma^2 = SS_E/(n-(p+1))\) est un estimateur sans biais de \(\sigma^2\). (Indice : utilisez le fait que \(SS_E = \boldsymbol{Y}^\top(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{Y}\) et la formule de l’espérance d’une forme quadratique.)

Solution

On a \(SS_E = \boldsymbol{e}^\top\boldsymbol{e} = \boldsymbol{Y}^\top(\boldsymbol{I}-\boldsymbol{H})^\top(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{Y} = \boldsymbol{Y}^\top(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{Y}\), car \(\boldsymbol{I}-\boldsymbol{H}\) est symétrique et idempotente. Pour une forme quadratique \(\boldsymbol{Y}^\top\boldsymbol{A}\boldsymbol{Y}\) avec \(E[\boldsymbol{Y}]=\boldsymbol{\mu}\) et \(Var(\boldsymbol{Y})=\sigma^2\boldsymbol{I}\), on a \(E[\boldsymbol{Y}^\top\boldsymbol{A}\boldsymbol{Y}] = \sigma^2\,\text{tr}(\boldsymbol{A}) + \boldsymbol{\mu}^\top\boldsymbol{A}\boldsymbol{\mu}\).

Ici \(\boldsymbol{\mu}=\boldsymbol{X\beta}\) et \((\boldsymbol{I}-\boldsymbol{H})\boldsymbol{X}=\boldsymbol{X}-\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{X} = \boldsymbol{X}-\boldsymbol{X}=\boldsymbol{0}\), donc le second terme s’annule. Il reste

\[E[SS_E] = \sigma^2\,\text{tr}(\boldsymbol{I}-\boldsymbol{H}) = \sigma^2\{n - \text{tr}(\boldsymbol{H})\}.\]

Or \(\text{tr}(\boldsymbol{H}) = \text{tr}\{\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\} = \text{tr}\{(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{X}\} = \text{tr}(\boldsymbol{I}_{p+1}) = p+1\) (invariance cyclique de la trace). Donc \(E[SS_E] = \sigma^2(n-(p+1))\), et \(\hat\sigma^2\) est sans biais.

Niveau difficile

Exercice 1. Démontrez le théorème de Gauss-Markov : montrez que, parmi tous les estimateurs linéaires sans biais de \(\boldsymbol{\beta}\), l’estimateur des moindres carrés \(\hat{\boldsymbol{\beta}}\) est celui de variance minimale (BLUE, best linear unbiased estimator). (Indice : considérez un estimateur linéaire quelconque \(\tilde{\boldsymbol{\beta}} = \boldsymbol{CY}\) sans biais, écrivez \(\boldsymbol{C} = (\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top + \boldsymbol{D}\) et montrez que \(\boldsymbol{DX}=\boldsymbol{0}\).)

Solution

Soit \(\tilde{\boldsymbol{\beta}}=\boldsymbol{CY}\) un estimateur linéaire quelconque. Le biais est \(E[\tilde{\boldsymbol{\beta}}]-\boldsymbol{\beta} = \boldsymbol{CX\beta}-\boldsymbol{\beta} = (\boldsymbol{CX}-\boldsymbol{I})\boldsymbol{\beta}\) ; pour que \(\tilde{\boldsymbol{\beta}}\) soit sans biais quel que soit \(\boldsymbol{\beta}\), il faut \(\boldsymbol{CX}=\boldsymbol{I}_{p+1}\).

Écrivons \(\boldsymbol{C} = (\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top + \boldsymbol{D}\). Alors \(\boldsymbol{CX} = \boldsymbol{I} + \boldsymbol{DX} = \boldsymbol{I}\) implique \(\boldsymbol{DX}=\boldsymbol{0}\).

La variance de \(\tilde{\boldsymbol{\beta}}\) est \(Var(\tilde{\boldsymbol{\beta}}) = \sigma^2\boldsymbol{CC}^\top\). En substituant et en développant :

\[\boldsymbol{CC}^\top = (\boldsymbol{X}^\top\boldsymbol{X})^{-1} + (\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{D}^\top + \boldsymbol{DX}(\boldsymbol{X}^\top\boldsymbol{X})^{-1} + \boldsymbol{DD}^\top.\]

Puisque \(\boldsymbol{DX}=\boldsymbol{0}\) (et donc \(\boldsymbol{X}^\top\boldsymbol{D}^\top=\boldsymbol{0}\)), les deux termes croisés s’annulent, et

\[Var(\tilde{\boldsymbol{\beta}}) = \sigma^2(\boldsymbol{X}^\top\boldsymbol{X})^{-1} + \sigma^2\boldsymbol{DD}^\top = Var(\hat{\boldsymbol{\beta}}) + \sigma^2\boldsymbol{DD}^\top.\]

Comme \(\boldsymbol{DD}^\top\) est semi-définie positive, \(Var(\tilde{\boldsymbol{\beta}}) \succeq Var(\hat{\boldsymbol{\beta}})\) (au sens des matrices), avec égalité si et seulement si \(\boldsymbol{D}=\boldsymbol{0}\). L’estimateur des moindres carrés est donc BLUE.

Exercice 2. Montrez que \((n-(p+1))\hat\sigma^2/\sigma^2 \sim \chi^2_{n-(p+1)}\) et que cette statistique est indépendante de \(\hat{\boldsymbol{\beta}}\). (Indice : utilisez le théorème de Cochran en exploitant le fait que \(\boldsymbol{H}\) et \(\boldsymbol{I}-\boldsymbol{H}\) sont des projecteurs orthogonaux de rangs \(p+1\) et \(n-(p+1)\), et que \(\hat{\boldsymbol{\beta}}\) et \(\boldsymbol{e}\) sont des transformations linéaires non corrélées, donc indépendantes, du vecteur gaussien \(\boldsymbol{Y}\).)

Solution

Loi du khi-carré. Posons \(\boldsymbol{z}=\boldsymbol{\epsilon}/\sigma \sim \mathcal{N}(\boldsymbol{0},\boldsymbol{I})\). Comme \((\boldsymbol{I}-\boldsymbol{H})\boldsymbol{X}=\boldsymbol{0}\), on a \(SS_E = \boldsymbol{Y}^\top(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{Y} = \boldsymbol{\epsilon}^\top(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{\epsilon} = \sigma^2\boldsymbol{z}^\top(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{z}\). Puisque \(\boldsymbol{I}-\boldsymbol{H}\) est symétrique et idempotente, c’est la matrice d’un projecteur orthogonal de rang \(\text{tr}(\boldsymbol{I}-\boldsymbol{H})=n-(p+1)\). Par décomposition spectrale, \(\boldsymbol{I}-\boldsymbol{H}=\boldsymbol{P}\boldsymbol{\Lambda}\boldsymbol{P}^\top\) avec \(n-(p+1)\) valeurs propres égales à \(1\) et les autres à \(0\), et \(\boldsymbol{P}^\top\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})\). Donc \(\boldsymbol{z}^\top(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{z}\) est une somme de \(n-(p+1)\) carrés de variables \(\mathcal{N}(0,1)\) indépendantes, soit \(\chi^2_{n-(p+1)}\). Ainsi \((n-(p+1))\hat\sigma^2/\sigma^2 = SS_E/\sigma^2 \sim \chi^2_{n-(p+1)}\).

Indépendance. \(\hat{\boldsymbol{\beta}}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{Y}\) et \(\boldsymbol{e}=(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{Y}\) sont deux transformations linéaires du même vecteur gaussien \(\boldsymbol{Y}\), donc conjointement gaussiennes ; il suffit de montrer qu’elles sont non corrélées pour conclure à l’indépendance. On calcule

\[Cov(\hat{\boldsymbol{\beta}},\boldsymbol{e}) = (\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\,\sigma^2\boldsymbol{I}\,(\boldsymbol{I}-\boldsymbol{H})^\top = \sigma^2(\boldsymbol{X}^\top\boldsymbol{X})^{-1}(\boldsymbol{X}^\top - \boldsymbol{X}^\top\boldsymbol{H}).\]

Or \(\boldsymbol{X}^\top\boldsymbol{H} = \boldsymbol{X}^\top\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top=\boldsymbol{X}^\top\), donc \(Cov(\hat{\boldsymbol{\beta}},\boldsymbol{e})=\boldsymbol{0}\). Comme \(\hat\sigma^2\) est fonction de \(\boldsymbol{e}\) (via \(SS_E=\boldsymbol{e}^\top\boldsymbol{e}\)), \(\hat\sigma^2\) est indépendant de \(\hat{\boldsymbol{\beta}}\).

Exercice 3. En partant du rapport des vraisemblances maximisées du modèle réduit (à \(q+1\) paramètres) et du modèle complet (à \(p+1\) paramètres, \(q<p\)), dérivez la statistique

\[F_{obs} = \frac{(SSE_{r\acute{e}duit}-SSE_{complet})/(p-q)}{SSE_{complet}/(n-(p+1))}\]

et montrez qu’elle suit une loi \(F_{p-q,\,n-(p+1)}\) sous l’hypothèse nulle que le modèle réduit est adéquat. Précisez quelles hypothèses de normalité et d’indépendance sont nécessaires.

Solution

Soit \(\boldsymbol{H}_c\) et \(\boldsymbol{H}_r\) les matrices chapeaux des modèles complet et réduit, de rangs respectifs \(p+1\) et \(q+1\). Comme le modèle réduit est emboîté dans le modèle complet, \(\boldsymbol{H}_r\boldsymbol{H}_c = \boldsymbol{H}_c\boldsymbol{H}_r = \boldsymbol{H}_r\), de sorte que \(\boldsymbol{H}_c-\boldsymbol{H}_r\) est elle-même symétrique et idempotente, de rang \((p+1)-(q+1)=p-q\), et orthogonale à \(\boldsymbol{I}-\boldsymbol{H}_c\) (c’est-à-dire \((\boldsymbol{H}_c-\boldsymbol{H}_r)(\boldsymbol{I}-\boldsymbol{H}_c)=\boldsymbol{0}\)).

Sous \(H_0\) (le modèle réduit est adéquat, c’est-à-dire \(E[\boldsymbol{Y}]=\boldsymbol{X}_r\boldsymbol{\beta}_r\)), on montre comme à l’exercice précédent que

\[\frac{SSE_{r\acute{e}duit}-SSE_{complet}}{\sigma^2} = \frac{\boldsymbol{Y}^\top(\boldsymbol{H}_c-\boldsymbol{H}_r)\boldsymbol{Y}}{\sigma^2} \sim \chi^2_{p-q}, \qquad \frac{SSE_{complet}}{\sigma^2}\sim\chi^2_{n-(p+1)},\]

et que ces deux statistiques sont indépendantes (car \(\boldsymbol{H}_c-\boldsymbol{H}_r\) et \(\boldsymbol{I}-\boldsymbol{H}_c\) sont des projecteurs orthogonaux l’un à l’autre, appliqués au même vecteur gaussien). Le rapport de deux variables \(\chi^2\) indépendantes, chacune divisée par son nombre de degrés de liberté, suit par définition une loi \(F\) :

\[F_{obs} = \frac{\{\boldsymbol{Y}^\top(\boldsymbol{H}_c-\boldsymbol{H}_r)\boldsymbol{Y}/\sigma^2\}/(p-q)}{\{SSE_{complet}/\sigma^2\}/(n-(p+1))} = \frac{(SSE_{r\acute{e}duit}-SSE_{complet})/(p-q)}{SSE_{complet}/(n-(p+1))} \sim F_{p-q,\,n-(p+1)}.\]

Les hypothèses nécessaires sont : erreurs indépendantes, identiquement distribuées selon \(\mathcal{N}(0,\sigma^2)\) (même variance dans les deux modèles), et modèles emboîtés (l’espace des moyennes du modèle réduit est inclus dans celui du modèle complet).

Exercice 4. La distance de Cook est définie par \[D_i \propto \sum_{k=1}^n (\hat Y_k - \hat Y_{k(-i)})^2,\] où \(\hat Y_{k(-i)}\) est la valeur ajustée en \(\boldsymbol{X}_k\) obtenue sans l’observation \(i\). À partir de cette définition, retrouvez la formule \[D_i = \frac{r_i^2}{p+1}\cdot\frac{h_i}{1-h_i}.\]

Indice : utilisez la formule de mise à jour des moindres carrés lorsqu’une observation est retirée, \[\hat{\boldsymbol{\beta}}_{(-i)} = \hat{\boldsymbol{\beta}} - \frac{(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}_i e_i}{1-h_i}.\]

Solution

Puisque \(\hat Y_k-\hat Y_{k(-i)}=\boldsymbol{X}_k^\top(\hat{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}_{(-i)})\), la définition se réécrit \[D_i = \frac{\sum_{k=1}^n(\hat Y_k-\hat Y_{k(-i)})^2}{(p+1)\hat\sigma^2} = \frac{(\hat{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}_{(-i)})^\top\boldsymbol{X}^\top\boldsymbol{X}(\hat{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}_{(-i)})}{(p+1)\hat\sigma^2}.\]

En substituant la formule de mise à jour \(\hat{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}_{(-i)} = \dfrac{(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}_i e_i}{1-h_i}\), on obtient

\[(\hat{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}_{(-i)})^\top\boldsymbol{X}^\top\boldsymbol{X}(\hat{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}_{(-i)}) = \frac{e_i^2}{(1-h_i)^2}\,\boldsymbol{X}_i^\top(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}_i = \frac{e_i^2}{(1-h_i)^2}\,\boldsymbol{X}_i^\top(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}_i.\]

Or \(\boldsymbol{X}_i^\top(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}_i = h_i\) (c’est l’élément diagonal \(i\) de \(\boldsymbol{H}\)), donc l’expression ci-dessus vaut \(\dfrac{e_i^2 h_i}{(1-h_i)^2}\). On obtient

\[D_i = \frac{e_i^2 h_i}{(p+1)\hat\sigma^2(1-h_i)^2} = \frac{e_i^2}{\hat\sigma^2(1-h_i)}\cdot\frac{h_i}{(p+1)(1-h_i)} = \frac{r_i^2}{p+1}\cdot\frac{h_i}{1-h_i},\]

en reconnaissant le résidu standardisé \(r_i = e_i/\{\hat\sigma\sqrt{1-h_i}\}\), donc \(r_i^2 = e_i^2/\{\hat\sigma^2(1-h_i)\}\).

Chapitre 3 : Modèles de régression pour données binaires

Niveau facile

Exercice 1. Pour chacune des trois fonctions de lien présentées dans ce chapitre (logit, probit, c-log-log), vérifiez que \(g^{-1}\) est bien l’inverse de \(g\), c’est-à-dire que \(g\{g^{-1}(\eta)\}=\eta\) pour tout \(\eta \in \mathbb{R}\).

Solution

Logit. On a \(g^{-1}(\eta) = \dfrac{e^\eta}{1+e^\eta} = \pi\), donc \(1-\pi = \dfrac{1}{1+e^\eta}\) et \[g(\pi) = \log\left(\frac{\pi}{1-\pi}\right) = \log(e^\eta) = \eta.\]

Probit. \(g^{-1}(\eta)=\Phi(\eta)\) et \(g(u)=\Phi^{-1}(u)\) sont par construction des fonctions réciproques l’une de l’autre (\(\Phi\) étant une bijection strictement croissante de \(\mathbb{R}\) vers \((0,1)\)), donc \(g\{g^{-1}(\eta)\} = \Phi^{-1}\{\Phi(\eta)\} = \eta\) directement.

C-log-log. Si \(g(u)=\log\{-\log(1-u)\}=\eta\), alors \(-\log(1-u)=e^\eta\), donc \(1-u=e^{-e^\eta}\) et \(u=g^{-1}(\eta)=1-e^{-e^\eta}=\pi\). En substituant, \[g(\pi) = \log\{-\log(1-\pi)\} = \log\{-\log(e^{-e^\eta})\} = \log(e^\eta) = \eta.\]

Exercice 2. Avec le lien logit, la cote s’écrit \(o_i = \pi_i/(1-\pi_i) = e^{\eta_i}\). Calculez \(o_i\) et \(\pi_i\) pour \(\eta_i = 0\), \(\eta_i = \ln 2\) et \(\eta_i = -\ln 3\).

Solution

Puisque \(o_i = e^{\eta_i}\) et \(\pi_i = o_i/(1+o_i)\) :

• \(\eta_i = 0\) : \(o_i = e^0 = 1\), donc \(\pi_i = 1/2\) (autant de chances de survivre que de ne pas survivre) ;
• \(\eta_i = \ln 2\) : \(o_i = e^{\ln 2} = 2\), donc \(\pi_i = 2/3\) ;
• \(\eta_i = -\ln 3\) : \(o_i = e^{-\ln 3} = 1/3\), donc \(\pi_i = (1/3)/(4/3) = 1/4\).

Exercice 3. Soit \(Y_i \sim \mathrm{Bin}(m_i,\pi_i)\), avec \(\mu_i = m_i\pi_i\). Montrez que \[Var[Y_i] = \frac{\mu_i(m_i-\mu_i)}{m_i}\] est maximale, pour \(m_i\) fixé, lorsque \(\pi_i = 1/2\). Quel est le lien avec le poids \(w_i = m_i\pi_i(1-\pi_i)\) utilisé dans l’estimation par IRLS ?

Solution

Puisque \(\mu_i = m_i\pi_i\), on a \(m_i - \mu_i = m_i(1-\pi_i)\), donc \[\frac{\mu_i(m_i-\mu_i)}{m_i} = \frac{m_i\pi_i \cdot m_i(1-\pi_i)}{m_i} = m_i\pi_i(1-\pi_i) = Var[Y_i].\]

En dérivant \(\pi_i(1-\pi_i)\) par rapport à \(\pi_i\) : \(\dfrac{d}{d\pi_i}\{\pi_i(1-\pi_i)\} = 1-2\pi_i\), qui s’annule en \(\pi_i=1/2\). La dérivée seconde vaut \(-2<0\), donc il s’agit bien d’un maximum, où \(Var[Y_i]=m_i/4\).

Le poids \(w_i = m_i\pi_i(1-\pi_i)\) est exactement \(Var[Y_i]\) : les observations dont la probabilité prédite est proche de \(1/2\) (les plus incertaines) reçoivent donc le plus grand poids dans l’algorithme IRLS, tandis que les observations dont \(\pi_i\) est proche de \(0\) ou \(1\) (presque certaines) reçoivent un poids proche de \(0\) — ce qui explique en partie l’instabilité de l’estimation en présence de séparation, où de nombreux \(\pi_i\) s’approchent de \(0\) ou \(1\).

Niveau moyen

Exercice 1. En partant de la log-vraisemblance avec lien logit, \[\ell(\boldsymbol{\beta}; \boldsymbol{y}) = \sum_{i=1}^n \left\{ \log\binom{m_i}{y_i} + y_i \eta_i - m_i \log\left(1+e^{\eta_i}\right) \right\}, \qquad \eta_i = \beta_0+\beta_1 X_{i1}+\cdots+\beta_p X_{ip},\] dérivez le score \(\boldsymbol{U}(\boldsymbol{\beta}) = \partial\ell/\partial\boldsymbol{\beta}\) et montrez qu’il s’écrit, sous forme matricielle, \(\boldsymbol{U}(\boldsymbol{\beta}) = \boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\).

Solution

Pour \(j=0,1,\ldots,p\), avec \(\partial\eta_i/\partial\beta_j = X_{ij}\) (et \(X_{i0}\equiv 1\)) : le terme \(\log\binom{m_i}{y_i}\) ne dépend pas de \(\boldsymbol\beta\), donc sa dérivée est nulle. Pour \(y_i\eta_i\), \(\partial(y_i\eta_i)/\partial\beta_j = y_iX_{ij}\). Pour le dernier terme, par dérivation en chaîne, \[\frac{\partial}{\partial\beta_j}\log(1+e^{\eta_i}) = \frac{e^{\eta_i}}{1+e^{\eta_i}}\cdot X_{ij} = \pi_i X_{ij},\] donc \(\partial\{m_i\log(1+e^{\eta_i})\}/\partial\beta_j = m_i\pi_i X_{ij}\). En combinant, \[\frac{\partial\ell}{\partial\beta_j} = \sum_{i=1}^n (y_i - m_i\pi_i)X_{ij} = \sum_{i=1}^n (y_i-\mu_i)X_{ij}.\]

En empilant ces \(p+1\) dérivées partielles et en reconnaissant que \(X_{ij}\) est l’élément \((i,j)\) de \(\boldsymbol{X}\), on obtient \(\boldsymbol{U}(\boldsymbol{\beta}) = \boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\).

Exercice 2. En dérivant le résultat de l’exercice précédent une fois de plus par rapport à \(\boldsymbol{\beta}\), montrez que la matrice d’information de Fisher pour le lien logit s’écrit \(I(\boldsymbol{\beta}) = \boldsymbol{X}^\top\boldsymbol{W}\boldsymbol{X}\), où \(\boldsymbol{W} = \mathrm{diag}(w_1,\ldots,w_n)\) et \(w_i = m_i\pi_i(1-\pi_i)\).

Solution

On a \(U_j(\boldsymbol\beta) = \sum_i (y_i-\mu_i)X_{ij}\), où \(\mu_i = m_i\pi_i = m_i\, e^{\eta_i}/(1+e^{\eta_i})\). Puisque \(y_i\) ne dépend pas de \(\boldsymbol\beta\), \[\frac{\partial U_j}{\partial\beta_k} = -\sum_i \frac{\partial\mu_i}{\partial\beta_k}X_{ij} = -\sum_i m_i\frac{\partial\pi_i}{\partial\eta_i}\cdot\frac{\partial\eta_i}{\partial\beta_k}\,X_{ij}.\]

Comme \(\pi_i = e^{\eta_i}/(1+e^{\eta_i})\) est la fonction logistique, \(\partial\pi_i/\partial\eta_i = \pi_i(1-\pi_i)\) (propriété classique de la fonction logistique), et \(\partial\eta_i/\partial\beta_k = X_{ik}\). Donc \[\frac{\partial U_j}{\partial\beta_k} = -\sum_i m_i\pi_i(1-\pi_i)\,X_{ij}X_{ik} = -\sum_i w_i X_{ij}X_{ik}.\]

L’élément \((j,k)\) de l’information de Fisher est \(I_{jk} = -\partial U_j/\partial\beta_k = \sum_i w_i X_{ij}X_{ik}\), ce qui correspond exactement à l’élément \((j,k)\) de \(\boldsymbol{X}^\top\boldsymbol{W}\boldsymbol{X}\). Donc \(I(\boldsymbol\beta) = \boldsymbol{X}^\top\boldsymbol{W}\boldsymbol{X}\).

Exercice 3. Avec un seul prédicteur continu \(x\) et le lien logit, \(\hat\pi(x) = g^{-1}(\hat\beta_0+\hat\beta_1 x)\). Montrez que \[\hat\pi'(x) = \hat\beta_1\,\hat\pi(x)\{1-\hat\pi(x)\}\] et que cette pente est maximale au point d’inflexion \(\hat\pi(x)=1/2\), où elle vaut \(\hat\beta_1/4\).

Solution

Posons \(t = \hat\beta_0+\hat\beta_1 x\) et \(s(t) = e^t/(1+e^t)\), de sorte que \(\hat\pi(x)=s(t)\). On vérifie d’abord que \(s'(t) = s(t)\{1-s(t)\}\) : \[s'(t) = \frac{e^t(1+e^t) - e^t\cdot e^t}{(1+e^t)^2} = \frac{e^t}{(1+e^t)^2} = \frac{e^t}{1+e^t}\cdot\frac{1}{1+e^t} = s(t)\{1-s(t)\}.\]

Par dérivation en chaîne, \(\hat\pi'(x) = s'(t)\cdot\dfrac{dt}{dx} = \hat\beta_1\, s(t)\{1-s(t)\} = \hat\beta_1\,\hat\pi(x)\{1-\hat\pi(x)\}\).

Comme montré à l’exercice 3 du niveau facile, \(u(1-u)\) (pour \(u=\hat\pi(x)\)) est maximal en \(u=1/2\), où il vaut \(1/4\). Donc \(\hat\pi'(x)\) est maximale lorsque \(\hat\pi(x)=1/2\), où elle vaut \(\hat\beta_1/4\).

Niveau difficile

Exercice 1. Montrez que, lorsque les \(m_i\to\infty\), la déviance individuelle \(d_i\) est asymptotiquement équivalente au carré du résidu de Pearson \(r_{P,i}^2\), c’est-à-dire que \[d_i = 2\left\{y_i\log\!\left(\frac{y_i}{\hat\mu_i}\right) + (m_i-y_i)\log\!\left(\frac{m_i-y_i}{m_i-\hat\mu_i}\right)\right\} \approx \frac{(y_i-\hat\mu_i)^2\,m_i}{\hat\mu_i(m_i-\hat\mu_i)} = r_{P,i}^2.\] En déduire que \(D \approx \chi^2_P\) pour des \(m_i\) grands. (Indice : posez \(\delta = y_i-\hat\mu_i\) et effectuez un développement de Taylor à l’ordre 2 de \(\log(1+\delta/\hat\mu_i)\) et \(\log(1-\delta/(m_i-\hat\mu_i))\) autour de \(\delta=0\).)

Solution

Posons \(\delta = y_i - \hat\mu_i\), de sorte que \(y_i = \hat\mu_i+\delta\) et \(m_i-y_i = (m_i-\hat\mu_i)-\delta\). En utilisant \(\log(1+t) \approx t - t^2/2\) pour \(t\) petit :

\[y_i\log\left(\frac{y_i}{\hat\mu_i}\right) = (\hat\mu_i+\delta)\log\left(1+\frac{\delta}{\hat\mu_i}\right) \approx (\hat\mu_i+\delta)\left(\frac{\delta}{\hat\mu_i}-\frac{\delta^2}{2\hat\mu_i^2}\right) \approx \delta + \frac{\delta^2}{2\hat\mu_i},\]

en ne conservant que les termes jusqu’à l’ordre \(\delta^2\). De même,

\[(m_i-y_i)\log\left(\frac{m_i-y_i}{m_i-\hat\mu_i}\right) = \{(m_i-\hat\mu_i)-\delta\}\log\left(1-\frac{\delta}{m_i-\hat\mu_i}\right) \approx -\delta + \frac{\delta^2}{2(m_i-\hat\mu_i)}.\]

En additionnant les deux approximations, les termes en \(\delta\) s’annulent :

\[y_i\log\left(\frac{y_i}{\hat\mu_i}\right) + (m_i-y_i)\log\left(\frac{m_i-y_i}{m_i-\hat\mu_i}\right) \approx \frac{\delta^2}{2}\left(\frac{1}{\hat\mu_i}+\frac{1}{m_i-\hat\mu_i}\right) = \frac{\delta^2}{2}\cdot\frac{m_i}{\hat\mu_i(m_i-\hat\mu_i)}.\]

En multipliant par \(2\) (définition de \(d_i\)), on obtient \[d_i \approx \frac{\delta^2\, m_i}{\hat\mu_i(m_i-\hat\mu_i)} = \frac{(y_i-\hat\mu_i)^2\,m_i}{\hat\mu_i(m_i-\hat\mu_i)} = r_{P,i}^2.\]

En sommant sur \(i=1,\ldots,n\), \(D = \sum_i d_i \approx \sum_i r_{P,i}^2 = \chi^2_P\). Cette approximation devient exacte lorsque \(\delta/\hat\mu_i\) et \(\delta/(m_i-\hat\mu_i)\) tendent vers \(0\), ce qui se produit quand \(m_i\to\infty\) (à \(\pi_i\) fixé).

Exercice 2. Montrez que, pour deux modèles emboîtés de déviances \(D_{r\acute{e}duit}\) et \(D_{complet}\), la statistique \(\chi^2 = D_{r\acute{e}duit}-D_{complet}\) est exactement égale à la statistique du test du rapport de vraisemblance \[G^2 = -2\left\{\ell(\text{modèle r\'eduit}) - \ell(\text{modèle complet})\right\}.\] Expliquez pourquoi le terme associé au modèle saturé n’apparaît pas dans cette différence, alors qu’il intervient dans la définition de chaque déviance prise individuellement.

Solution

Par définition, \(D_{r\acute{e}duit} = 2\{\ell(\boldsymbol y,\boldsymbol y) - \ell(\hat{\boldsymbol\mu}_r;\boldsymbol y)\}\) et \(D_{complet} = 2\{\ell(\boldsymbol y,\boldsymbol y)-\ell(\hat{\boldsymbol\mu}_c;\boldsymbol y)\}\), où \(\ell(\boldsymbol y,\boldsymbol y)\) est la log-vraisemblance du modèle saturé (la même dans les deux expressions, puisqu’elle ne dépend que des données observées \(\boldsymbol y\), pas du modèle ajusté). En soustrayant,

\[D_{r\acute{e}duit}-D_{complet} = 2\{\ell(\boldsymbol y,\boldsymbol y)-\ell(\hat{\boldsymbol\mu}_r;\boldsymbol y)\} - 2\{\ell(\boldsymbol y,\boldsymbol y)-\ell(\hat{\boldsymbol\mu}_c;\boldsymbol y)\} = 2\{\ell(\hat{\boldsymbol\mu}_c;\boldsymbol y) - \ell(\hat{\boldsymbol\mu}_r;\boldsymbol y)\},\]

ce qui est exactement \(-2\{\ell(\text{r\'eduit})-\ell(\text{complet})\} = G^2\). Le terme \(\ell(\boldsymbol y,\boldsymbol y)\) du modèle saturé s’annule dans la différence car il est commun aux deux déviances : il ne dépend que des données, identiques pour les deux modèles ajustés sur le même échantillon. C’est précisément pour cette raison que comparer des déviances revient à comparer des vraisemblances, sans jamais avoir à évaluer explicitement \(\ell(\boldsymbol y,\boldsymbol y)\).

Exercice 3. Considérez un modèle logistique à un seul prédicteur \(x\), sans ordonnée à l’origine (\(\eta_i = \beta x_i\)), et supposez qu’il y a séparation complète au sens où \(y_i=1\) si \(x_i>0\) et \(y_i=0\) si \(x_i<0\) (aucun \(x_i\) n’est nul). Montrez que le score \(U(\beta) = \sum_i (y_i-\pi_i)x_i\) est strictement positif pour toute valeur finie de \(\beta\), et concluez que l’estimateur du maximum de vraisemblance \(\hat\beta\) n’existe pas (c’est-à-dire que \(\ell(\beta)\) n’atteint son supremum qu’à la limite \(\beta\to\infty\)).

Solution

Fixons \(\beta\in\mathbb{R}\) et examinons le terme \((y_i-\pi_i)x_i\) pour chaque observation \(i\), où \(\pi_i = e^{\beta x_i}/(1+e^{\beta x_i}) \in (0,1)\) pour tout \(\beta\) fini.

• Si \(x_i>0\), alors \(y_i=1\) par hypothèse de séparation, donc \((y_i-\pi_i)x_i = (1-\pi_i)x_i\). Comme \(1-\pi_i>0\) et \(x_i>0\), ce terme est strictement positif.
• Si \(x_i<0\), alors \(y_i=0\), donc \((y_i-\pi_i)x_i = -\pi_i x_i\). Comme \(\pi_i>0\) et \(x_i<0\), on a \(-\pi_i x_i>0\) : ce terme est aussi strictement positif.

Chaque terme de la somme \(U(\beta)=\sum_i(y_i-\pi_i)x_i\) est donc strictement positif, quel que soit \(\beta\in\mathbb{R}\) fini, donc \(U(\beta)>0\) pour tout \(\beta\).

Puisque \(U(\beta) = \ell'(\beta)\), la log-vraisemblance \(\ell\) est strictement croissante sur tout \(\mathbb{R}\) : elle n’a donc aucun point stationnaire (\(U(\beta)=0\) n’a pas de solution finie), et \(\sup_{\beta\in\mathbb{R}} \ell(\beta) = \lim_{\beta\to\infty}\ell(\beta)\). Ce supremum n’est donc jamais atteint pour une valeur finie de \(\beta\) : l’estimateur du maximum de vraisemblance \(\hat\beta\) n’existe pas au sens usuel, ce qui confirme formellement, dans ce cas simple, le phénomène décrit qualitativement dans le chapitre.

Exercice 4. On considère un modèle logistique à un seul prédicteur \(x\) (sans ordonnée à l’origine, données individuelles \(m_i=1\)), pour lequel l’information de Fisher est \(I(\beta) = \sum_{i=1}^n w_i x_i^2\), avec \(w_i=\pi_i(1-\pi_i)\). Supposons que les données présentent une séparation complète comme à l’exercice précédent, et que l’algorithme d’estimation s’arrête à une valeur \(\hat\beta\) grande mais finie (à cause du critère d’arrêt numérique). Montrez que, pour \(x_i\neq 0\) fixés, \(w_i\) décroît exponentiellement vite vers \(0\) quand \(\hat\beta\to\infty\), et expliquez pourquoi cela implique que la statistique de Wald \(z_{obs} = \hat\beta\sqrt{I(\hat\beta)}\) tend vers \(0\) plutôt que vers l’infini — c’est-à-dire le phénomène de Hauck-Donner.

Solution

Pour \(x_i>0\) fixé, \(\pi_i = \dfrac{1}{1+e^{-\hat\beta x_i}}\), donc \(1-\pi_i = \dfrac{e^{-\hat\beta x_i}}{1+e^{-\hat\beta x_i}} \approx e^{-\hat\beta x_i}\) quand \(\hat\beta\to\infty\) (puisque le dénominateur tend vers \(1\)). Comme \(\pi_i\to 1\), on a \[w_i = \pi_i(1-\pi_i) \approx e^{-\hat\beta x_i} \longrightarrow 0\] à un taux exponentiel en \(\hat\beta\). Un calcul symétrique pour \(x_i<0\) donne \(w_i \approx e^{-\hat\beta|x_i|}\). Dans les deux cas, \(w_i\) décroît exponentiellement vite vers \(0\) quand \(\hat\beta\to\infty\), beaucoup plus vite que n’importe quelle décroissance polynomiale (en \(1/\hat\beta\), \(1/\hat\beta^2\), etc.).

Par conséquent, \[I(\hat\beta) = \sum_{i=1}^n w_i x_i^2 \approx \sum_{i=1}^n x_i^2\, e^{-\hat\beta|x_i|} \longrightarrow 0\] également à un taux exponentiel. La statistique de Wald s’écrit \[z_{obs} = \frac{\hat\beta}{\sqrt{\hat{\text{V}}(\hat\beta)}} = \hat\beta\sqrt{I(\hat\beta)}.\] Même si \(\hat\beta\) croît (de façon au plus polynomiale, voire logarithmique, puisqu’il est borné par le critère d’arrêt numérique), le facteur \(\sqrt{I(\hat\beta)}\) décroît exponentiellement vite vers \(0\). Une décroissance exponentielle domine toujours une croissance polynomiale (ou plus lente) : le produit \(\hat\beta\sqrt{I(\hat\beta)}\) tend donc vers \(0\), et non vers l’infini, même si l’effet réel \(\beta\) est très grand. C’est exactement le phénomène de Hauck-Donner : plus la séparation est marquée (plus \(\hat\beta\) est grand), plus la statistique de Wald \(z_{obs}\) s’effondre vers \(0\), donnant l’illusion d’un effet non significatif.

Chapitre 4 : Modèles linéaires généralisés (GLM)

Niveau facile

Exercice 1. Soit \(Y\sim\mathrm{Exp}(\lambda)\) (loi exponentielle de taux \(\lambda>0\)), de densité \(f(y)=\lambda e^{-\lambda y}\), \(y>0\). Mettez cette densité sous la forme de la famille exponentielle, identifiez \(\theta\), \(\phi\), \(a(y,\phi)\) et \(b(\theta)\), puis retrouvez \(E[Y]=1/\lambda\) et \(Var[Y]=1/\lambda^2\) à l’aide de \(b'(\theta)\) et \(b''(\theta)\).

Solution

On écrit \(f(y) = \exp\{-\lambda y + \log\lambda\}\). En posant \(\theta=-\lambda\) (donc \(\lambda=-\theta\)), le terme \(-\lambda y = y\theta\) et \(\log\lambda = \log(-\theta)\). On a donc \[f(y) = \exp\{y\theta - b(\theta)\}, \qquad b(\theta) = -\log(-\theta), \qquad \phi=1, \qquad a(y,\phi)=1.\]

On vérifie que \(b'(\theta) = -1/\theta\) et \(b''(\theta) = 1/\theta^2\). En substituant \(\theta=-\lambda\) : \[E[Y] = b'(\theta) = -\frac{1}{-\lambda} = \frac{1}{\lambda}, \qquad Var[Y] = \phi\,b''(\theta) = \frac{1}{(-\lambda)^2} = \frac{1}{\lambda^2},\] ce qui correspond bien aux moments usuels de la loi exponentielle. (On notera que \(b(\theta)=-\log(-\theta)\) est la même fonction cumulante que pour la loi gamma, puisque la loi exponentielle est le cas particulier \(\alpha=1\) de la gamma.)

Exercice 2. Pour la loi binomiale, on a \(b(\theta) = m\log(1+e^\theta)\), où \(\theta=\mathrm{logit}(\pi)\). En calculant \(b'(\theta)\) et \(b''(\theta)\), puis en substituant \(\pi=\mu/m\) et \(1-\pi=(m-\mu)/m\), retrouvez la fonction de variance \(V(\mu) = \mu(m-\mu)/m\).

Solution

On a \(b'(\theta) = m\cdot\dfrac{e^\theta}{1+e^\theta} = m\pi = \mu\), et en dérivant une seconde fois (dérivée de la fonction logistique) : \[b''(\theta) = m\pi(1-\pi).\]

En substituant \(\pi = \mu/m\) et \(1-\pi = (m-\mu)/m\) : \[V(\mu) = b''(\theta) = m\cdot\frac{\mu}{m}\cdot\frac{m-\mu}{m} = \frac{\mu(m-\mu)}{m}.\]

Exercice 3. Pour la loi de Poisson, \(b(\theta) = e^\theta\). Un lien \(g\) est dit canonique s’il satisfait \(g\{b'(u)\}=u\) pour tout \(u\). Montrez que le lien canonique de la loi de Poisson est bien \(g(\mu)=\log(\mu)\).

Solution

On a \(b'(\theta) = e^\theta\). En posant \(\mu = b'(u) = e^u\), la condition de lien canonique \(g\{b'(u)\}=u\) devient \(g(e^u) = u\). En posant \(\mu=e^u\), c’est-à-dire \(u = \log(\mu)\), on obtient \(g(\mu) = \log(\mu)\), ce qui est bien le lien canonique annoncé pour la loi de Poisson.

Niveau moyen

Exercice 1. On admet les deux identités générales de la théorie de la vraisemblance, appliquées ici au paramètre canonique \(\theta\) : \(E\!\left[\dfrac{\partial\log f}{\partial\theta}\right]=0\) et \(Var\!\left[\dfrac{\partial\log f}{\partial\theta}\right] = -E\!\left[\dfrac{\partial^2\log f}{\partial\theta^2}\right]\). En partant de \(\log f(y) = \log a(y,\phi) + (y\theta-b(\theta))/\phi\), utilisez ces deux identités pour retrouver \(E[Y]=b'(\theta)\) et \(Var[Y]=\phi\,b''(\theta)\).

Solution

On calcule \(\dfrac{\partial\log f}{\partial\theta} = \dfrac{y-b'(\theta)}{\phi}\). La première identité donne \[E\left[\frac{Y-b'(\theta)}{\phi}\right] = 0 \quad\Longrightarrow\quad E[Y] = b'(\theta).\]

En dérivant une seconde fois, \(\dfrac{\partial^2\log f}{\partial\theta^2} = -\dfrac{b''(\theta)}{\phi}\), donc \(-E\!\left[\dfrac{\partial^2\log f}{\partial\theta^2}\right] = \dfrac{b''(\theta)}{\phi}\). La seconde identité donne \[Var\left[\frac{Y-b'(\theta)}{\phi}\right] = \frac{Var[Y]}{\phi^2} = \frac{b''(\theta)}{\phi} \quad\Longrightarrow\quad Var[Y] = \phi\, b''(\theta).\]

Exercice 2. En partant de la log-vraisemblance individuelle \(\ell_i(\theta_i) = \log a(y_i,\phi) + (y_i\theta_i-b(\theta_i))/\phi\), où \(\theta_i\) dépend de \(\boldsymbol\beta\) via \(\mu_i=g^{-1}(\eta_i)\) puis \(\theta_i=(b')^{-1}(\mu_i)\), dérivez le score général \[\frac{\partial\ell}{\partial\beta_j} = \frac{1}{\phi}\sum_{i=1}^n \frac{y_i-\mu_i}{V(\mu_i)\,g'(\mu_i)}\,X_{ij}\] à l’aide de la règle de dérivation en chaîne \(\dfrac{\partial\ell_i}{\partial\beta_j} = \dfrac{\partial\ell_i}{\partial\theta_i}\cdot\dfrac{\partial\theta_i}{\partial\mu_i}\cdot\dfrac{\partial\mu_i}{\partial\eta_i}\cdot\dfrac{\partial\eta_i}{\partial\beta_j}\). Montrez ensuite que, pour le lien canonique, cette expression se simplifie en \(U(\boldsymbol\beta) = \boldsymbol{X}^\top(\boldsymbol y-\boldsymbol\mu)/\phi\).

Solution

Chaque facteur de la règle en chaîne se calcule comme suit :

• \(\dfrac{\partial\ell_i}{\partial\theta_i} = \dfrac{y_i-b'(\theta_i)}{\phi} = \dfrac{y_i-\mu_i}{\phi}\) ;
• puisque \(\mu_i = b'(\theta_i)\), on a \(\dfrac{\partial\mu_i}{\partial\theta_i} = b''(\theta_i) = V(\mu_i)\), donc \(\dfrac{\partial\theta_i}{\partial\mu_i} = \dfrac{1}{V(\mu_i)}\) ;
• puisque \(g(\mu_i)=\eta_i\), on a \(\dfrac{\partial\eta_i}{\partial\mu_i}=g'(\mu_i)\), donc \(\dfrac{\partial\mu_i}{\partial\eta_i} = \dfrac{1}{g'(\mu_i)}\) ;
• \(\dfrac{\partial\eta_i}{\partial\beta_j} = X_{ij}\).

En multipliant les quatre facteurs, \[\frac{\partial\ell_i}{\partial\beta_j} = \frac{y_i-\mu_i}{\phi}\cdot\frac{1}{V(\mu_i)}\cdot\frac{1}{g'(\mu_i)}\cdot X_{ij},\] et en sommant sur \(i\), on obtient la formule annoncée.

Cas du lien canonique. Le lien canonique satisfait \(g=(b')^{-1}\), c’est-à-dire \(\theta = g(\mu)\). En dérivant \(\mu=b'\{g(\mu)\}=\mu\) par rapport à \(\mu\) (identité triviale), on obtient \(1 = b''\{g(\mu)\}\,g'(\mu) = V(\mu)\,g'(\mu)\), donc \(V(\mu)g'(\mu)=1\). En substituant dans le score général, le facteur \(1/\{V(\mu_i)g'(\mu_i)\}\) vaut \(1\) pour chaque \(i\), d’où \[U(\boldsymbol\beta) = \frac{1}{\phi}\sum_{i=1}^n (y_i-\mu_i)X_{ij} \quad\Longrightarrow\quad \boldsymbol U(\boldsymbol\beta) = \frac{1}{\phi}\boldsymbol X^\top(\boldsymbol y-\boldsymbol\mu).\]

Ceci explique pourquoi le score s’est simplifié de façon identique (à \(\phi\) près) pour la régression logistique et la régression de Poisson vues aux chapitres précédents : les liens logit et log sont les liens canoniques de leurs lois respectives.

Exercice 3. Pour la loi gamma, le lien canonique est \(g(\mu) = -1/\mu\) et la fonction de variance est \(V(\mu)=\mu^2\). Vérifiez que \(V(\mu)\,g'(\mu) = 1\), conformément à la propriété générale établie à l’exercice précédent pour tout lien canonique.

Solution

On a \(g(\mu) = -1/\mu = -\mu^{-1}\), donc \(g'(\mu) = \mu^{-2} = 1/\mu^2\). En multipliant par \(V(\mu)=\mu^2\) : \[V(\mu)\,g'(\mu) = \mu^2 \cdot \frac{1}{\mu^2} = 1,\] ce qui confirme que le lien \(g(\mu)=-1/\mu\) est bien le lien canonique de la loi gamma, et que le score associé a la forme simplifiée \(U(\boldsymbol\beta) = \boldsymbol X^\top(\boldsymbol y-\boldsymbol\mu)/\phi\).

Niveau difficile

Exercice 1. Le chapitre affirme, sans démonstration, que la loi asymptotique de \(\hat{\boldsymbol\beta}\) a pour matrice de variance \(\phi\,(\boldsymbol X^\top\boldsymbol W\boldsymbol X)^{-1}\), où \(W_i = 1/\{V(\hat\mu_i)[g'(\hat\mu_i)]^2\}\). En partant du score général \[U_j(\boldsymbol\beta) = \frac{1}{\phi}\sum_{i=1}^n \frac{y_i-\mu_i}{V(\mu_i)\,g'(\mu_i)}\,X_{ij}\] dérivé à l’exercice moyen 2, et en utilisant le fait que \(E[Y_i-\mu_i]=0\) pour annuler les termes où la dérivée tombe sur \(\mu_i\) (information de Fisher espérée, plutôt qu’observée), montrez que \(I(\boldsymbol\beta) = -E\!\left[\dfrac{\partial\boldsymbol U}{\partial\boldsymbol\beta^\top}\right] = \dfrac{1}{\phi}\boldsymbol X^\top\boldsymbol W\boldsymbol X\).

Solution

Posons \(D_i = 1/\{V(\mu_i)g'(\mu_i)\}\), de sorte que \(U_j(\boldsymbol\beta) = \frac1\phi\sum_i (y_i-\mu_i)D_i X_{ij}\). En dérivant par rapport à \(\beta_k\), et puisque \(\mu_i\) (et donc \(D_i\)) dépend de \(\boldsymbol\beta\), la règle du produit donne deux termes : \[\frac{\partial U_j}{\partial\beta_k} = \frac1\phi\sum_i\left[-\frac{\partial\mu_i}{\partial\beta_k}D_i X_{ij} + (y_i-\mu_i)\frac{\partial D_i}{\partial\beta_k}X_{ij}\right].\]

En prenant l’espérance, le second terme s’annule car \(E[y_i-\mu_i]=0\) (et \(\partial D_i/\partial\beta_k\) ne dépend pas de \(y_i\)) : \[E\left[\frac{\partial U_j}{\partial\beta_k}\right] = -\frac1\phi\sum_i D_i X_{ij}\,\frac{\partial\mu_i}{\partial\beta_k}.\]

Par la règle en chaîne, \(\dfrac{\partial\mu_i}{\partial\beta_k} = \dfrac{\partial\mu_i}{\partial\eta_i}\cdot\dfrac{\partial\eta_i}{\partial\beta_k} = \dfrac{X_{ik}}{g'(\mu_i)}\) (comme à l’exercice moyen 2). En substituant, \[E\left[\frac{\partial U_j}{\partial\beta_k}\right] = -\frac1\phi\sum_i \frac{1}{V(\mu_i)g'(\mu_i)}\cdot\frac{1}{g'(\mu_i)}\,X_{ij}X_{ik} = -\frac1\phi\sum_i \frac{X_{ij}X_{ik}}{V(\mu_i)[g'(\mu_i)]^2}.\]

En posant \(w_i = 1/\{V(\mu_i)[g'(\mu_i)]^2\}\), l’élément \((j,k)\) de l’information de Fisher est \[I_{jk}(\boldsymbol\beta) = -E\left[\frac{\partial U_j}{\partial\beta_k}\right] = \frac1\phi\sum_i w_i X_{ij}X_{ik},\] ce qui correspond exactement à l’élément \((j,k)\) de \(\frac1\phi \boldsymbol X^\top\boldsymbol W\boldsymbol X\), où \(\boldsymbol W=\mathrm{diag}(w_1,\ldots,w_n)\). Donc \(I(\boldsymbol\beta) = \frac1\phi\boldsymbol X^\top\boldsymbol W\boldsymbol X\), et \(\widehat{\text{Var}}(\hat{\boldsymbol\beta}) = I(\boldsymbol\beta)^{-1} = \phi\,(\boldsymbol X^\top\boldsymbol W\boldsymbol X)^{-1}\).

Exercice 2. Démontrez formellement la convergence de la loi binomiale vers la loi de Poisson : si \(m\to\infty\) et \(\pi\to 0\) de sorte que \(\mu=m\pi\) demeure fixe, montrez que, pour tout \(y\) fixé, \[\binom{m}{y}\pi^y(1-\pi)^{m-y} \longrightarrow \frac{e^{-\mu}\mu^y}{y!}.\]

Solution

En substituant \(\pi=\mu/m\) : \[\begin{align} \binom{m}{y}\pi^y(1-\pi)^{m-y} &=\frac{m(m-1)\cdots(m-y+1)}{y!}\left(\frac{\mu}{m}\right)^y\left(1-\frac{\mu}{m}\right)^{m-y}\\ &= \frac{\mu^y}{y!}\cdot\underbrace{\frac{m(m-1)\cdots(m-y+1)}{m^y}}_{(A)}\cdot\underbrace{\left(1-\frac{\mu}{m}\right)^{m-y}}_{(B)}. \end{align}\]

Terme (A). Pour \(y\) fixé, \[\frac{m(m-1)\cdots(m-y+1)}{m^y} = \prod_{k=0}^{y-1}\left(1-\frac{k}{m}\right) \longrightarrow 1 \quad\text{quand } m\to\infty,\] puisque chaque facteur tend vers \(1\) et qu’il y en a un nombre fini (\(y\) facteurs).

Terme (B). On écrit \((1-\mu/m)^{m-y} = (1-\mu/m)^m \cdot (1-\mu/m)^{-y}\). Le second facteur tend vers \(1\) (car \(\mu/m\to 0\)), et le premier est la limite classique \[\left(1-\frac{\mu}{m}\right)^m \longrightarrow e^{-\mu} \quad\text{quand } m\to\infty.\] Donc \((B) \to e^{-\mu}\).

En combinant (A) et (B), \[\binom{m}{y}\pi^y(1-\pi)^{m-y} \longrightarrow \frac{\mu^y}{y!}\cdot 1 \cdot e^{-\mu} = \frac{e^{-\mu}\mu^y}{y!},\] ce qui est bien la fonction de masse de la loi \(\mathcal{P}(\mu)\).

Exercice 3. Le test de surdispersion de \(H_0:\alpha=0\) (Poisson) contre \(H_1:\alpha>0\) (binomiale négative) utilise la statistique \(Q=D_0-D_1\), dont la loi asymptotique sous \(H_0\) est le mélange \(\frac12\chi^2_0+\frac12\chi^2_1\) plutôt qu’une \(\chi^2_1\) usuelle. Esquissez un argument justifiant ce résultat, sachant que : (i) l’estimateur du maximum de vraisemblance non contraint de \(\alpha\), \(\hat\alpha_{MV}\), serait asymptotiquement \(\mathcal N(0,v)\) sous \(H_0\) (par la théorie usuelle de la vraisemblance) ; (ii) \(\alpha\) est en réalité contraint à \([0,\infty)\), de sorte que l’estimateur contraint est \(\hat\alpha^+=\max(\hat\alpha_{MV},0)\) ; (iii) pour un paramètre scalaire non contraint, \(2\{\ell(\hat\alpha_{MV})-\ell(0)\} \overset{\cdot}{\sim} \chi^2_1\) sous \(H_0\).

Solution

Puisque \(\hat\alpha_{MV} \overset{\cdot}{\sim} \mathcal N(0,v)\) sous \(H_0\), et que la loi normale est symétrique autour de \(0\), on a \(P(\hat\alpha_{MV}>0) = P(\hat\alpha_{MV}<0) = 1/2\) (asymptotiquement).

Cas \(\hat\alpha_{MV} \leq 0\) (probabilité \(\to 1/2\)). L’estimateur contraint est alors \(\hat\alpha^+ = 0\), c’est-à-dire que le maximum de la log-vraisemblance sur \([0,\infty)\) est atteint exactement à la frontière. La statistique \(Q = 2\{\ell(\hat\alpha^+)-\ell(0)\} = 2\{\ell(0)-\ell(0)\} = 0\) dans ce cas : \(Q\) prend la valeur \(0\) avec probabilité \(\to 1/2\), ce qui correspond à la composante \(\chi^2_0\) (loi dégénérée en \(0\)) du mélange.

Cas \(\hat\alpha_{MV} > 0\) (probabilité \(\to 1/2\)). La contrainte \(\alpha\geq 0\) n’est alors pas active : \(\hat\alpha^+ = \hat\alpha_{MV}\), et \(Q = 2\{\ell(\hat\alpha_{MV})-\ell(0)\}\) coïncide avec la statistique non contrainte, qui suit (conditionnellement à \(\hat\alpha_{MV}>0\), mais la loi \(\chi^2_1\) émerge aussi de cette moitié de la distribution normale) une loi \(\chi^2_1\) par le résultat (iii).

En combinant les deux cas, chacun de probabilité asymptotique \(1/2\), la loi de \(Q\) sous \(H_0\) est bien le mélange \[Q \ \overset{\cdot}{\sim}\ \tfrac12\chi^2_0 + \tfrac12\chi^2_1,\] c’est-à-dire une masse ponctuelle de poids \(1/2\) en \(0\), et pour le reste, une loi \(\chi^2_1\) de poids \(1/2\). Ceci explique pourquoi le seuil observé se calcule par \(\alpha^* = P(\chi^2_1>Q_{obs})/2\) plutôt que \(P(\chi^2_1>Q_{obs})\) : on ne retient que la moitié de la probabilité, celle associée au cas où la contrainte n’est pas active.

Exercice 4. Avec le lien logit, on considère une observation future dont le prédicteur linéaire estimé est \(\hat\eta^*=3\) avec \(\widehat{\text V}(\hat\eta^*)=4\). On utilise \(z_{0{,}025}=1{,}96\).

1. Construisez l’intervalle de confiance à 95% pour \(\eta^*\), puis transformez ses bornes par \(g^{-1}(\eta)=e^\eta/(1+e^\eta)\) pour obtenir un intervalle pour \(\mu^*\). Vérifiez qu’il est contenu dans \([0,1]\).

2. À l’aide de la méthode delta, montrez que \(\widehat{\text V}(\hat\mu^*) \approx [g^{-1\prime}(\hat\eta^*)]^2\,\widehat{\text V}(\hat\eta^*)\), calculez cette variance, puis construisez l’intervalle de Wald directement sur \(\hat\mu^*\). Que remarquez-vous ?

Solution

a) L’intervalle pour \(\eta^*\) est \([3-1{,}96\times 2;\ 3+1{,}96\times 2] = [-0{,}92;\ 6{,}92]\). En appliquant \(g^{-1}(\eta)=e^\eta/(1+e^\eta)\) aux deux bornes : \[g^{-1}(-0{,}92) \approx 0{,}285, \qquad g^{-1}(6{,}92) \approx 0{,}999.\] L’intervalle obtenu pour \(\mu^*\), \([0{,}285;\ 0{,}999]\), est bien contenu dans \([0,1]\) — et ce sera toujours le cas, puisque \(g^{-1}\) prend ses valeurs dans \((0,1)\) par construction, quelles que soient les bornes de l’intervalle sur \(\eta^*\).

b) La méthode delta donne, pour une transformation \(h\) dérivable, \(Var[h(\hat\eta^*)] \approx [h'(\hat\eta^*)]^2\,Var[\hat\eta^*]\) ; avec \(h=g^{-1}\), on obtient \(\widehat{\text V}(\hat\mu^*) \approx [g^{-1\prime}(\hat\eta^*)]^2\,\widehat{\text V}(\hat\eta^*)\). Puisque \(g^{-1\prime}(\eta) = \pi(1-\pi)\) avec \(\pi=g^{-1}(\eta)\) (propriété de la fonction logistique, voir chapitre 3), et \(\hat\mu^*=g^{-1}(3)\approx 0{,}953\) : \[g^{-1\prime}(3) = 0{,}953\times(1-0{,}953) \approx 0{,}0448.\] D’où \(\widehat{\text V}(\hat\mu^*) \approx (0{,}0448)^2\times 4 \approx 0{,}00803\), et \(\sqrt{\widehat{\text V}(\hat\mu^*)}\approx 0{,}0896\).

L’intervalle de Wald direct sur \(\hat\mu^*\) est alors \[[0{,}953 - 1{,}96\times 0{,}0896;\ 0{,}953+1{,}96\times 0{,}0896] \approx [0{,}777;\ 1{,}129].\]

La borne supérieure, \(1{,}129\), dépasse \(1\) : cet intervalle n’est pas valide pour une probabilité. Ceci illustre concrètement pourquoi on préfère toujours construire l’intervalle sur l’échelle du prédicteur linéaire \(\eta^*\) (sans borne) avant de transformer ses bornes par \(g^{-1}\), plutôt que d’appliquer la méthode delta directement sur \(\hat\mu^*\) : la première approche respecte automatiquement le domaine \((0,1)\) de \(\mu^*\), contrairement à la seconde.