8 Pour aller plus loin
8.1 Sélection des variables
Un modèle complet, dans le sens utilisé ici, contient toutes les variables candidates et éventuellement des interactions. Un modèle réduit n’en contient qu’un sous-ensemble.
Le “meilleur” modèle dépend de l’objectif :
- prédiction : performance hors échantillon, validation croisée, pénalisation ;
- explication : contrôle des facteurs de confusion, interprétabilité, cohérence scientifique ;
- description : modèle simple, lisible et fidèle aux données observées.
En pratique, nous avons souvent beaucoup de variables candidates, mais nous ne pouvons pas toutes les inclure dans le modèle final. Un choix qui est alors fait en épidémiologie, est de faire comme un “tri” de certaines variables, et ceci se fait en deux étapes.
8.1.1 Choix initial des variables
Ce premier tri permet de trouver des variables à exclure pour l’analyse multivariée. Pour cela, on fait des analyses univariées pour chaque variable candidate, et on ne retient que celles qui sont associées à la variable d’intérêt avec un seuil de significativité plus large que celui habituellement utilisé (p. ex., 20 % ou 25 %). On ne testera pas, ni l’exposition d’intérêt (car cela n’aurait pas de sens), ni les facteurs de confusion connus (si on en connaît).
Dans l’exemple birthweight, on étudie l’association entre le tabagisme maternel pendant la grossesse et le poids de naissance. Les variables qu’on a à disposition sont :
- l’âge maternel
AGE; - le poids maternel avant grossesse
LWT; - l’origine maternelle
Origine; - l’irritabilité utérine
UI; - l’hypertension
HT; - les visites médicales au premier trimestre, recodées en
FTV2; - les antécédents d’accouchement prématuré, recodés en
PTL2.
Le tabagisme maternel étant notre exposition d’intérêt, il ne sert à rien de la tester dans les analyses univariées. Ici, on va naïvement considérer qu’on n’a pas trouvé de facteurs de confusion connus : les variables candidates sont alors toutes celles listées ci-dessus.
Le tableau suivant résume les résultats des analyses univariées :
| Variable | Estimation | Erreur standard | p-valeur |
|---|---|---|---|
| Âge maternel (années) | 12,36 | 10,02 | 0,219 |
| Poids maternel (livres) | 4,43 | 1,71 | 0,010 |
| Origine maternelle (ref. Caucasienne) | 0,008 | ||
| Afro-américaine vs Caucasienne | -384,05 | 157,87 | 0,016 |
| Autre vs Caucasienne | -299,72 | 113,68 | 0,009 |
| Irritabilité utérine (oui vs. non) | -580,18 | 143,53 | < 0,001 |
| Hypertension (oui vs. non) | -435,56 | 215,71 | 0,045 |
| Au moins une visite au 1er trimestre (oui vs. non) | 168,58 | 105,81 | 0,113 |
| Antécédent d’accouchement prématuré (oui vs. non) | -434,17 | 142,00 | 0,003 |
Ici, avec un seuil de 25 %, toutes les variables candidates sont retenues pour l’étape suivante. Avec un seuil de 20 %, on aurait éliminé l’âge maternel.
8.1.2 Sélection manuelle
Gardons le seuil de 20 % : on enlève alors l’âge maternel. On construit maintenant notre modèle multivarié, avec toutes les variables candidates retenues à l’étape précédente. A partir de ce modèle, on va successivement retirer les variables non significatives, en vérifiant à chaque étape que le coefficient de l’exposition d’intérêt ne change pas fortement : c’est la sélection manuelle.
Classiquement :
- on estime le modèle avec toutes les variables candidates ;
- s’il y a des variables non significatives (hors exposition et facteurs de confusion connus), on retire celle qui a la p-valeur la plus élevée ;
- on réestime le modèle sans cette variable ;
- on vérifie que le coefficient de l’exposition d’intérêt n’a pas changé fortement (donc que la variation relative est inférieure à 10 % - 20 %), sinon la variable retirée est considérée comme un facteur de confusion et on la réintègre dans le modèle ;
- on recommence cela jusqu’à ce qu’on ne puisse plus retirer des variables ;
- pour toutes les variables retenues, on vérifie qu’elles ne modifient pas l’effet de l’exposition d’intérêt (donc que l’interaction n’est pas significative) ; si c’est le cas, on peut envisager de garder l’interaction dans le modèle final.
Le résultat de cette procédure est présenté dans le tableau ci-dessous :
| Terme | M1 | M2 | M3 | |||
|---|---|---|---|---|---|---|
| Estimation (ES) | p-valeur | Estimation (ES) | p-valeur | Estimation (ES) | p-valeur | |
| Intercept | 2 859 (252) | < 0,001 | 2 872 (244) | < 0,001 | 2 838 (244) | < 0,001 |
| SMOKE | -319 (108) | 0,003 | -324 (105) | 0,002 | -355 (103) | < 0,001 |
| LWT | 4,03 (1,68) | 0,017 | 4,04 (1,67) | 0,017 | 4,24 (1,68) | 0,012 |
| Origine (réf. Caucasienne) | 0,001a | < 0,001a | < 0,001a | |||
| Afro-américaine | -464 (146) | 0,002 | -466 (145) | 0,002 | -476 (146) | 0,001 |
| Autre | -330 (116) | 0,005 | -336 (112) | 0,003 | -350 (112) | 0,002 |
| UI | -488 (137) | < 0,001 | -490 (136) | < 0,001 | -524 (135) | < 0,001 |
| HT | -570 (200) | 0,005 | -574 (199) | 0,004 | -585 (200) | 0,004 |
| FTV2 | 21,4 (98,3) | 0,828 | ||||
| PTL2 | -211 (134) | 0,118 | -208 (134) | 0,120 | ||
| Variable retirée | - | FTV2 | PTL2 | |||
| Variation relative de \(\widehat{\beta}_{\mathrm{SMOKE}}\) | - | 1,6 % | 9,6 % | |||
Notes : La variation relative de \(\widehat{\beta}_{\mathrm{SMOKE}}\) est calculée par rapport au modèle précédent. ES : erreur standard.
a Test global de l’effet de l’origine maternelle, obtenu par test F partiel.
Le modèle final retenu est : \[ \begin{aligned} \mathbb{E}\!\left(\mathrm{BWT}\mid\mathrm{SMOKE},\mathrm{LWT},\mathrm{Origine},\mathrm{UI},\mathrm{HT}\right) ={}& \beta_0 +\beta_1\mathrm{SMOKE} +\beta_2\mathrm{LWT} +\beta_3\mathbb{I} \left(\mathrm{Origine}=\text{Afro-américaine}\right)\\ &+\beta_4\mathbb{I} \left(\mathrm{Origine}=\text{Autre}\right) +\beta_5\mathrm{UI} +\beta_6\mathrm{HT}. \end{aligned} \]
8.1.3 Limites
Déjà, la limite la plus flagrante : on utilise des seuils arbitraires. Dans l’exemple précédent, on a choisi 20 % pour le tri initial et 10 % pour la variation relative du coefficient de l’exposition d’intérêt. On pourrait tout aussi bien choisir 25 % et 15 %, ou 30 % et 5 %. Le modèle final serait alors différent…
Il y en a bien sûr d’autres :
- on teste beaucoup de modèles, sans tenir compte de la multiplicité des tests ;
- le modèle final dépend fortement de l’échantillon (on ne l’a pas dit, mais cette méthode ne fonctionne pas si on a des valeurs manquantes) ;
- les p-valeurs et les IC ne tiennent pas compte de la sélection (on a pu propager de l’erreur) ;
- une variable peut être non significative et pourtant importante pour contrôler la confusion ;
- ou encore une variable peut être significative mais être un médiateur ou un collider (on ne l’a pas non plus abordé dans cette UE, mais il ne faut pas ajuster sur une telle chose).
Notez qu’il existe d’autres méthodes de sélection automatique des variables. Ceci mériterait un cours entier, mais sachez que des méthodes pénalisées existent - comme le LASSO ou l’elastic-net (pareil, mais plus adapté quand on a des variables corrélées entre elles) - et sont plus communes que la sélection automatique.
Le LASSO a été proposé par Robert Tibshirani en 1996 (Tibshirani 1996), et l’elastic-net par Hui Zou et Trevor Hastie en 2005 (Zou et Hastie 2005).
8.2 Colinéarité
La colinéarité survient lorsque deux covariables portent presque la même information.
En cas de colinéarité exacte, le logiciel ne peut pas estimer séparément tous les coefficients : il en supprime généralement un ou renvoie une estimation non définie. En cas de colinéarité forte mais non exacte, le modèle peut être estimé, mais les coefficients deviennent instables et leurs erreurs standards augmentent.
8.3 Non-linéarité et transformations
Si la relation linéaire est peu plausible, plusieurs options existent :
- transformer la variable explicative (p. ex., si \(X>0\), on peut utiliser \(\log(X)\), en adaptant alors l’interprétation du coefficient) ;
- transformer la réponse si l’objectif et l’interprétation le permettent ;
- ajouter un terme quadratique \(X^2\) ;
- utiliser des splines ;
- catégoriser uniquement si les seuils ont un sens clinique clair.
Catégoriser une variable continue uniquement pour obtenir une p-valeur plus favorable est une mauvaise pratique. On perd de l’information et on crée des seuils artificiels.
8.4 Estimation des coefficients de régression
Nous ne l’avons pas présenté dans le cours, mais les coefficients \(\hat\beta\) sont obtenus par la formule : \[ \hat\beta = (\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top\mathbf y \] où \[ \mathbf X = \begin{pmatrix} 1 & x_{11} & \cdots & x_{1p}\\ \vdots & \vdots & & \vdots\\ 1 & x_{n1} & \cdots & x_{np} \end{pmatrix} \] est une matrice \(n\times(p+1)\), avec \(p\) le nombre de variables explicatives, et “+1” la colonne de 1 (correspondant à l’intercept). Pour qu’on puisse faire \((\mathbf X^\top\mathbf X)^{-1}\) (qu’on appelle “inversion”), il faut que \(\mathbf X\) soit de rang complet : il ne faut pas avoir deux colonnes identiques, ou une colonne qui est une combinaison linéaire exacte d’autres colonnes.
C’est pour cela que, si on essaie d’estimer un modèle avec “age” et “2*age” par exemple, le logiciel renverra une erreur !
Pour les férus de maths qui souhaitaient leur dose, on propose de démontrer ce résultat ci-dessous. On note : \[ y = \begin{pmatrix} y_1\\ \vdots\\ y_n \end{pmatrix} \] la matrice des valeurs observées de la variable dépendante \(Y\). En utilisant les notations précédentes, on peut réécrire le modèle de régression linéaire sous forme matricielle : \[ \mathbf y=\mathbf X\boldsymbol\beta+\boldsymbol\varepsilon, \] où \(\mathbf X\) contient une colonne de 1 pour l’intercept. On a deux manières d’obtenir les coefficients \(\hat\beta\) : par la méthode des moindres carrés, ou par la méthode du maximum de vraisemblance. Sous l’hypothèse d’erreurs normales indépendantes et homoscédastiques, les deux méthodes donnent le même estimateur de \(\boldsymbol\beta\).
8.4.1 Méthode des moindres carrés
Adrien-Marie Legendre publie une première présentation de la méthode des moindres carrés en 1805 pour déterminer des orbites de comètes (Legendre 1805). Carl Friedrich Gauss la présente à son tour en 1809 et indique l’avoir utilisée dès 1795 (Gauss 1809).
Les moindres carrés choisissent les coefficients qui minimisent la somme des carrés des résidus : \[ Q(\boldsymbol\beta) =\sum_{i=1}^n \left(y_i-\beta_0-\beta_1x_{i1}-\cdots-\beta_px_{ip}\right)^2 =(\mathbf y-\mathbf X\boldsymbol\beta)^\top (\mathbf y-\mathbf X\boldsymbol\beta) \]
Cette quantité est notée \(SCR\) dans ce cours, et souvent \(SSE\) ou \(RSS\) en anglais.
8.4.1.1 Régression linéaire simple
Ici, \(p=1\) et le modèle est \[ y_i=\beta_0+\beta_1x_i+\varepsilon_i, \; i=1,\ldots,n \] On souhaite commettre le moins d’erreur possible sur la prédiction de \(y_i\) à partir de \(x_i\). On cherche alors \((\beta_0,\beta_1)\) de façon à minimiser la fonction \(Q(\beta_0,\beta_1)\). Autrement dit : \[ (\hat\beta_0,\hat\beta_1)=\arg\min_{\beta_0,\beta_1} Q(\beta_0,\beta_1) \]
On veut donc minimiser \[ Q(\beta_0,\beta_1) = \sum_{i=1}^n(y_i-\beta_0-\beta_1x_i)^2 \] En mathématiques, on peut trouver un extremum d’une fonction en annulant son gradient (le vecteur des dérivées partielles). Ici, ses dérivées partielles sont \[ \begin{aligned} \frac{\partial Q}{\partial\beta_0} &=-2\sum_{i=1}^n(y_i-\beta_0-\beta_1x_i),\\ \frac{\partial Q}{\partial\beta_1} &=-2\sum_{i=1}^n x_i(y_i-\beta_0-\beta_1x_i) \end{aligned} \] donc le gradient est \[ \nabla Q(\beta_0,\beta_1) =\begin{pmatrix} \frac{\partial Q}{\partial\beta_0} \\ \frac{\partial Q}{\partial\beta_1} \end{pmatrix} \]
A l’extremum, ces deux dérivées doivent être nulles. Ceci se traduit par le système d’équations suivant, appelé équations normales : \[ \sum_{i=1}^n y_i-n\widehat\beta_0-\widehat\beta_1\sum_{i=1}^n x_i=0 \iff \bar y-\widehat\beta_0-\widehat\beta_1\bar x=0 \iff \widehat\beta_0=\bar y-\widehat\beta_1\bar x \tag{1}\]
Une conséquence de ce résultat est que la droite des moindres carrés passe donc par le point \((\bar x,\bar y)\).
La deuxième équation normale s’écrit \[ \sum_{i=1}^n x_i y_i -\widehat\beta_0\sum_{i=1}^n x_i -\widehat\beta_1\sum_{i=1}^n x_i^2=0 \] En injectant \(\widehat\beta_0\) par l’estimation de l’équation 1, on a : \[ \begin{aligned} & 0= \sum_{i=1}^n x_i y_i - (\bar y-\widehat\beta_1\bar x)\sum_{i=1}^n x_i - \widehat\beta_1\sum_{i=1}^n x_i^2 \\ &= \sum_{i=1}^n x_i y_i-n\bar x\bar y +\widehat\beta_1n\bar x^2-\widehat\beta_1\sum_{i=1}^n x_i^2 \\ \implies \widehat\beta_1 \left(\sum_{i=1}^n x_i^2-n\bar x^2\right) &=\sum_{i=1}^n x_i y_i-n\bar x\bar y \end{aligned} \] Par conséquent, \[ \boxed{ \widehat\beta_1 =\frac{\sum_{i=1}^n(x_i-\bar x)(y_i-\bar y)} {\sum_{i=1}^n(x_i-\bar x)^2}}, \] et, en utilisant l’équation 1, \[ \boxed{\widehat\beta_0=\bar y-\widehat\beta_1\bar x} \]
Il se trouve que, comme
\[ \sum_{i=1}^n(x_i-\bar x)^2 = \sum_{i=1}^n x_i^2-n\bar x^2 \]
et
\[ \sum_{i=1}^n(x_i-\bar x)(y_i-\bar y) = \sum_{i=1}^n x_i y_i-n\bar x\bar y \]
alors
\[ \widehat\beta_1 =\frac{\sum_{i=1}^n(x_i-\bar x)(y_i-\bar y)}{\sum_{i=1}^n(x_i-\bar x)^2} =\frac{\frac1n\sum_{i=1}^n(x_i-\bar x)(y_i-\bar y)}{\frac1n\sum_{i=1}^n(x_i-\bar x)^2} =\frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}(X)} \]
En mathématiques, pour vérifier qu’un extremum est bien un minimum, on peut regarder la matrice hessienne (la matrice des dérivées secondes). Si elle est définie positive - c’est-à-dire que tous ses mineurs principaux sont positifs (critère de Sylvester) - alors on a un minimum.
La matrice hessienne de \(Q\) est \[ \mathbf H = \begin{pmatrix} \frac{\partial^2 Q}{\partial\beta_0^2} & \frac{\partial^2 Q}{\partial\beta_0\partial\beta_1} \\ \frac{\partial^2 Q}{\partial\beta_1\partial\beta_0} & \frac{\partial^2 Q}{\partial\beta_1^2} \end{pmatrix} = \begin{pmatrix} 2n & 2\sum_i x_i\\ 2\sum_i x_i & 2\sum_i x_i^2 \end{pmatrix} =2\begin{pmatrix} n & \sum_i x_i\\ \sum_i x_i & \sum_i x_i^2 \end{pmatrix} \] Le premier mineur principal dominant correspond à la sous-matrice \(1\times 1\) en haut à gauche, soit \(2n>0\) (car \(n>0\)). Le second mineur principal dominant correspond à la matrice complète, et son déterminant vaut \[ \left(2n\times2\sum_i x_i^2-2\sum_i x_i \times 2\sum_i x_i\right) = 4n\sum_i(x_i-\bar x)^2 \] qui est strictement positif si les \(x_i\) ne sont pas tous égaux. Ainsi, on a bien un minimum.
8.4.1.2 Régression linéaire multiple
C’est globalement le raisonnement précédent, mais avec plus de variables explicatives. Le plus simple est alors de travailler avec la notation matricielle de \(Q\). En développant le produit scalaire, on obtient \[ \begin{aligned} Q(\boldsymbol\beta) &=(\mathbf y-\mathbf X\boldsymbol\beta)^\top (\mathbf y-\mathbf X\boldsymbol\beta)\\ &=\mathbf y^\top\mathbf y -2\boldsymbol\beta^\top\mathbf X^\top\mathbf y + \boldsymbol\beta^\top\mathbf X^\top\mathbf X\boldsymbol\beta \end{aligned} \] Le gradient par rapport à \(\boldsymbol\beta\) est \[ \nabla_{\boldsymbol\beta}Q = -2\mathbf X^\top\mathbf y +2\mathbf X^\top\mathbf X\boldsymbol\beta \] En l’annulant, on obtient les équations normales : \[ \mathbf X^\top\mathbf X\widehat{\boldsymbol\beta} = \mathbf X^\top\mathbf y \] Si \(\mathbf X\) est de rang colonne complet, alors \(\mathbf X^\top\mathbf X\) est inversible et \[ \boxed{ \widehat{\boldsymbol\beta} =(\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top\mathbf y.} \] La hessienne vaut \[ \nabla^2_{\boldsymbol\beta}Q=2\mathbf X^\top\mathbf X \] Elle est définie positive lorsque \(\mathbf X\) est de rang colonne complet, car, pour tout vecteur non nul \(\mathbf a\), \[ \mathbf a^\top\mathbf X^\top\mathbf X\mathbf a =\|\mathbf X\mathbf a\|^2>0 \] L’estimateur est donc l’unique minimum global !
8.4.2 Maximum de vraisemblance
L’autre possibilité consiste à travailler avec la vraisemblance, qui correspond au produit des densités de probabilité des observations si elles sont indépendantes.
En 1922, Ronald A. Fisher développe la méthode du maximum de vraisemblance dans un cadre général d’estimation statistique (Fisher 1922).
Pour rappel, on suppose que les erreurs \(\varepsilon_i\) sont indépendantes et identiquement distribuées selon une loi normale centrée réduite, de variance \(\sigma^2\) inconnue : \[ \varepsilon_i\mid\mathbf X \overset{\mathrm{iid}}{\sim}\mathcal N(0,\sigma^2) \] Ceci permet d’écrire, conditionnellement à \(\mathbf X\), \[ Y_i\sim\mathcal N(\mathbf x_i^\top\boldsymbol\beta,\sigma^2) \]
La densité d’une observation \(y_i\) est donc \[ f(y_i\mid\mathbf x_i\, ; \,\boldsymbol\beta,\sigma^2) =\frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left[-\frac{(y_i-\mathbf x_i^\top\boldsymbol\beta)^2} {2\sigma^2}\right] \] Comme les \(Y_i\) sont conditionnellement indépendants, la vraisemblance est le produit des densités de probabilité de chaque \(Y_i\) : \[ \begin{aligned} L(\boldsymbol\beta,\sigma^2) &=\prod_{i=1}^n f(y_i\mid\mathbf x_i\, ; \,\boldsymbol\beta,\sigma^2)\\ &=(2\pi\sigma^2)^{-n/2} \exp\!\left[-\frac{1}{2\sigma^2} \sum_{i=1}^n(y_i-\mathbf x_i^\top\boldsymbol\beta)^2\right]\\ &=(2\pi\sigma^2)^{-n/2} \exp\!\left[-\frac{Q(\boldsymbol\beta)}{2\sigma^2}\right] \end{aligned} \]
Il est difficile de travailler avec cette expression. On préfère donc travailler avec la log-vraisemblance, qui est la fonction logarithme de la vraisemblance. Comme le logarithme est strictement croissant, maximiser la log-vraisemblance revient à maximiser la vraisemblance.
C’est cette quantité qu’on utilise pour le test du rapport de vraisemblance !
La log-vraisemblance vaut \[ \ell(\boldsymbol\beta,\sigma^2)= -\frac n2\log(2\pi) - \frac n2\log(\sigma^2) - \frac{Q(\boldsymbol\beta)}{2\sigma^2} \]
Comme nous souhaitons maximiser la (log-)vraisemblance, on peut chercher les quantités \((\widehat{\boldsymbol\beta},\widehat\sigma^2)\) qui annulent le gradient de la log-vraisemblance.
En dérivant par rapport à \(\boldsymbol\beta\), on obtient \[ \frac{\partial\ell}{\partial\boldsymbol\beta} =\frac{1}{\sigma^2}\mathbf X^\top (\mathbf y-\mathbf X\boldsymbol\beta) = 0 \iff \mathbf X^\top\mathbf X\widehat{\boldsymbol\beta} =\mathbf X^\top\mathbf y \implies \hat{\boldsymbol\beta}=(\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top\mathbf y \]
Similairement, en dérivant par rapport à \(\sigma^2\), on obtient \[ \frac{\partial\ell}{\partial\sigma^2} =-\frac{n}{2\sigma^2} +\frac{Q(\boldsymbol\beta)}{2\sigma^4} =0 \iff \boxed{\widehat\sigma^2_{\mathrm{MV}} =\frac{Q(\widehat{\boldsymbol\beta})}{n}} \]
CQFD !