Cours de Régression Logistique Appliquée
Patrick Taffé, PhD Institut Universitaire de Médecine Sociale et Préventive (IUMSP) et Centre d’épidémiologie Clinique (CepiC) Lausanne, Août 2004
i
Table des matières Introduction .............................................................................................................. 1 Pourquoi la statistique ? .......................................................................................... 1 Pourquoi la régression logistique ? ......................................................................... 1 1) La modélisation d’une variable qualitative dichotomique................................ 3 Exercice 1................................................................................................................ 5 2) Formulation mathématique du modèle de régression logistique (*) ............... 7 2.1) Le modèle de régression linéaire Normal......................................................... 7 2.2) Le modèle de régression logistique.................................................................. 8 2.3) Y-a-t’il d’autres modèles !? ............................................................................... 9 Exercice 2.............................................................................................................. 10 3) Estimation et tests (*)......................................................................................... 13 3.1) L’estimation du modèle .................................................................................. 13 3.2) Test de significativité des coefficients ............................................................ 13 Exercice 3.............................................................................................................. 14 4) La transformation logit ...................................................................................... 17 Exercice 4.............................................................................................................. 18 5) Le succès du modèle Logit : l’Odds Ratio ....................................................... 21 5.1) L’Odds Ratio comme mesure d’association ................................................... 21 5.2) L’Odds Ratio comme mesure du risque relatif (RR)....................................... 22 Exercice 5.............................................................................................................. 24 6) L’interprétation des coefficients ....................................................................... 27 6.1) Le cas d’un modèle additif, i.e. sans interactions ........................................... 27 a) La constante du modèle ................................................................................. 28 b) Coefficient d’une variable explicative dichotomique ....................................... 29 c) Coefficient d’une variable explicative polytomique ......................................... 30 d) Coefficient d’une variable explicative continue............................................... 31 e) L’Odds ratio associé à la variation de plusieurs co-variables.......................... 32 6.2) Le cas d’un modèle non additif, i.e. avec interactions .................................... 32 Exercice 6.............................................................................................................. 34 7) Stratégie de modélisation.................................................................................. 39 Pourquoi construire un modèle ?........................................................................... 39 Existe-t-il une stratégie de modélisation conduisant à un « bon » modèle ?.......... 39
ii
7.1) Le choix des co-variables............................................................................... 40 7.2) Le choix de la forme fonctionnelle des co-variables ...................................... 40 7.3) L’adéquation du modèle aux données « Goodness of fit » (*) ....................... 41 a) La notion de « covariate pattern » ................................................................. 42 b) Evaluation de la calibration du modèle : le test de Hosmer et Lemeshow ..... 42 c) L’analyse des résidus..................................................................................... 43 c.1) Le résidu de Pearson............................................................................................... 44 c.2) Le résidu de déviance ............................................................................................. 46 d) Détection des « covariate patterns » mal ajustés .......................................... 47 e) Détection des points influants (effet de levier) ............................................... 48 f) Evaluation du pouvoir discriminant du modèle : sensibilité, spécificité et courbe ROC ................................................................................................................... 49 g) La validation du modèle ................................................................................. 51 7.4) Limitations et biais (*)..................................................................................... 52 a) Le problème de la séparabilité ou quasi-séparabilité (*) ................................ 52 b) Le problème de « l’overfitting » ...................................................................... 53 c) Le biais de sélection....................................................................................... 53 d) Le problème de surdispersion « overdispersion ».......................................... 54 e) Extensions ..................................................................................................... 54 e.1) Le cas de données répétées..................................................................................... 54 e.2) Le cas de données agrégées « cluster » .................................................................. 54 Exercice 7 ............................................................................................................. 54 8) Le logiciel statistique STATA............................................................................ 55 Bibliographie .......................................................................................................... 59 Livres :................................................................................................................... 59 Articles: ................................................................................................................. 59 Pour l’utilisation de STATA se référer aux manuels suivants :.............................. 60
iii
Avant propos Ce cours a pour but d’introduire le lecteur à la problématique de la modélisation des variables qualitatives dichotomiques (i.e. comportant deux catégories comme « sain » et « malade ») au moyen de la régression logistique. L’analyse de régression logistique est plus complexe que celle de régression linéaire, car le modèle logistique est non-linéaire. Nous allons, autant que possible, faire un parallèle entre les deux types d’analyses et illustrer les différences fondamentales. Il s’agit d’un cours de régression logistique appliquée de sorte que nous n’insisterons pas sur les détails mathématiques, mais plutôt sur les concepts fondamentaux. Néanmoins, la statistique est avant tout une discipline faisant appel aux mathématiques et même si le programme statistique prend en charge tous les aspects formels, un minimum de formalisme est nécessaire pour bien illustrer les concepts. Nous avons donc décidé de ne pas occulter complètement les mathématiques de ce cours et les sections d’un caractère plus technique seront indiquées par un astérisque « * ». Les données pour les exercices peuvent être téléchargées depuis le web aux adresses : ftp://ftp.wiley.com/public/sci_tech_med/logistic/ http://www.ats.ucla.edu/stat/stata/examples/alr2/default.htm
iv
Introduction Le but de ce cours est d’exposer les fondements de la régression logistique de manière intuitive et aussi peu formelle que possible, et d’illustrer les étapes de la modélisation des variables qualitatives binaires.
Pourquoi la statistique ? En général, le but de la plupart des recherches est de déterminer des relations entre un ensemble de variables. Les techniques « multivariables » ont été développées à cette fin. Souvent on considère une variable dépendante que l’on veut prédire et des variables indépendantes ou explicatives. Remarquons que bien souvent le terme « multivarié » est confondu avec « multivariables », ce qui peut porter à confusion étant donné que le premier se réfère à la situation où l’on considère plusieurs variables dépendantes à la fois, tandis que le deuxième plus vague correspond peutêtre mieux à la situation la plus fréquente en épidémiologie où l’on considère une seule variable dépendante et plusieurs variables explicatives. Il est difficile de donner une définition consensuelle de la statistique, mais certainement cette discipline traite de l’incertitude, de la variabilité, de l’inférence (test d’hypothèses, intervalles de confiance, prédiction, …). On retiendra qu’elle a pour but de quantifier un phénomène d’intérêt et d’apporter une information concernant la précision avec laquelle les résultats ont été établis. Par exemple, pour estimer la taille moyenne des jeunes de 15 ans en Suisse on considère un échantillon d’élèves dans une école et l’on calcule leur taille moyenne. Cette estimation ne sera certainement pas parfaite puisqu’elle repose sur un petit collectif dont on espère qu’il soit suffisamment représentatif de l’ensemble de cette population en Suisse. Un intervalle de confiance nous permettra d’apprécier le degré d’incertitude de notre évaluation. L’analyse de régression est une technique statistique permettant d’établir une relation entre une variable dépendante et des variables explicatives, afin d’étudier les associations et de faire des prévisions. On peut, par exemple, s’intéresser à quantifier la relation entre le risque de décès et la quantité de cigarettes fumées quotidiennement, tout en ajustant pour l’âge, le sexe, et éventuellement d’autres facteurs de risque.
Pourquoi la régression logistique ? Lorsque la variable dépendante n’est pas quantitative mais qualitative ou catégorielle le modèle de régression linéaire n’est pas approprié. Ce qui distingue le modèle de régression logistique du modèle de régression linéaire est que dans le premier la variable dépendante est qualitative, i.e. cette variable prend comme valeur un attribut et non pas une valeur numérique : par exemple la variable état de santé prend les attributs « sain » ou « malade », la variable sexe « mâle » ou « femelle », une autre variable les attributs « rouge » ou « noir », etc. Lorsque le nombre d’attributs est deux l’on parle de variable dichotomique, e.g. le sexe « mâle » ou « femelle », tandis que s’il est supérieur à deux l’on a une variable polytomique, e.g. une pression « haute », « normale » ou « basse ».
1
Dans le modèle de régression linéaire la variable dépendante est, en revanche, quantitative, car elle admet une échelle de mesure naturelle : par exemple la pression systolique 50-200 mmHg, le poids 30-200 kg, la taille 1-2 m, le niveau de CD4 0-2000 cell/ìL, etc. Lorsque la variable dépendante est quantitative l’hypothèse de normalité de la distribution de cette variable ou d’une transformation est généralement plausible, tandis que lorsqu’elle est qualitative elle n’admet pas de valeur numérique naturelle (puisqu’elle ne peut prendre que des attributs) et le modèle normal n’est pas approprié. Une variable aléatoire qualitative est décrite par les probabilités des différents attributs qu’elle peut prendre et pour évaluer l’influence de différents facteurs sur cette variable il est d’usage de modéliser les probabilités des différents attributs. Un modèle décrivant la probabilité avec laquelle la variable qualitative dichotomique sexe prend les attributs « femelle » ou « mâle » est le modèle « binomial » (avec n = 11). Lorsque le nombre d’attributs que peu prendre cette variable est supérieur à deux on a une variable polytomique et un modèle décrivant cette situation est le modèle « multinomial ». On a représenté, ci-dessous, différents graphes illustrant les différences fondamentales entre variable qualitative et variable quantitative. Dans le premier graphe la variable dépendante est la maladie coronarienne. Cette variable peut prendre les attributs « oui » ou « non » de sorte qu’il n’est pas possible d’écrire une relation directement entre la maladie coronarienne et l’âge. Dans le second graphe la variable dépendante est quantitative, il s’agit de la taille, de sorte qu’il est possible d’établir directement une relation (linéaire ou pas) entre la taille et l’âge. Le troisième graphe illustre l’hypothèse de Normalité souvent adoptée lorsque la variable dépendante est quantitative. Maladie coronarienne
Relation entre taille et âge chez les enfants taille
oui
non âge
âge
Relation entre taille et âge chez les enfants: hypothèse de Normalité taille
âge
figures 1 à 3 1
Lorsque n=1 le modèle binomial se réduit au modèle de Bernoulli.
2
1) La modélisation d’une variable qualitative dichotomique Nous avons vu que lorsque la variable dépendante était qualitative elle n’admettait pas d’échelle de mesure naturelle et que l’on modélisait, par conséquent, sa probabilité de prendre tel ou tel attribut. Voyons comment cela s’applique dans notre exemple de maladie coronarienne en fonction de l’âge. Dans le graphique suivant l’on a regroupé les données concernant l’âge en catégories et calculé dans chacune de ces catégories le pourcentage de personnes souffrant d’une maladie coronarienne :
Pourcentage de personnes souffrant d’une maladie coronarienne par catégorie d’âge 1
0.5
0 âge
figure 4
On constate que l’on a une relation sigmoïdale, i.e. en forme de S, entre la proportion de maladie coronarienne et l’âge. On en déduit, ainsi, que pour modéliser la probabilité de maladie coronarienne en fonction de l’âge il faudra utiliser une relation sigmoïdale. En effet, une probabilité étant par définition comprise entre 0 et 1 le modèle linéaire n’est bien entendu pas approprié (puisqu’il ne limite pas les valeurs de notre probabilité au domaine compris entre 0 et 1) et la relation est forcément non-linéaire :
Pourcentage de personnes souffrant d’une maladie coronarienne par catégorie d’âge: relation linéaire
Pourcentage de personnes souffrant d’une maladie coronarienne par catégorie d’âge: relation non linéaire (sigmoïdale)
>1 1
1
0.5
0.5
0
0 <0
âge
âge
figures 5 & 6
3
Remarquons qu’une probabilité est une caractéristique d’une population, tandis qu’une proportion est calculée à partir d’un échantillon. Cette dernière s’approche d’autant plus de la probabilité (inconnue en général) que l’échantillon est grand. Un choix intuitif pour modéliser une probabilité est d’utiliser une fonction de répartition ou fonction cumulative. Illustrons ce point avec l’exemple des fonctions de répartition des lois Normale et Logistique. Pour rappeler la différence, nous illustrons aussi les fonctions de densité correspondantes :
Fonctions de densité de diverses lois Normales
0
0
.2
.1
.4
.2
.6
.3
.8
1
.4
Fonctions cumulatives de diverses lois Normales
-10
-5
0 x N(0,1) N(2,3)
5
10
-10
-5
N(0,3)
0 x N(0,1) N(2,3)
10
Fonctions de densité de diverses lois Logistiques
.4 .3 .2 .1 0
0
.2
.4
.6
.8
1
.5
Fonctions cumulatives de diverses lois Logistiques
5 N(0,3)
-10
-5
0 x Logistique(0,1) Logistique(2,3)
5
10
-10
Logistique(0,3)
-5
0 x Logistique(0,1) Logistique(2,3)
5
10
Logistique(0,3)
figures 7 à 10
On constate qu’en fonction de la moyenne (1er paramètre) les courbes se déplacent le long de l’abscisse et qu’en fonction de la variance (2e paramètre) la pente de la fonction de répartition change. On peut déjà anticiper que se seront ces deux paramètres qu’il faudra estimer (au moyen d’une technique statistique) pour obtenir un bon ajustement de notre courbe aux données. Remarque (*)
(
)
La fonction cumulative d’une loi Normale de moyenne µ et de variance σ 2 , i.e. N µ ,σ 2 , s’écrit :
F ( x) = ∫
x
−∞
1 t − µ 2 dt exp − 2 σ 2π σ 1
4
tandis que celle d’une loi Logistique de moyenne µ et de variance σ 2 s’écrit : π x−µ exp 3 σ F ( x) = π x−µ 1 + exp 3 σ En résumé, nous avons donc dit que lorsque la variable dépendante était qualitative l’on modélisait la probabilité de ses attributs, qu’un modèle mathématique adéquat avait une forme sigmoïdale comme une fonction de répartition et que la forme de cette sigmoïde changeait en fonction des paramètres caractérisant cette fonction de répartition. Il s’agit, ensuite, d’établir un lien entre ces paramètres (donc la forme et la position de notre courbe sigmoïdale), la probabilité de maladie coronarienne (la variable dépendante d’intérêt) et l’âge (la variable explicative). Pour cela, dans le prochain chapitre, nous allons formuler un modèle de régression (non-linéaire). Un modèle très utilisé en épidémiologie est le modèle Logistique.
Exercice 1 Le but de cet exercice est d’illustrer la différence fondamentale entre variable qualitative et variable quantitative. Nous allons montrer, en particulier, qu’il n’est pas approprié de traiter une variable qualitative comme si elle était quantitative : par exemple de régresser directement une variable dépendante qualitative codée « 0 » et « 1 » en fonction d’une variable explicative, comme on le fait en régression linéaire. A cette fin, nous allons utiliser des données rapportant la présence ou l’absence d’une maladie coronarienne. 1) Représentation graphique des données Dans ce fichier de données nous avons une seule variable explicative l’âge. Afin de « visualiser » la relation entre la présence ou l’absence d’une maladie coronarienne en fonction de l’âge nous allons représenter les données sur un graphe. Pour cela l’option « jitter(2) » de STATA s’avère utile. Essayez la commande suivante avec et sans cette option : scatter chd age, jitter(2) ylabel(0(1)1) coronarienne en fonction de l'âge)
ytitle(chd
0/1)
title(Maladie
On constate que plus on est âgé plus le risque d’avoir un problème coronarien semble élevé. 2) Un exemple à ne pas suivre : estimation d'une relation linéaire entre la variable dépendante dichotomique chd et l’age La régression linéaire de la variable chd en fonction de la variable age fournit une droite n’ayant pas de sens, car chd ne peut prendre que deux valeurs 0 ou 1 tandis que la droite de régression linéaire prédit des valeurs impossibles. * régression linéaire
5
regress chd age cap drop fit predict fit, xb * graphe de la relation linéaire entre les variables chd et age scatter chd age, jitter(2) ylabel(0(1)1) ytitle(chd 0/1) title(Maladie coronarienne en fonction de l'âge , size(medium)) subtitle(régression linéaire) || scatter fit age, c(l) sort saving(g1, replace)
3) Relation fonctionnelle entre la probabilité de maladie coronarienne et l'âge Pour représenter la relation fonctionnelle entre la probabilité de maladie coronarienne et l'âge nous allons définir des catégories d’âge et calculer le pourcentage de maladie coronarienne dans chacune de ces catégories : * calcul des percentiles de la variable age centile age, centile(10 20 30 40 50 60 70 80 90) * génération des percentiles de la variable age cap drop pct_age xtile pct_age=age, nquantiles(9) tab pct_age * calcul des proportions de maladie coronarienne dans les catégories d'âge sort pct_age cap drop p_chd by pct_age: egen p_chd=mean(chd) * graphe de la relation entre les catégories d'âge et la proportion de maladie coronarienne scatter p_chd pct_age, ylabel(0(0.2)1) ytitle(P(chd)) title(Proportion de maladie coronarienne en fonction de la catégorie d'âge, size(medium)) saving(g2, replace) * un graphe plus joli en utilisant une régression non paramétrique twoway scatter chd age, jitter(2) ylabel(0(0.2)1) ytitle(P(chd)) title(Proportion de maladie coronarienne en fonction de l'âge, size(medium)) || lowess chd age, sort legend(off) saving(g3, replace) graph combine g2.gph g3.gph, iscale(.55)
6
2) Formulation mathématique du modèle de régression logistique (*) Dans « modèle de régression logistique » nous avons les termes « régression » et « logistique ». Dans cette section, nous allons en illustrer la raison. Ceci nous permettra de bien comprendre les différences fondamentales entre les modèles de régression linéaire et logistique. Néanmoins, d’emblée on peut remarquer que le terme « régression » impliquera qu’on considérera un ensemble de variables explicatives et que le terme « logistique » fera référence à une hypothèse de distribution (du même nom).
2.1) Le modèle de régression linéaire Normal En statistique, le terme de « régression » de « y » par rapport à « x » fait référence à l’espérance mathématique de y conditionnelle à x, E ( y | x ) . Concrètement, cette espérance mathématique établit une relation entre x et y : connaissant la valeur prise par la variable x on prédira que y prendra en moyenne la valeur E ( y | x ) . Par exemple, E ( y | x = 5) = 33 veut dire que lorsque x vaut 5 la valeur espérée (attendue, moyenne) de y vaut 33. En régression linéaire l’on modélise l’espérance mathématique de y conditionnelle à x au moyen d’une équation linéaire : E ( y | x, β 0 , β 1 ) = β 0 + β 1 x
et le modèle s’écrit : y = β 0 + β1 x + ε où å est un résidu que l’on suppose d’espérance nulle E (ε ) = 0 et de variance constante ou homoscédastique V (ε ) = σ 2 . Souvent, on fait aussi l’hypothèse de normalité du résidu å, on dit que l’on adopte le modèle « Normal » ou « Gaussien », afin de procéder à des tests sur les paramètres β 0 et β 1 :
(
)
ε ≅ N 0, σ 2 .
On parle, alors, du modèle linéaire classique.
7
2.2) Le modèle de régression logistique Nous allons voir qu’en régression logistique l’on modélise aussi l’espérance mathématique de y conditionnelle à x, mais cette fois la relation est non-linéaire et les résidus ne peuvent pas être distribués « Normalement ». Rappelons que lorsque la variable dépendante était qualitative elle n’admettait pas de valeur numérique naturelle. On peut, néanmoins, introduire un codage quantitatif permettant de représenter les différents attributs. Par exemple, on codera « 1 » si l’attribut est « sain » et « 0 » sinon. A partir de ce codage quantitatif, on établit un lien entre l’espérance mathématique de y conditionnelle à x et la probabilité de y : 1 y= 0
(i.e." sain" ) avec probabilité
P = F ( x, β 0 , β 1 )
(i.e." malade" ) avec probabilité 1 − P = 1 − F ( x, β 0 , β 1 )
L’espérance mathématique de y conditionnelle à x (i.e. la régression de y par rapport à x), s’écrit : E ( y | x, β 0 , β 1 ) = 1 × P + 0 × (1 − P) = P = F ( x, β 0 , β 1 ) En ayant adopté le codage 0/1 la probabilité de y correspond à son espérance conditionnelle. Cette relation justifie l’utilisation du terme « régression » logistique. Il nous reste à expliquer la raison du terme « logistique ». Nous avons vu qu’un choix intuitif pour modéliser une probabilité était d’utiliser une fonction de répartition. Lorsque cette fonction de répartition est celle de la loi Logistique on obtient le modèle de régression logistique ou plus simplement le modèle Logit. Remarques : 1) Le codage en 0/1 est arbitraire mais n’a aucune influence sur les résultats des estimations, car la vraisemblance s’exprime en fonction des probabilités P et pas de l’espérance conditionnelle E ( y | x, β 0 , β 1 ) . 2) On peut écrire le modèle de régression logistique sous la même forme que le modèle de régression linéaire : y = F ( x, β 0 , β 1 ) + ε Cependant, cette fois le modèle est non-linéaire et le résidu å ne peut pas être distribué selon une loi Normale.
8
En effet, il ne peut prendre que deux valeurs ε = 1 − F ( x, β 0 , β 1 ) si
y = 1 ou
ε = − F ( x, β 0 , β 1 ) si y = 0 . De plus, sa variance n’est pas σ 2 mais V (ε ) = F ( x, β 0 , β 1 ) [1 − F ( x, β 0 , β 1 )] . On constate que la variance dépend de la variable x et, par conséquent, elle n’est pas constante mais hétéroscédastique. Formellement, appliqué à notre exemple de la maladie coronarienne le modèle Logit s’écrit :
P(maladie coronarienne | âge) = F (âge β 0 , β 1 ) =
exp( β 0 + β 1 âge) 1 + exp( β 0 + β 1âge)
Remarque : Dans cette expression la probabilité de maladie coronarienne est modélisée au moyen de la fonction de répartition d’une loi Logistique d’espérance µ = − β 0 / β 1 et d’écarttype σ = π /( 3β 1 ) . En définitive En définitive on notera que le modèle de régression logistique se distingue du modèle de régression linéaire de part 1) la distribution de la variable dépendante n’est pas Normale mais Binomiale 2) le modèle de régression est non-linéaire 3) la variance est hétéroscédastique.
2.3) Y-a-t’il d’autres modèles !? Nous avons vu qu’un choix intuitif pour modéliser une probabilité était d’utiliser une fonction de répartition. Il en existe, bien évidemment, un choix quasiment infini. Pour des raisons historiques (existence d’une tabulation, simplicité, ...) ce choix s’est porté souvent sur les fonctions de répartition des lois Normale et Logistique, la première conduisant à un modèle appelé Probit et la deuxième comme on l’a vu au modèle Logit. Ainsi, si l’on choisi la fonction de répartition de la loi Normale pour modéliser notre probabilité l’on obtient le modèle Probit :
P(maladie coronarienne | âge) = F (âge β 0 , β 1 ) = ∫
β 0 + β1âge
−∞
t 1 exp(− 2 ) dt 2 2π
Remarque : Dans cette expression la probabilité de maladie coronarienne est modélisée au moyen de la fonction de répartition d’une loi Normale d’espérance µ = − β 0 / β 1 et d’écart-type σ = 1 / β1 .
9
Les lois Normale et Logistique se distinguent, en particulier, en fonction de l’épaisseur de la queue de probabilité de la fonction de densité correspondante, ce qui a une influence sur la « vitesse » avec laquelle la fonction de répartition s’éloigne de 0 ou s’approche de 1 :
Fonctions de densité des lois Normale(0,1) et Logistique(0,1)
.4 .3 .2 .1 0
0
.2
.4
.6
.8
1
.5
Fonctions cumulatives des lois Normale(0,1) et Logistique(0,1)
-5
0 x N(0,1)
5
-5
Logistique(0,1)
0 x N(0,1)
5 Logistique(0,1)
figures 11 & 12
Néanmoins, comme on le constate sur ces figures, la différence entre les deux modèles est infime de sorte qu’en pratique l’on peut choisir indifféremment l’une ou l’autre des lois. Toutefois le modèle Logit permet une interprétation plus habituelle en épidémiologie car elle fait intervenir des Odds Ratio. Remarquons que ce résultats est valable uniquement dans le cas de la modélisation d’une variable qualitative dichotomique et que dans le cas polytomique la différence est importante.
Exercice 2 Dans cet exercice nous allons estimer une relation sigmoïdale entre les variables chd et age au moyen d’un modèle de régression logistique. Nous comparerons cette estimation avec celle fournie par un modèle Probit. On en conclura que la différence entre les deux modèles est, ici, infime. Le modèle de régression logistique est très utilisé, surtout en épidémiologie, principalement à cause de l’interprétation du coefficient d’une co-variable comme le logarithme de son Odds Ratio. Autrement dit, l’exponentiel du coefficient d’une co-variable correspond à un Odds Ratio. 1) Estimation d'une relation sigmoïdale entre les variables chd et age Pour cela nous allons utiliser la commande « logistic » de STATA. * régression logistique logistic chd age * calcul des probabilités estimées cap drop p predict p
10
* graphe de la relation sigmoïdale entre chd et age scatter chd age, jitter(2) ylabel(0(1)1) ytitle(P(chd)) title(Maladie coronarienne en fonction /// de l'âge, size(medium)) subtitle(régression logistique) || scatter p age, c(l) sort saving(g4, replace)
2) Comparaison des modèles linéaire et logistique graph combine g1.gph g4.gph, iscale(.75)
On vérifie sur ces graphes que le modèle logistique fournit une probabilité estimée de maladie coronarienne comprise entre 0 et 1, tandis que la régression linéaire fournit des valeurs aberrantes de la variable chd. 3) Estimation d'un modèle probit de la relation entre chd et age probit chd age cap drop p_probit predict p_probit * graphe de la relation entre chd et age scatter chd age, jitter(2) ylabel(0(1)1) ytitle(P(chd)) title(Maladie coronarienne en fonction /// de l'âge, size(medium)) subtitle(régression probit) || scatter p age, c(l) sort saving(g5, replace) * comparaison des modèles logit et probit graph combine g4.gph g5.gph, iscale(.75)
En conclure que, dans le cas dichotomique, la différence entre les modèles de régression logistique et probit est infime.
11
12
3) Estimation et tests (*) 3.1) L’estimation du modèle L’estimation du modèle de régression logistique se fait généralement par la méthode du maximum de vraisemblance. Pour cela on écrit la vraisemblance de l’échantillon. Lorsque les observations individuelles yi, i=1,…,n, sont supposées indépendantes, cette vraisemblance s’écrit comme le produit des probabilités : L( β 0 , β 1 ) = ∏ [P( y = 1 x, β 0 , β 1 )] i [1 − P( y = 1 x, β 0 , β 1 )] n
y
1− yi
i =1
Ensuite, on maximise cette vraisemblance par rapport aux paramètres β 0 , β 1 au moyen d’un algorithme numérique (par ex. une méthode de gradient). Remarques : 1) Quand on fait l’hypothèse d’indépendance des observations on entend qu’elles sont conditionnellement indépendantes. C’est-à-dire que les probabilités individuelles sont supposées indépendantes après ajustement pour les facteurs de risques. Ainsi, deux individus présentant les mêmes facteurs de risque ne sont pas indépendants, mais conditionnellement à ces facteurs on suppose qu’il le sont. Autrement dit, une fois que l’on a ajusté pour l’effet des différents facteurs de risque les observations peuvent être considérées comme indépendantes (mathématiquement ε | x ≈ iid (0, σ 2 ) ). 2) Lorsqu’on est en présence de mesures répétées pour chaque individu ou que les données présentent une « structure hiérarchiques », comme c’est le cas lorsqu’on échantillonne des familles et que l’on s’intéresse aux caractéristiques des membres de ces familles, l’hypothèse d’indépendance des données n’est pas plausible. En effet, les mesures répétées d’un même individu ou des membres d’une même famille sont plus semblables qu’entre individus ou familles. Dans ce cas, il faut utiliser d’autres méthodes qui prennent en compte la corrélation des données (ex : modèle marginal avec GEE, modèle logistique conditionnel, modèle mixte).
3.2) Test de significativité des coefficients Pour tester la significativité d’un ou plusieurs coefficients, par ex. Ho : β k = 0 versus Ha : β k ≠ 0 , on utilisera soit le test de Wald W, soit le test du rapport de vraisemblance LR. Dans le cas où l’on veut tester la significativité d’un seul coefficient ces statistiques s’écrivent :
W=
βˆ k SEˆ ( βˆ k )
→ N (0,1)
13
LR = −2 log (
Lc ) → χ 2 (1) Lc
tandis que si l’on veut tester la significativité de plusieurs coefficients, par ex. Ho : β 1 = β 2 = L = β M = 0 , alors elles s’écrivent :
(
W = βˆ ′ Vˆ ( βˆ )
)
LR = −2 log (
−1
βˆ
→ χ 2 (M )
Lc ) → χ 2 (M ) Lc
où Lc est la vraisemblance évaluée sous la contrainte Ho et Lc la vraisemblance non contrainte. NB : La statistique de Wald fait intervenir les expressions matricielles suivantes : 1 x11 L x1 p Pˆ1 (1 − Pˆ1 ) L 0 M M −1 ˆ ˆ ˆ ˆ V ( β ) = X ′V X , V = M O M et X = ˆ (1 − Pˆ ) 0 L P n n 1 x n1 L x np
(
)
Exercice 3 Dans cet exercice nous allons introduire un nouveau jeu de données qui nous servira jusqu’à la fin de ce cours afin d’illustrer les propos. Il s’agit des données « Low birth weight » issues d’une étude des facteurs de risque liés à la mise au monde d’un bébé de petit poids de naissance, i.e.< 2500g. L’échantillon concerne 189 femmes dont 59 ont eu un bébé pesant < 2500g. Les facteurs de risque potentiels évalués sont l’âge de la mère age (en années), son poids lors de ses dernières règles lwt (en livres), la race race (blanc, noir, autre), la fumée durant la grossesse smoke (oui/non), le nombre d’épisodes de contractions importantes avant terme ptl (0,1,2, etc), un antécédent de problème d’hypertension ht (oui/non), la présence d’une irritation utérine ui (oui/non) et le nombre de visites au médecin durant les trois premiers mois de grossesse ftv (0,1,2,etc.). Remarquons qu’on pourrait aussi étudier la relation entre le poids de naissance bwt (en grammes) et ces facteurs au moyen, cette fois, d’une régression linéaire puisque bwt est une variable continue. Eventuellement, il faudra au préalable transformer cette variable pour rendre sa distribution plus symétrique et sa variance plus stable.
14
1) Description des données Nous allons commencer par décrire nos données : fréquences, données manquantes, etc. Describe tab low summarize age summarize lwt tab race, missing tab smoke, missing tab ptl, missing tab ht, missing tab ui, missing tab ftv, missing
2) Analyse bivariable Avant d’analyser nos données au moyen d’un modèle de régression logistique multivariables il est d’usage de procéder à des analyses bivariables, en particulier lorsque le nombre de variables candidates à introduire dans le modèle est élevé. Ces analyses bivariables nous permettront d’appréhender les facteurs de risque potentiellement associés avec l’outcome. Sur la base de ces résultats, on procédera à un tri préalable de ces facteurs selon leur degré d’évidence (p-value) et nos connaissances théoriques, afin de ne pas tous les introduire dans le modèle (risque de multicolinéarité, difficulté d’interprétation des résultats, overfitting, etc.). Lorsque la variable explicative est continue on peut former des catégories afin de représenter graphiquement sa relation avec la variable dépendante. Lorsqu’une variable explicative catégorielle comporte des catégories n’ayant pas assez d’observations (e.g. <5) on procède à leur regroupement, afin d’obtenir des fréquences suffisamment élevées. 2.1) lorsque la variable explicative est continue * génération des percentiles de la variable age cap drop pct_age xtile pct_age=age, nquantiles(9) sort pct_age tab pct_age * calcul de la proportion de petits poids dans les catégories d'âge by pct_age: egen p_low=mean(low) * graphe de la relation entre les catégories d'age et la proportion de petits poids scatter p_low pct_age, ylabel(0(0.2)1) ytitle(P(low)) title(Proportion de petits poids /// en fonction de l'âge) saving(g6, replace) * une autre représentation scatter low age, jitter(2) ylabel(0(1)1) title(Petit poids de naissance en fonction de l'âge) /// || lowess low age, sort bwidth(1) * génération des percentiles de la variable lwt
15
cap drop pct_lwt xtile pct_lwt=lwt, nquantiles(9) sort pct_lwt * calcul de la proportion de petits poids dans les catégories de lwt cap drop p_low by pct_lwt: egen p_low=mean(low) tab pct_lwt * graphe de la relation entre les catégories d'lwt et la proportion de petits poids scatter p_low pct_lwt, ylabel(0(0.2)1) ytitle(P(low)) title(Proportion de petits poids /// en fonction du poids de la mère) saving(g7, replace) * une autre représentation scatter low lwt, jitter(2) ylabel(0(1)1) title(Petit poids de naissance en fonction du poids de la mère) /// || lowess low lwt, sort bwidth(1)
2.2) lorsque la variable explicative est catégorielle tab low race, chi2 row col tab low smoke, chi2 row col tab low ptl, chi2 row col * Lorsqu’il y a des catégories qui sont très peu représentées on procède à un regroupement recode ptl (0=0) (1 2 3=1), gen(ptl_g) tab low ptl_g, chi2 row col tab low ht, chi2 row col tab low ui, chi2 row col tab low ftv, chi2 row col recode ftv (0=0) (1=1) (2=2) (*=3), gen(ftv_g) tab low ftv_g, chi2 row col
Sur la base de ces résultats l’on pourra pré-selectionner les variables candidates pour l’analyse multivariables. Les variables ayant une p-value supérieure à 0.2 auront peu de chance d’être retenues dans le modèle multivariables. S’il y a beaucoup de co-variables et que l’on ne peut pas toutes les introduire à la fois dans le modèle, on donnera une préférence à celles dont la pvalue est la plus petite. Il faudra, néanmoins, ultérieurement ré-introduire une à une ces variables dans le modèle multivariables pour ré-évaluer leur association.
16
4) La transformation logit Une transformation centrale dans l’analyse de régression logistique est la transformation « logit ». En effet, cette transformation permet d’établir une relation entre la probabilité de l’outcome et le prédicteur linéaire β 0 + β 1 x : P ( y = 1 | x) logit [P( y = 1 | x)] = log = β 0 + β1 x 1 − P( y = 1 | x) Elle s’interprète comme le logarithme du rapport des cotes p/(1-p). La transformation « logit » ou plus simplement le « logit » permet d’interpréter les résultats d’une estimation sur l’échelle « logit ». L’intérêt de raisonner sur l’échelle « logit » réside avant tout dans la possibilité d’évaluer approximativement d’un coup d’oeil la probabilité associée à une combinaison des co-variables, ainsi que l’importance relative de celles-ci. Voyons cela : la probabilité de y s’exprime à partir du « logit » comme suit : e log it [P ( y =1| x ) ] P ( y = 1 x) = 1 + e log it [P ( y =1| x ) ] Le « logit » peut prendre des valeurs entre -inf et +inf, mais la zone d’intérêt se situe entre -5 et +5, car au delà de ces limites la probabilité est soit 0 soit 1 :
0
.1
.2
.3
.4
P(y=1) .5 .6
.7
.8
.9
1
Probabilité de y en fonction du logit
-5
-4
-3
-2
-1
0 logit
1
2
3
4
5
figure 13
Par exemple, lorsque le logit vaut 0 la probabilité de y est de 0.5, tandis que lorsqu’il vaut +5 elle est de 0.993 et lorsqu’il vaut -5 de 0.007. Ainsi, à partir des résultats de l’estimation des coefficients il est facile de calculer le « logit » et d’évaluer approximativement la probabilité de l’outcome. Considérons l’exemple fictif suivant du résultat de l’estimation d’une régression logistique comportant les variables âge et sexe :
17
logit [P( y = 1 | âge, sexe)] = −5 + 0.1 × âge + 2 × sexe où la variable sexe prend la valeur 0 pour les femmes et 1 pour les hommes. Pour une femme d’âge 50 ans le « logit » est égal à 0 et, en se référant à la figure 13, on évalue la probabilité de l’outcome à 0.5 . Pour un homme, par contre, le « logit » prendrait la valeur 2 et la probabilité serait pratiquement de 0.9 . On constate, d’autre part, que l’effet du veillissement d’une année est vingt fois moins important (en terme d’augmentation du risque) que le changement de catégorie pour le genre. On remarque que, d’une part, plus le coefficient d’une co-variable est grand plus l’effet d’une variation unitaire de cette variable est important sur la probabilité de y, d’autre part, que lorsqu’on se situe sur l’échelle « logit » proche de 0 cette variation aura un effet plus marqué que lorsqu’on est proche de 3 ou -3. On peut se poser la question : à partir de quelle amplitude d’un coefficient un changement unitaire de la co-variable a un effet sensible sur la probabilité. D’après la figure 13 on sait que cet effet sera maximum lorsque le « logit » est proche de 0. Ainsi, un coefficient d’amplitude 0.2 engendrera au plus un changement de 5% de la probabilité, tandis que si il vaut 0.5 alors le changement est au plus de 12%. Un autre intérêt du « logit », comme on le verra dans le prochain chapitre, est sa relation avec une mesure d’association très utilisée en épidémiologie entre un facteur explicatif et l’outcome : l’Odds Ratio.
Exercice 4 Dans cet exercice, nous allons apprendre à raisonner sur l’échelle « logit ». Autrement dit, à évaluer directement le niveau de la probabilité associée à une combinaison des co-variables, ainsi que l’impact d’un accroissement unitaire d’une de ces co-variables. Pour cela, il est utile de bien avoir en tête la figure 13, en particulier les niveaux de probabilités associées à différent points entre -5 et +5 sur l’échelle « logit ». Afin de donner un sens à la constante du modèle, nous allons voir qu’il est utile de centrer les co-variables continues. Pour estimer les coefficients du modèle nous utiliserons la commande « logit ». La commande « logistic » permet aussi d’estimer le modèle mais fournit les résultats sous forme d’Odds Ratios). Remarque : les variables smoke, ht et ui ont été codées 0/1 de sorte qu’on peut directement les utiliser dans le modèle sans les préfixer par « i. ». La variable race, en revanche, a été codée 1, 2 et 3 et il faut créer des variables binaires 0/1 pour représenter les différentes catégories. Pour cela, STATA possède une commande automatique « xi : » qui créera les variables binaires nécessaires pour toutes les co-variables catégorielles préfixées par « i. ». Cette fois, nous allons considérer un modèle multivariables afin d’étudier l’effet conjoint de plusieurs co-variables sur la probabilité de petit poids de naissance.
18
1) Analyse multivariables Dans cet exercice nous allons illustrer l’effet du centrage des co-variables continue. * estimation sans centrage des covariables continues xi: logit low age lwt i.race smoke ptl_g ht ui i.ftv_g * estimation avec centrage des covariables continues egen mean_age=mean(age) gen age_c=age-mean_age egen mean_lwt=mean(lwt) gen lwt_c=lwt-mean_lwt xi: logit low age_c lwt_c i.race smoke ptl_g ht ui i.ftv_g
En comparant ces deux estimations l’on peut constater que seule la constante change. Dans le premier cas, i.e. sans centrage, la constante n’a pas de sens puisqu’elle correspond à une femme d’âge 0 et de poids 0 kg aux dernières règles… Dans, le 2ième modèle, en revanche, la constante a l’honorable rôle de représenter une femme d’âge moyen et de poids moyen. 2) Déterminez une combinaison de co-variables de sorte que la probabilité prédite de poids de naissance <2500 kg soit d’au moins 0.5 Pour cela, il vous faut calculer la valeur du « logit » pour différents niveaux des co-variables. A vous de proposer des valeurs… 3) Effet du changement d’unités de mesure Le poids des femmes aux dernières règles est mesuré en [livres]. Afin d’interpréter les résultats en [kilogrammes] nous allons recoder lwt et ré-estimer le modèle. * recodage du poids aux dernières règles en kg gen lwt_kg_c=lwt_c/2 xi: logit low age_c lwt_kg_c i.race smoke ptl_g ht ui i.ftv_g * recodage du poids aux dernières règles en 10 kg gen lwt_10kg_c=lwt_c/20 xi: logit low age_c lwt_10kg_c i.race smoke ptl_g ht ui i.ftv_g
On constate que sur l’échelle « logit » le changement d’unités affecte le coefficient estimé de manière proportionnelle. Remarquons que la probabilité estimée, par contre, change de manière non proportionnelle puisque la relation entre le « logit » et la probabilité de l’outcome est non linéaire.
19
4) Test de « significativité » des coefficients (test de Wald) Pour les variables explicatives dichotomiques ou continues STATA nous fournit directement la p-value du test de « significativité » du coefficient, tandis que pour les variables explicatives polytomiques il faut invoquer le test de Wald au moyen de la commande « test ». Nous allons tester si les variables race et ftv_g sont significatives. Rappelons que dans l’exercice précédent nous avions regroupé les catégories de la variable ftv en créant la nouvelle variable ftv_g afin d’augmenter les effectifs dans les cellules : * test de significativité de la variable ftv_g test _Iftv_g_1 _Iftv_g_2 _Iftv_g_3 * test de significativité de la variable race test _Irace_2 _Irace_3
5) Sensibilité de la probabilité à un changement unitaire d’une co-variable Nous allons évaluer l’impact sur la probabilité d’un changement unitaire d’une co-variable en fonction de la position sur l’échelle « logit ». Pour illustrer ceci, nous allons calculer la probabilité pour différents accroissements et valeurs du « logit » : * lorsque le logit est proche de 0 disp "prob="exp(0)/(1+exp(0)) disp "prob="exp(0+0.2)/(1+exp(0+0.2)) disp "prob="exp(0+0.5)/(1+exp(0+0.5)) * lorsque le logit est proche de -2 disp "prob="exp(-2)/(1+exp(-2)) disp "prob="exp(-2+0.2)/(1+exp(-2+0.2)) disp "prob="exp(-2+0.5)/(1+exp(-2+0.5))
En conclure que pour avoir un effet unitaire suffisamment sensible il faut que le coefficient d’une co-variable soit au moins d’amplitude 0.5 . Remarquer que le choix des unités de mesure est primordial pour cette interprétation.
Dans le prochain chapitre nous introduirons la notion d’Odds Ratio qui est intimement liée à la transformation « logit ».
20
5) Le succès du modèle Logit : l’Odds Ratio Si le modèle Logit est très utilisé en épidémiologie c’est avant tout à cause de l’interprétation de l’exponentielle du coefficient d’une co-variable comme un Odds Ratio. Pour comprendre ce que représente un Odds Ratio voyons comment il est défini. Pour cela, considérons un modèle avec une seule variable explicative dichotomique comme le sexe (le cas plus général du modèle de régression multiple incorporant plusieurs co-variables ainsi que des interactions sera abordé dans la section suivante) et adoptons le codage suivant : « 0 » pour les femmes et « 1 » pour les hommes, de sorte qu’on écrira la probabilité P(y = 1) = p0 pour les femmes et P(y = 1) = p1 pour les hommes.
5.1) L’Odds Ratio comme mesure d’association Un Odds est défini comme le rapport des cotes :
Odds =
p 1− p
où p est par exemple la probabilité de gagner. On définit l’Odds Ratio (OR) associé à la variable sexe comme suit : p1 1 − p1 OR = p0 1 − p0 Si p0 représente la probabilité d’être malade pour une femme et p1 celle pour un homme, alors un Odds Ratio de 1 signifie que la probabilité d’être malade est la même chez les hommes et chez les femmes. Autrement dit, le risque de maladie n’est pas associé au sexe. En revanche, un Odds Ratio différent de 1 signifie qu’il y a une association entre la maladie et le genre. Si cet Odds Ratio est >1 cela signifie que le numérateur est plus grand que le dénominateur et, par conséquent, que les hommes ont un plus grand risque d’être malade que les femmes. C’est le contraire s’il est <1. Revenons à notre modèle Logit comportant comme variable explicative uniquement le sexe : logit [P( y = 1 | sexe)] = β 0 + β 1 sexe
Pour les hommes on a : 21
logit [P( y = 1 sexe = 1)] = β 0 + β 1
et pour les femmes : logit [P( y = 1 sexe = 0)] = β 0 En utilisant la relation entre la probabilité de y et le « logit » vue dans la section précédente on obtient : p1 1 − p1 e log it [P ( y =1 sexe=1) ] e β 0 + β1 OR = = log it [P ( y =1 sexe=0 ) ] = β 0 = e β1 p0 e e 1 − p0
de sorte que dans un modèle Logistique l’exponentielle du coefficient d’une variable explicative s’interprète comme son Odds Ratio.
5.2) L’Odds Ratio comme mesure du risque relatif (RR) De façon analogue à la définition de l’Odds Ratio dans la section précédente, on définit le Risque Relatif (RR) associé à la variable sexe comme :
RR =
p1 p0
Cette grandeur a une interprétation intuitive claire, ce qui n’est pas le cas de l’Odds Ratio. Lorsque la prévalence de l’événement à expliquer est faible, i.e. p0 et p1 sont petites, l’Odds Ratio fournit une approximation du risque relatif :
OR =
p1 (1 − p 0 ) ≅ 1 p1 × ≅ = RR p 0 (1 − p1 ) ≅ 1 p 0
Cependant, lorsque ces prévalences ne sont pas tout petites on a, en général , OR ≠ RR :
22
OR =
(1 − p 0 ) p1 (1 − p 0 ) × = RR × p 0 (1 − p1 ) (1 − p1 )
Afin d’illustrer ce dernier point, considérons un exemple issu d’une étude transversale portant sur 170 enfants âgés de 24 à 36 mois d’une région rurale africaine où l’on s’intéressait à l’association entre le retard de croissance staturale et le petit poids de naissance (< 2500g). Les données sont récapitulées dans le tableau suivant : | retard de croissance staturale row |
1
0 |
Total
-----------+----------------------+---------<2500g
|
18
13 |
31
|
58.06
41.94 |
100.00
|
37.50
10.66 |
18.24
-----------+----------------------+---------=2500g
|
30
109 |
139
|
21.58
78.42 |
100.00
|
62.50
89.34 |
81.76
-----------+----------------------+---------Total |
48
122 |
170
|
28.24
71.76 |
100.00
|
100.00
100.00 |
100.00
Dans cet exemple, le risque de retard de croissance staturale chez les petits poids de naissance est 18/31=0.58, tandis qu’il est de 30/139=0.22 chez les autres, des sorte que :
RR =
18 P(retard croissance | < 2500 g ) 31 = 2.69 = 30 P(retard croissance | ≥ 2500 g ) 139 P(retard croissance | < 2500g )
18 13 = 5.03 OR = = P(retard croissance | ≥ 2500g ) 30 109 1 − P(retard croissance | ≥ 2500 g ) 1 − P(retard croissance | < 2500g )
On constate, ici, que l’Odds Ratio sur-estime le risque relatif de façon importante, ce qui n’est pas surprenant puisque les prévalences p0=0.22 et p1=0.58 ne sont pas petites.
23
Exercice 5 Dans cet exercice, nous allons reconsidérer l’estimation du modèle de l’exercice 4 avec les données « Low birth weight » et interpréter les coefficients estimés en terme d’Odds Ratios. L’Odds Ratio a une interprétation claire et intuitive uniquement lorsqu’il fournit une « bonne » approximation du risque relatif, ce qui est le cas lorsque la prévalence de l’outcome est petite dans les deux catégories considérées. Autrement, il fournit une mesure d’association qui, ma fois, n’est pas interprétable clairement : que signifie concrètement le ratio de deux autres ratios ?!. Pour obtenir les résultats de l’estimation du modèle Logistique en terme d’Odds Ratios on utilisera la commande « logistic » de STATA. Remarquons que, cette fois, cela n’a pas d’importance si les variables continues ont été centrées ou pas, car le calcul de l’Odds Ratio ne fait pas intervenir la constante du modèle. Aussi, dans l’exercice 4 on testait la « significativité » d’une covariable en testant si son coefficient était significativement différent de 0. Lorsqu’on travaille avec les Odds Ratios, le test porte alors sur la valeur 1. Autrement dit, une variable sera « significativement » associée à l’outcome si son Odds Ratio est « significativement » différent de 1. 1) Comparaison des p-values et OR avec et sans centrage xi: logistic low age_c lwt_c i.race smoke ptl_g ht ui i.ftv_g xi: logistic low age lwt i.race smoke ptl_g ht ui i.ftv_g
On vérifie que le centrage des co-variables continue n’affecte pas les Odds Ratios. xi: logit low age lwt i.race smoke ptl_g ht ui i.ftv_g xi: logistic low age lwt i.race smoke ptl_g ht ui i.ftv_g
On vérifie que les p-values des paramètres estimés avec la commande « logit » sont bien les mêmes que celles estimées avec la commande « logistic » . 2) Interprétation des OR Pour commencer, nous allons estimer un modèle avec uniquement la variable explicative smoke : logistic low smoke
L’ OR s’interprète comme une mesure d’association. S’il est supérieur à 1 la relation est croissante, et décroissante s’il est inférieur à 1. Lorsqu’il est égal à 1 il n’y a pas d’association. Afin de d’anticiper si l’on peut espérer que l’OR associé à la variable smoke fournit, ici, une « bonne » approximation du RR nous allons calculer la prévalence du risque de petit poids de naissance dans les deux catégories de la variable d’exposition. * prévalences de l'outcome dans les 2 catégories de smoke tab low smoke, row col * alternativement logistic
en
utilisant
les
résultats
fournis
par
la
commande
24
disp "p0=" exp(_b[_cons])/(1+exp(_b[_cons])) disp "p1=" exp(_b[_cons]+_b[smoke])/(1+exp(_b[_cons]+_b[smoke]))
En conclure que les prévalences sont élevées et que l’approximation du RR par l’OR est susceptible d’être très imprécise. * calcul du RR associé à la variable smoke et comparaison avec son OR disp "RR=" exp(_b[smoke])*(1+exp(_b[_cons]))/(1+exp(_b[_cons]+_b[smoke])) disp "OR=" exp(_b[smoke])
En conclure que l’approximation du RR associé à la variable smoke par son OR conduit, ici, à une surestimation. Les OR et RR que nous venons d’estimer sont non ajustés puisque le modèle comporte uniquement la variable d’exposition. Comparez l’OR non ajusté avec l’OR ajusté et en conclure, qu’ici, ils diffèrent peu. * OR non ajusté/ajusté logistic low smoke xi: logistic low age lwt i.race smoke ptl_g ht ui i.ftv_g
Le calcul du RR ajusté est plus complexe et sera abordé dans le prochain chapitre.
25
26
6) L’interprétation des coefficients Nous avons vu dans la section précédente que dans le cas d’un modèle comportant une seule variable explicative dichotomique l’exponentielle du coefficient de cette variable s’interprétait comme un Odds Ratio. Voyons ce qui se passe lorsque la variable explicative admet plusieurs catégories, i.e. elle est polytomique, ou qu’elle est continue, ou encore que le modèle incorpore d’autres co-variables ainsi que des interactions.
6.1) Le cas d’un modèle additif, i.e. sans interactions Un modèle est additif2 lorsque les co-variables x1, x2, …, xp entrent dans le modèle de manière additive sans faire intervenir le produit d’une variable avec une autre :
[
]
logit P( y = 1 | x1 , L , x p ) = β 0 + β 1 x1 + β 2 x 2 + ... + β p x p où β 0 est la constante du modèle. Pour illustrer, considérons le modèle suivant : logit [P( y = 1| âge, sexe)] = β 0 + β 1 âge + β 2 sexe où les variables explicatives sont l’âge et le sexe. Il s’agit d’un modèle additif car il n’y a pas d’interaction (de produit) entre les variables âge et sexe. Autrement dit, dans ce modèle on postule que l’effet de l’âge et du sexe sont indépendants (sur l’échelle logit). Graphiquement, cette hypothèse implique que la droite représentant l’effet de l’âge est simplement translatée sur une distance β 2 lorsqu’on passe d’un genre à l’autre.
Relation entre le logit et l’âge chez les femmes et les hommes dans un modèle additif 2
femmes
logit
hommes 0
β2
-2 âge
figure 14
2
Dans le cas de la régression logistique, le modèle est additif sur l’échelle « logit », mais multiplicatif lorsqu’on considère la probabilité.
27
Dans cet exemple, le vieillissement a le même effet chez les hommes et chez les femmes, mais le niveau absolu du risque est différent (les deux droites ne sont pas superposées). Autrement dit, un accroissement unitaire de l’âge augmentera le logit du même montant quel que soit le genre, et l’ Odds Ratio associé à la variable âge sera le même pour les hommes et les femmes. Remarques (*) 1) Même si dans ce cas l’Odds Ratio associé à la variable âge est le même pour les femmes et les hommes, le risque relatif est différent si β 2 ≠ 0 . En effet, cela est dû au fait que le niveau absolu du risque est plus bas, dans cet exemple, chez les hommes et l’effet de l’accroissement d’une année d’âge n’augmente pas la probabilité P(y = 1) du même montant (même si le logit change de la même quantité). La raison provient de la relation non linéaire entre le logit et la probabilité P(y). Par exemple, si le logit passe de 1 à 2, la probabilité passe de 0.73 à 0.88, tandis que si le logit passe de 3 à 4, alors la probabilité passe de 0.95 à 0.98. 2) Dans les modèles non-linéaires, comme le modèle Logit, même si l’on introduit pas de terme produit croisé de deux co-variables celles-ci présentent, en général, une interaction (1). Soit le modèle de régression non linéaire : E ( y = 1| âge, sexe) = f (β 0 + β 1 âge + β 2 sexe + β 12 âge × sexe ) L’effet de l’interaction entre âge et sexe se calcule comme : ∂ 2 f (.) = f ′′(.) × ( β 2 + β 12 âge) × ( β 1 + β 12 sexe) + f ′(.) × β 12 ∂ âge ∂ sexe de sorte que même si β 12 = 0 cette expression ne s’annule pas. Ce phénomène est propre aux modèles de régression non-linéaire. 3) La définition que nous avons adoptée d’une interaction est justifiée si l’on travaille sur l’échelle logit et que l’on s’intéresse à l’effet conjoint de deux co-variables sur l’Odds Ratio et non pas sur le Risque Relatif. a) La constante du modèle La constante du modèle s’interprète comme « l’effet » de la catégorie de référence. Autrement dit, β 0 permet de calculer la probabilité de y lorsque toutes les co-variables x1, x2, …, xp sont nulles. Si l’on revient à notre exemple d’un modèle contenant l’âge et le sexe comme variables explicatives : logit [P( y = 1| âge, sexe)] = β 0 + β 1 âge + β 2 sexe
28
Nous avons arbitrairement choisi de coder les valeurs de la variable sexe = 0 pour les femmes et sexe = 1 pour les hommes, de sorte que β 0 s’interprète comme le logit de la probabilité d’une femme d’âge 0. En effet, la probabilité P(y = 1| âge et sexe), e.g. d’être malade en fonction de son âge et sexe, s’écrit :
P( y = 1 | âge, sexe) =
e β 0 + β1 âge + β 2 sexe 1 + e β 0 + β1 âge + β 2 sexe
de sorte que pour une femme d’âge 0 on obtient :
P( y = 1 | âge = 0, sexe = 0) =
e β 0 + β1 × 0 + β 2 × 0 e β0 = 1 + e β 0 + β1 × 0 + β 2 × 0 1 + e β 0
sa probabilité ne dépend que de β 0 . Pour un homme d’âge 0, en revanche, la probabilité dépend aussi de β 1 : e β 0 + β1 ×0+ β 2 ×1 e β 0 + β1 P( y = 1 | âge = 0, sexe = 1) = = 1 + e β 0 + β1 ×0+ β 2 ×1 1 + e β 0 + β1 Remarque Pour que la constante du modèle admette une interprétation plus honorable que le logit pour une femme d’âge 0, il est préférable de centrer la variable âge, âge_c = âge-moyenne(âge). Dans ce cas, la constante s’interprète comme le logit d’une femme d’âge égal à l’âge moyen dans l’échantillon. b) Coefficient d’une variable explicative dichotomique Lorsque la variable explicative est dichotomique l’exponentielle du coefficient de cette variable s’interprète comme l’Odds Ratio (OR) associé au passage de la catégorie de référence 0 à la catégorie 1. Ainsi, dans notre exemple, lorsque la variable sexe passe de 0 à 1, on a : log it [P ( y =1 âge , sexe =1) ]
e e β 0 + β 1 âge + β 2 OR = log it [P ( y =1 âge , sexe = 0 ) ] = β 0 + β 1 âge = e β 2 e e
29
Il s’agit d’un Odds Ratio ajusté puisque modèle comporte en plus de la variable d’exposition sexe la variable explicative âge. Remarquons que l’Odds Ratio ajusté est en général différent de celui non ajusté, même si son calcul ne fait pas intervenir directement la variable âge, car l’estimation de β 2 dépend de celle de β 1 . Remarque Le calcul du Risque Relatif (RR) est plus complexe et fait intervenir toutes les co-variables du modèle : p1 (1 − p1 ) 1 − e β 0 + β1 âge + β 2 (1 + e β 0 + β1 âge + β 2 ) β2 RR = = OR × =e × p0 (1 − p 0 ) 1 − e β 0 + β1 âge (1 + e β 0 + β1 âge ) = e β2 ×
1 + e β 0 + β1 âge 1 + e β 0 + β1 âge + β 2
Il s’agit d’un RR ajusté puisque modèle comporte en plus de la variable d’exposition sexe la variable explicative âge. Remarquons que le RR ajusté est non seulement différent de l’OR ajusté, mais qu’en plus il n’est pas constant, il dépend des valeurs de la co-variable âge. c) Coefficient d’une variable explicative polytomique Lorsque la variable explicative est polytomique, i.e. elle admet plus de deux catégories, on choisi l’une des catégories comme référence et l’on calcule des Odds Ratios pour les autres catégories par rapport à cette référence. Considérons par exemple la variable éducation comportant 3 niveaux : 1 pour niveau « fin de scolarité », 2 pour « apprentissage » et 3 pour « études supérieures ». Pour représenter une telle variable l’on considérera un modèle avec, en plus de la constante, deux variables « indicatrice » ou « dummy » prenant la valeur 1 si l‘individu possède l’attribut et 0 sinon : 1 si apprentissage D1 = 0 sinon 1 si études supérieures D2 = 0 sinon et le logit s’écrit : logit [P( y = 1| éducation)] = β 0 + β 1 D1 + β 2 D2
L’Odds Ratio associé au passage de la catégorie 1 « fin de scolarité » à la catégorie 2 « apprentissage » est : 30
OR =
log it [P ( y =1 éducation = 2 ) ]
e e
log it [P ( y =1
e β 0 + β1 = β 0 = e β1 éducation =1) ] e
Tandis que celui associé au passage de la catégorie 1 « fin de scolarité » à la catégorie 3 « études supérieures » est :
OR =
log it [P ( y =1 éducation =3 ) ]
e e
log it [P ( y =1
e β0 + β2 = = e β2 éducation =1) ] e β0
Ces 2 Odds Ratio sont directement fournis par le programme. Si, en revanche, l’on désire l’Odds Ratio associé au passage de la catégorie 2 « apprentissage » à la catégorie 3 « études supérieures » il faut calculer :
OR =
log it [P ( y =1 éducation = 3) ]
e β0 + β2 = = e β 2 − β1 β 0 + β1 log it [P ( y =1 éducation = 2 ) ] e e e
d) Coefficient d’une variable explicative continue Lorsque la variable explicative est continue on calcule un Odds Ratio associé à un accroissement unitaire. Par exemple, considérons la variable âge mesurée en années et supposons que la personne soit d’âge x. Le vieillissement d’une année est associé à un Odds Ratio donné par l’expression :
OR =
e
log it [P ( y =1 âge = x +1) ]
e
log it [P ( y =1
e β 0 + β1 ( x +1) = = e β1 β 0 + β1 ( x ) âge = x ) ] e
Remarques 1) On notera que cet Odds Ratio dépend des unités et si l’on avait mesuré l’âge en décades on aurait obtenu un Odds Ratio plus élevé :
OR =
e
log it [P ( y =1 âge = x +10 ) ]
e
log it [P ( y =1
e β 0 + β1 ( x +10) = β 0 + β1 ( x ) = e β1×10 = e β1 âge = x ) ] e
( )
10
Ainsi, pour calculer un Odds Ratio associé à un accroissement de « z » années il suffit d’élever à la puissance « z » l’Odds Ratio calculé pour une année.
31
2) L’Odds Ratio ne fait pas intervenir la valeur x prise par la variable âge, tandis que c’est le cas du Risque Relatif : p1 (1 − p1 ) 1 − e β 0 + β1 ( x +1) (1 + e β 0 + β1 ( x +1) ) 1 + e β 0 + β1 ( x ) β1 = OR × = e β1 × = e × p0 (1 − p0 ) 1 − e β 0 + β1 ( x ) (1 + e β 0 + β1 ( x ) ) 1 + e β 0 + β1 ( x +1)
RR =
Le Risque Relatif associé à un accroissement unitaire de la variable âge dépend du niveau de référence x de l’âge, tandis que ce n’est pas le cas de l’Odds Ratio. Le Risque relatif et l’Odds ratio sont deux mesures différentes d’associations entre l’exposition et l’outcome. Parfois, ces deux grandeurs sont proches mais, en général, elles n’ont pas la même interprétation. 3) Comme le RR dépend des valeurs des co-variables dans une étude donnée l’on a autant de valeurs que de participants. Pour simplifier la présentation des résultats on peut calculer le RR moyen. La formule suivante fournit une approximation (2) du RR moyen et permet d’effectuer le calcul directement à partir de l’OR estimé et la prévalence de l’outcome chez les non exposés p0 :
RRmoyen ≅
OR (1 − p 0 ) + ( p 0 × OR)
e) L’Odds ratio associé à la variation de plusieurs co-variables Revenons à notre modèle comportant les variables âge et sexe et calculons l’Odds Ratio associé à une variation simultanée des deux co-variables : logit [P( y = 1| âge, sexe)] = β 0 + β 1 âge + β 2 sexe où, comme avant, sexe prend la valeur 0 pour les femmes et 1 pour les hommes. L’Odds Ratio associé au passage de la catégorie 0 à 1 pour le genre et à l’augmentation de l’âge de ∆ unités s’écrit :
OR =
e
log it [P ( y =1 âge = x + ∆ , sexe =1) ]
e
log it [P ( y =1
e β 0 + β1 ( x + ∆ ) + β 2 = = e β1 × ∆ + β 2 β 0 + β1 ( x ) âge = x , sexe = 0 ) ] e
6.2) Le cas d’un modèle non additif, i.e. avec interactions Un modèle n’est pas additif (sur l’échelle logit) lorsqu’il contient non seulement les covariables x1, x2, …, xp mais en plus le produit de deux ou plusieurs co-variables :
[
]
logit P( y = 1 | x1 , L , x p ) = β 0 + β 1 x1 + β 2 x 2 + ... + β p x p + β p +1 x1 x 2 + β p + 2 x1 x3 + ...
32
Pour illustrer, considérons le modèle suivant : logit [P( y = 1| âge, sexe)] = β 0 + β 1 âge + β 2 sexe + β 12 âge × sexe où les variables explicatives sont l’âge et le sexe, et où l’on postule qu’il peut exister une interaction entre les variables âge et sexe. Concrètement, cela signifie que sur l’échelle logit ou en terme d’Odds ratio l’effet de l’âge n’est pas indépendant de celui du sexe (i.e. l’Odds Ratio des femmes n’est pas le même que celui des hommes). Autrement dit, pour étudier l’effet de l’âge il faut aussi spécifier pour quelle niveau de la variable sexe, puisque cet effet est différent pour les deux sexes. Graphiquement, l’existence d’une interaction implique que la droite représentant l’effet de l’âge pour les femmes n’a pas la même pente que pour les hommes.
Relation entre le logit et l âge chez les femmes et les hommes dans un modèle non additif (avec interaction) hommes
logit
2
femmes 0
-2 âge
figure 15
Le calcul d’un Odds Ratio dans un modèle non additif se traite de la même manière que lorsqu’on considère le calcul d’un Odds Ratio pour la variation simultanée de plusieurs covariables. Pour illustrer, considérons toujours notre exemple d’un modèle contenant les co-variables âge et sexe, ainsi que l’interaction âge × sexe . Alors, l’Odds Ratio associé à l’incrément d’une unité de la variable âge s’écrit :
OR =
e
log it [P ( y =1 âge +1, sexe, ( âge +1) × sexe ) ]
e
log it [P ( y =1 âge , sexe , âge × sexe ) ]
e β 0 + β1 ( âge +1)+ β 2 sexe+ β12 ( âge +1) ×sexe = β 0 + β1 ( âge )+ β 2 sexe+ β12 ( âge ) × sexe = e β1 + β12 sexe e
On constate que s’il s’agit d’une femme on a sexe = 0 et OR = e β1 + β12 ×0 = e β1 , tandis que s’il s’agit d’un homme on a sexe = 1 et OR = e β1 + β12 ×1 = e β1 + β12 . En ce qui concerne le Risque Relatif associé à l’accroissement unitaire de l’âge on a :
33
RR = OR ×
1 − p1 1 + e β 0 + β1âge + β 2 sexe+ β12 âge × sexe = e β1 + β12 sexe × 1 − p0 1 + e β 0 + β1 ( âge +1) + β 2 sexe+ β12 ( âge +1) ×sexe
Contrairement à l’OR, le RR dépend non seulement de la variable sexe, mais aussi de la variable âge. Remarque On distingue les interactions d’ordre 1 de celles d’ordre supérieures. Lorsque l’interaction fait intervenir uniquement le produit de deux co-variables elle est d’ordre 1, tandis que si elle définie comme le produit de trois co-variables il s’agit d’une interaction à l’ordre 2, etc. En général, pour des questions d’interprétation on se limite à des interactions d’ordre 1.
Exercice 6 Dans cet exercice nous allons calculer et comparer les Odds Ratios et les Risques Relatifs afin d'évaluer empiriquement la différence entre ces deux mesures d’association. Nous continuons notre présentation avec les données « Low birth weight ». 1) Modèle additif Pour commencer, nous allons considérer un modèle additif comportant les variables explicatives age et smoke. logistic low age smoke
1.1) Comparaison de l’OR et du RR associés à la variable smoke (modèle additif) * NB: le RR dépend non seulement de smoke mais aussi de l'âge cap drop logit* p1* p0* RR* OR * logit chez les fumeurs gen logit1=_b[_cons]+_b[age]*age+_b[smoke]*1 * logit chez les non fumeurs gen logit0=_b[_cons]+_b[age]*age+_b[smoke]*0 gen p0=exp(logit0)/(1+exp(logit0)) gen p1=exp(logit1)/(1+exp(logit1)) gen RR=p1/p0 summarize RR disp "0R=" exp(_b[smoke]) gen OR=exp(_b[smoke]) scatter RR age, ytitle(OR/RR) ylabel(1(0.2)2) title(Modèle additif: OR et RR associés à la cigarette, size(medsmall)) || scatter OR age, saving(g8, replace)
34
On constate que le RR associé à la cigarette dépend de l’âge, tandis que ce n’est pas le cas de l’OR. Conclure que dans notre cas l’OR sur-estime le RR. On peut obtenir une approximation du RR moyen à partir de la prévalence p0 : tab low smoke, row col disp "RR0=" exp(_b[smoke])/((1-0.2522)+(0.2522*exp(_b[smoke])))
1.2) OR et RR associés à la variable age (modèle additif) * NB: le RR dépend non seulement de l'âge mais aussi de smoke. On calculera le RR pour un fumeur et pour un non fumeur. cap drop logit* p1* p0* RR* OR * logit chez les non fumeurs gen logit1=_b[_cons]+_b[age]*(age+1)+_b[smoke]*0 gen logit0=_b[_cons]+_b[age]*age+_b[smoke]*0 gen p0=exp(logit0)/(1+exp(logit0)) gen p1=exp(logit1)/(1+exp(logit1)) gen RR_non_f=p1/p0 summarize RR_non_f * logit chez les fumeurs gen logit1f=_b[_cons]+_b[age]*(age+1)+_b[smoke]*1 gen logit0f=_b[_cons]+_b[age]*age+_b[smoke]*1 gen p0f=exp(logit0f)/(1+exp(logit0f)) gen p1f=exp(logit1f)/(1+exp(logit1f)) gen RR_f=p1f/p0f summarize RR_f gen OR=exp(_b[age]) scatter RR_non_f RR_f age, ytitle(OR/RR) title(Modèle additif: OR et RR associés à l'âge, size(medsmall)) || scatter OR age, saving(g9, replace)
Remarquer que le RR associé à l’accroissement unitaire de la variable continue age dépend non seulement des autres co-variables, mais aussi du niveau de l’âge. 1.3) OR associé aux variables age et smoke (modèle additif) Calculer l’OR associé à un accroissement unitaire de la variable age ainsi qu’au passage de la catégorie non fumeur à fumeur. disp "0R=" exp(_b[age]+_b[smoke]) * avec la commande lincom lincom age + smoke
1.4) OR associé à la variable race Nous allons calculer des OR dans le cas d’une variable explicative polytomique. gen black=cond(race==2,1,0) gen other=cond(race==3,1,0)
35
xi:logistic low age smoke i.race * OR associé au passage de la catégorie 2 à 3 disp "0R=" exp(_b[_Irace_3]-_b[_Irace_2]) lincom _Irace_3-_Irace_2
2) Modèle avec interaction Nous allons considérer, cette fois, un modèle non additif comportant une interaction entre les variables age et smoke. gen age_smoke=age*smoke gen age_c_smoke=age_c*smoke logit low age smoke age_smoke
2.1) OR et RR associés à la variable smoke (modèle avec interaction) * NB: dans un modèle avec interaction l'OR dépend de l'âge cap drop logit* p1* p0* RR* OR * logit chez les fumeurs gen logit1=_b[_cons]+_b[age]*age+_b[smoke]*1+_b[age_smoke]*age*1 * logit chez les non fumeurs gen logit0=_b[_cons]+_b[age]*age+_b[smoke]*0+_b[age_smoke]*age*0 gen p0=exp(logit0)/(1+exp(logit0)) gen p1=exp(logit1)/(1+exp(logit1)) gen RR=p1/p0 summarize RR gen OR=exp(_b[smoke]*1+_b[age_smoke]*age*1) summarize OR scatter OR RR age, ytitle(OR et RR) title(Modèle avec interaction: OR et RR associés à la cigarette, size(medsmall)) saving(g10, replace) * ex: pour un fumeur de 20, 50 ans disp "0R=" exp(_b[smoke]*1+_b[age_smoke]*20*1) lincom 20*age_smoke + smoke, or disp "0R=" exp(_b[smoke]*1+_b[age_smoke]*50*1) lincom 50*age_smoke + smoke, or
Constater que dans ce modèle avec interaction l’OR dépend de l’âge. 2.2) OR et RR associés à la variable age (modèle avec interaction) cap drop logit* p1* p0* RR* OR* * logit chez les non fumeurs gen logit1=_b[_cons]+_b[age]*(age+1)+_b[smoke]*0+_b[age_smoke]*(age+1)*0 gen logit0=_b[_cons]+_b[age]*age+_b[smoke]*0+_b[age_smoke]*age*0 gen OR_non_f=exp(_b[age]*1+_b[age_smoke]*1*0)
36
gen p0=exp(logit0)/(1+exp(logit0)) gen p1=exp(logit1)/(1+exp(logit1)) gen RR_non_f=p1/p0 summarize RR_non_f * logit chez les fumeurs gen logit1f=_b[_cons]+_b[age]*(age+1)+_b[smoke]*1+_b[age_smoke]*(age+1)*1 gen logit0f=_b[_cons]+_b[age]*age+_b[smoke]*1+_b[age_smoke]*age*1 gen OR_f=exp(_b[age]*1+_b[age_smoke]*1*1) gen p0f=exp(logit0f)/(1+exp(logit0f)) gen p1f=exp(logit1f)/(1+exp(logit1f)) gen RR_f=p1f/p0f summarize RR_f scatter RR_non_f RR_f age, ytitle(OR/RR) title(Modèle avec interaction: OR et RR associés à l'âge, size(medsmall)) || scatter OR_non_f OR_f age, saving(g11, replace) graph combine g8.gph g9.gph g10.gph g11.gph, iscale(.55)
Remarquer que le RR associé à l’âge est pratiquement constant chez les fumeurs. Ceci provient du fait que le coefficient associé à la variable age est très petit et que lorsqu’il s’agit d’un fumeur l’effet total est encore plus petit. La figure comportant les 4 graphes est instructive et permet de comparer les OR et RR estimés dans un modèle additif avec ceux estimés dans un modèle avec interaction. Le but de cet exercice était d’illustrer la différence entre OR et RR. On constate que la différence est parfois importante, mais que les deux mesures d’association vont dans le même sens lorsqu’il y a augmentation ou diminution du risque. On retiendra, donc, que l’OR qui est plus simple à calculer que le RR fournit une « bonne » approximation du RR pour autant que la prévalence de l’outcome soit faible (<10%). Si cette prévalence est élevée il peut y avoir une grande différence entre les deux mesures d’association et il faudra être prudent dans l’interprétation de l’OR. Il sera raisonnable de dire « il y a une association croissante… » si l’OR>1, mais il ne faudra pas considérer qu’il s’agit d’une « bonne » approximation du RR.
37
38
7) Stratégie de modélisation Sous le vocable de « stratégie de modélisation » on entend la problématique de choisir un modèle mathématique/statistique, choisir les co-variables, estimer les paramètres du modèle, analyser l’adéquation du modèle aux données et le valider. Il s’agit d’un processus complexe nécessitant de l’expérience.
Pourquoi construire un modèle ? Revenons, brièvement, sur le « pourquoi » construire un modèle !? La raison est qu’à partir d’un modèle décrivant nos données nous espérons pouvoir inférer des caractéristiques de la population d’intérêt, comme l’association entre l’âge et le risque de maladies coronariennes, montrer que le traitement A est plus efficace que le B, que le risque est différent selon le genre, etc. Afin de pouvoir extrapoler les résultats de l’analyse de l’échantillon de données dont on dispose à l’entier de la population il faut : 1. Disposer d’un échantillon représentatif de la population d’intérêt. C’est, en général, le cas lorsqu’il y a eu tirage aléatoire des observations. Parfois les observation ne proviennent pas exactement d’un tirage aléatoire mais elles ont été récoltées de façon plus ou moins aléatoire et l’échantillon est supposé représentatif (ex : étude de cohorte). Toutefois, dans ce cas, il y a possiblement un biais de sélection. 2. S’assurer que le modèle décrive « bien » les données et qu’il n’y a pas de biais systématique. Si ces deux conditions sont vérifiées ont considère qu’il est raisonnable d’inférer les résultats des analyses à l’ensemble de la population d’intérêt. Un modèle est par définition une représentation simplifiée de la réalité. Par conséquent, il est forcément, dans un certain sens, approximatif et incorrect. Néanmoins, sous les deux conditions que nous avons évoquées ci-dessus il permet d’étudier les relations entre des variables et de tirer des conclusions valables pour la population.
Existe-t-il une stratégie de modélisation conduisant à un « bon » modèle ? Un « bon » modèle est un modèle qui, à priori, fournit une description raisonnable. Comment parvenir à un tel modèle ? Il n’existe pas de stratégie de modélisation optimale mais des principes. On procédera par étapes : •
Tout d’abord, la modélisation (mathématique) requière le choix d’une variable dépendante dont on aimerait connaître les déterminants ou variables explicatives qui l’influencent. La nature de cette variable (quantitative, qualitative, positive ou prenant n’importe quelle valeur sur l’axe des réels, etc.) conditionnera le choix du modèle.
39
•
Sur la base de ses connaissances théoriques et de la littérature on choisira les variables explicatives.
•
Le nombre de variables que l’on peut raisonnablement introduire dans un modèle va dépendre du nombre d’observations que l’on a disposition. Une règle du pouce pour la régression logistique est d’avoir au moins 10 outcomes positifs par co-variable (3).
•
Le nombre de modèles que l’on peut formuler à partir d’un ensemble de k co-variables est quasiment infini, même pour k = 3 ou 4. En effet, avec k co-variables on peut formuler 2 k − 1 modèles contenant les variables x1 à xk. Ensuite, ce nombre est décuplé lorsqu’on introduit les interactions possibles entre les co-variables, et finalement on obtient un nombre incalculable de modèles en variant la forme fonctionnelle de la régression et de chacune des co-variables (4).
•
Sur la base de toute une batterie de tests et d’outils statistiques, le modélisateur sélectionnera un modèle parcimonieux et interprétable, maximisant ses chances de reproduire ses résultats dans une nouvelle étude (5-9).
Dans ce chapitre, nous allons présenter des outils d’aide permettant au modélisateur de sélectionner un modèle et d’évaluer sa qualité.
7.1) Le choix des co-variables Le choix des co-variables à introduire dans le modèle de régression repose non seulement sur les connaissances biomédicales, mais aussi sur la finalité du modèle. En effet, il faut distinguer entre une analyse « pronostic » et une analyse « étiologique », car la finalité n’est pas tout à fait la même. Dans une analyse « pronostic » on cherche avant tout à construire un modèle permettant de prédire (discriminer dans la régression logistique) le mieux possible les « outcomes » (0 et 1) à partir des co-variables, tandis que dans une analyse « étiologique » on s’intéresse plus particulièrement à évaluer le risque associé à un facteur. Dans ce cas le choix des facteurs confondants est primordial pour éliminer les biais autant que possible. Pour les analyses étiologiques une théorie assez sophistiquée basée sur des graphes a été développée par les chercheurs de divers disciplines (épidémiologie, économie, sciences du management, etc.). Nous n’allons pas aborder ce thème, ici, mais renvoyons le lecteur intéressé à la littérature (10,11). On retiendra que, en général, l’approche qui consiste à inclure dans le modèle toutes les covariables que l’on peut trouver peut introduire un biais et que le choix approprié des régresseurs doit se faire de manière judicieuse. La théorie des graphes nous enseigne, en particulier, qu’il ne faut pas introduire dans le modèle un facteur intermédiaire se trouvant sur le chemin de causalité entre « l’exposition » et « l’outcome ».
7.2) Le choix de la forme fonctionnelle des co-variables Une fois effectuée la sélection des variables explicatives candidates à l’analyse, il faut déterminer la forme fonctionnelle avec laquelle chaque variable continue entre dans le modèle.
40
En effet, il se peut que la transformation appropriée pour la variable « âge » par exemple soit le logarithme ou qu’il faille introduire un terme quadratique. Pour ce qui est des variables catégorielles le problème ne se pose pas puisqu’elles sont codées au moyen de variables dichotomiques 0/1. On recherchera donc la transformation la plus appropriée soit en comparant les modèles avec la variable « âge » versus « âge+âge2 » ou « âge » versus « log(âge) », soit en utilisant une méthode plus générale comme l’approche par les polynômes fractionnels ou la méthode des splines (12-14). Une méthode pour obtenir une indication sur la forme fonctionnelle à utiliser peut aussi être basée sur la définition de catégories d’âges, par exemple des quintiles, et l’estimation des coefficients associés à chacune de ces catégories. Ensuite, on représente ces coefficients en fonction du milieu des catégories d’âge sur un graphe. Remarque La méthode qui consiste à utiliser les catégories plutôt que la variable continue est utile pour rechercher la transformation adéquate, mais il ne faudrait pas remplacer la variable continue dans le modèle final par les catégories, car l’on perd non seulement de la puissance statistique, mais aussi on introduit un biais de « confounding résiduel » dans l’estimation de chacun des coefficients.
7.3) L’adéquation du modèle aux données « Goodness of fit » (*) Une fois que les étapes du choix des co-variables et de leurs formes fonctionnelles ont été effectuées, on peut déterminer la qualité de l’ajustement du modèle aux données ou, en anglais, le « Goodness of fit ». Pour fixer les idées, notons les valeurs observées de l’outcome y ' = ( y1 , y 2 , L, y n ) et les valeurs prédites par le modèle yˆ ' = ( yˆ1 , yˆ 2 , L, yˆ n ) , où n est la taille de l’échantillon. On considérera que l’ajustement est satisfaisant si : 1. La distance entre l’outcome observé y et l’outcome prédit par le modèle yˆ est petite. 2. Le modèle est bien « calibré », i.e. les fréquences prédites sont proches de celles observées. 3. Le modèle permet de bien discriminer entre les valeurs de y = 0 et y = 1 en fonction des variables explicatives x1, x2, …, xp, i.e. on obtient de bonnes sensibilités et spécificités. Pour nous aider dans cette tâche, nous allons nous appuyer sur des tests de « Goodness of fit » comme, le test de Hosmer et Lemeshow, la statistique de Pearson et la déviance, sur l’analyse des résidus comme, les résidus de déviance et de Pearson, ainsi que sur l’évaluation de la capacité à discriminer les outcomes positifs y = 1 des outcomes négatifs y = 0 par l’inspection des courbes de sensibilité et spécificité, et la courbe ROC. La démarche que nous allons adopter consiste à évaluer, d’abord, globalement l’adéquation du modèle au moyen des différents tests de « Goodness of fit », puis, en principe lorsqu’on est 41
satisfait de la qualité de l’ajustement global, à déterminer s’il n’y a pas localement des observations très mal ajustées et ayant possiblement un effet important sur l’estimation des coefficients. Le but des ces évaluations globale et locale est de s’assurer que l’ajustement du modèle soit satisfaisant pour toutes les valeurs observées dans l’échantillon des variables explicatives. Finalement, l’évaluation du pouvoir discriminant du modèle nous permettra d’appréhender si nous avons choisi les « bonnes » variables explicatives ou s’il manque d’importants régresseurs pour arriver à prédire avec suffisamment de précision l’outcome. Avant de procéder, il nous faut définir la notion de « covariate pattern » qui est essentielle dans l’analyse du « Goodness of fit » du modèle logistique. a) La notion de « covariate pattern » Supposons que notre modèle contienne p variables explicatives, que l’on notera par le vecteur x' = ( x1 , x 2 ,L , x p ) . Chaque individu dans l’échantillon aura une vecteur x ' caractérisant son âge, son sexe, etc. On appelle un « covariate pattern » une valeur caractéristique du vecteur x ' de sorte que chaque individu caractérisé par un même vecteur x ' aura le même « covariate pattern ». Ainsi, si dans un échantillon de taille n, J individus seulement ont des vecteurs x ' différents, on aura J
42
| 7 | 0.3872 | 6 | 6.5 | 13 | 12.5 | 19 | | 8 | 0.4828 | 7 | 8.2 | 12 | 10.8 | 19 | | 9 | 0.5941 | 10 | 10.3 | 9 | 8.7 | 19 | | 10 | 0.8391 | 13 | 12.8 | 5 | 5.2 | 18 | +--------------------------------------------------------+ number of observations number of groups Hosmer-Lemeshow chi2(8) Prob > chi2
= = = =
189 10 9.65 0.2904
On constate que dans notre exemple le test de Hosmer et Lemeshow est passé et que, par conséquent, l’ajustement global du modèle aux données est satisfaisant. Néanmoins, ce test est basé sur un regroupement des données en catégories et certains « covariate patterns » très mal ajustés peuvent avoir échappé. Avant d’accepter ce modèle il faut analyser les résidus pour confirmer que l’ajustement est « bon » pour tous les « covariate patterns ». Remarque La statistique de Hosmer et Lemeshow est calculée à partir de la table de contingence (g x 2) (lorsqu’on forme g groupes) des fréquences observées et espérées. Hosmer et Lemeshow ont montré que leur statistique suivait approximativement une loi du Chi2 à g-2 degrés de liberté, mais cette approximation est bonne pour autant que les fréquences espérées soient ≥ 5 , hormis une. La puissance du test est relativement faible lorsque la taille de l’échantillon est ≤ 400 . c) L’analyse des résidus Un résidu est une mesure de la distance entre l’outcome observé y et l’outcome prédit par le modèle yˆ . Comme on va le voir, il existe plusieurs définitions de résidus et chacune d’elle correspond à un concept particulier de distance. Le but de l’analyse des résidus est multiple : 1) il s’agit de vérifier qu’il n’y a pas des erreurs systématiques, 2) de déterminer s’il y a des observations très mal expliquées (résidus extrêmes) et 3) si certaines observations ont un effet important de levier sur les résultats des estimations. Comme chaque observation a son résidu associé il y a autant de résidus que d’observations. L’on considérera donc des mesures globales résumant l’ensemble des résidus par un seul chiffre et permettant ainsi d’apprécier l’ajustement global du modèle aux données (autrement dit on résume la distance entre y et yˆ ), ainsi que des mesures locales fournies par chacun des résidus et permettant de vérifier que la contribution à la mesure globale de chacune des observations est plus ou moins équivalente. Ce dernier point fait référence à la qualité de l’ajustement pour chaque « covariate pattern ». L’on désire, en effet, que le modèle ajuste bien les observations pour l’ensemble des « covariate patterns » et si certains sont très mal ajustés il faudra en déterminer, si possible, la raison. Eventuellement, certains « covariate patterns » peuvent être éliminés de l’analyse s’ils sont jugés trop « loin » du nuage de points (il s’agit d’outliers) et qu’ils ont un effet de levier important sur l’estimation des coefficients. En effet, les « outliers » peuvent avoir un effet catastrophique sur les estimations et biaiser les analyses. On cherche, en définitive, à ajuster le modèle sur le centre de gravité du nuage de points et il n’est pas désirable que quelques valeurs extrêmes (qui peuvent être des erreurs de mesure ou des cas complètement atypiques) modifient sensiblement les estimations.
43
Outliers ayant un fort effet de levier
Outliers n’ayant pas un effet de levier important Droite de régression avec outliers
y
y
Droite de régression avec outliers
Droite de régression sans outliers
x
Droite de régression sans outliers
x
figure 16 & 17
c.1) Le résidu de Pearson Nous avons noté mj le nombre d’individus avec le même « covariate pattern » x' = x ′j . De même, nous noterons le nombre de réponses positives (y = 1) au sein d’un « covariate pattern » J y j , de sorte que ∑1 y j = n1 où n1 représente le nombre total d’individus avec y = 1. Ainsi, le nombre de réponses y = 1 au sein d’un « covariate pattern » prédites par le modèle s’exprime par : yˆ j = m j Pˆ j ,
j = 1, …, J
où Pˆ j = Pˆ ( y = 1 | x j ) est la probabilité prédite par le modèle. Le résidu de Pearson est défini comme :
r ( y j , yˆ j ) = r j =
( y j − yˆ j ) m j Pˆ j (1 − Pˆ j )
.
Ce résidu sera d’autant plus grand que le nombre de cas positifs prédit est différent du nombre observé et que le dénominateur est petit. Remarque Remarquons que V ( y j | x j ) = m j Pj (1 − Pj ) , mais V ( y j − yˆ j | x j ) ≠ m j Pj (1 − Pj ) de sorte que le résidu de Pearson n’est pas standardisé (i.e. de variance constante égale à 1). La difficulté avec des résidus non standardisés est que leur amplitude dépend non seulement de la
44
différence y j − yˆ j mais aussi de sa variance, ce qui complique leur interprétation. Pour cela, on a défini le résidu de Pearson standardisé. Le résidu de Pearson standardisé est défini comme :
rs ( y j , yˆ j ) = rs j =
rj 1− hj
où h j = m j Pˆ j (1 − Pˆ j ) x ′j ( X ′VX ) −1 x j est une fonction compliquée des observations. On appelle h j le levier, car cette grandeur mesure dans un certain sens la distance du jème « covariate pattern » au centre de gravité des « covariate patterns ». Or, l’effet d’un « covariate pattern » éloigné du centre de gravité peut être important sur l’estimation des paramètres. Remarquons qu’en régression logistique l’interprétation de l’effet de levier d’un « covariate pattern » dépend de sa probabilité associée et un « covariate pattern » éloigné du centre de gravité n’est pas forcément un levier important. A partir des résidus de Pearson on définit la statistique de Pearson : X 2 = ∑ j =1 rs ( y j , yˆ j ) 2 J
Il s’agit d’une mesure globale résumant la distance entre y et yˆ . Si le nombre d’observations mj , j = 1, …, J, dans chacun des « covariate patterns » est assez grand la statistique de Pearson suit approximativement une loi du Chi2 à J-p degrés de liberté (p est le nombre de paramètres estimés dans le modèle). Lorsque le modèle contient une ou plusieurs variables continue, en général, J ≅ n de sorte que la distribution de cette statistique n’est pas une loi du Chi2 à J-p degrés de liberté et l’on retiendra que l’ajustement est d’autant meilleur que cette statistique est petite. Pour déterminer si une co-variable est mal ajustée par le modèle on peut représenter les résidus de Pearson standardisés en fonction de celle-ci. Si les résidus ne sont pas distribués de manière plus ou moins égale entre les niveaux bas et haut de cette co-variable on en déduira que l’on a un problème d’ajustement pour cette dernière. Cette démarche est illustrée dans le graphe cidessous pour un modèle expliquant la probabilité de mettre au monde en enfant de petit poids en fonction de l’âge de la maman :
45
3 Résidus standardisés de Pearson -1 0 1 2 -2 10
20
30 age of mother
40
50
figure 18
On ne constate pas de gros déséquilibre dans la distribution des résidus pour les mamans jeunes et les moins jeunes, ni de valeurs très grandes des résidus de sorte que l’on est satisfait de la qualité de l’ajustement pour cette co-variable. Remarque Si les « covariate patterns » contiennent assez d’observations les résidus sont distribués approximativement Normalement et l’on s’attend à trouver environ le 95% des points dans l’intervalle -2 à +2. c.2) Le résidu de déviance On définit le résidu de déviance comme : 1/ 2
yj m j − y j d ( y j , yˆ j ) = ± 2 y j log + (m j − y j ) log m j Pˆ j m j (1 − Pˆ j )
,
et la déviance : D = ∑ j =1 d ( y j , yˆ j ) 2 . J
Tout comme la statistique de Pearson, la déviance est une mesure globale résumant la distance entre y et yˆ . Si le nombre d’observations mj , j = 1, …, J, dans chacun des « covariate patterns » est assez grand la déviance suit approximativement une loi du Chi2 à J-p degrés de liberté (p est le nombre de paramètres estimés dans le modèle).
46
Remarque En régression logistique, maximiser la vraisemblance revient à minimiser la déviance, soit minimiser la somme des carrés des résidus de déviance. Dans ce sens, les résidus de déviance sont équivalent aux résidus du modèle linéaire.
-2
-1
Deviance Residual 0 1
2
Comme avec les résidus de Pearson, pour déterminer si une co-variable est mal ajustée par le modèle on peut représenter les résidus de déviance en fonction de cette celle-ci. En reprenant l’exemple d’avant concernant la probabilité de mettre au monde en enfant de petit poids en fonction de l’âge de la maman :
10
20
30 age of mother
40
50
figure 19
d) Détection des « covariate patterns » mal ajustés On peut repérer un « covariate pattern » mal ajusté en considérant sa contribution à la statistique du Chi2 de Pearson ou à la déviance. La réduction de la statistique de Pearson due à l’élimination des sujets du jème « covariate pattern » est donnée par : ∆X 2j = rs2j ,
tandis que pour la déviance on a :
∆D j =
d 2j (1 − h j )
.
Typiquement, un graphe de ∆X 2j (respectivement de ∆D j ) versus la probabilité prédite a l’allure suivante :
47
logistics regression dioagnostic
5
10
logistics regression dioagnostic 36
36 77
77
4
8
188 133 132 16 57 33 10
3 2
71 81 35
101 100
.2
.4
.6
154
160104 166 111 155 199 216 115 202 128 161 187 85 9396 87 143 181 146 163 121 130 141 67 201 140 113 205 170 142 224 208 150 95 107 156 103 149 127 91 176 177 179 105 135 124 125 139 206 94 109 186 106 209 210 214 97 212 123 167 189
75 8411 7883 4 5037 442313
193 192
211 145 218 200 196 219 86 92 213 225 136 203 215 184 220 174 182 168 185 222 221 159 169 190 195 151 134 223 217 131 120 226 114 191 126 204 112 129 175 173 183 207 108
.8
0
P(Y=1)
.2
154
117 116 56 46 18 43 31 69 52
148 147
71 81 35 56 119 46 137 31 69 18 43 117 116 162 99 172 52 101 100 34 62 1789 15 197 138 180 88 118 208260 42 164 148 147 160104 102 67189 140 40 5132 167 61 19 166 111 155 199 115 216 193 192 202 96 128 161 187 93 85 87 143 181 146 163 121 130 141 142 113 205 170 224 150 208 201 95 107 103 156 91 149 127 176 177 105 135 179 125 124 206 186 94 109 139 106 209 210 214 97 212 211 123 145 218 196 219 200 213 225 136 92 86 203 215 184 182 174 220 159 168 185 222 221 169 217 195 151 134 190 226 114 120 131 223 126 204 191 175 129 112 207 183 173 108
144
1
2
68 59 47 54 79 76 3045 26
133 132 68 47 59 54 79 76 3045 26
0
22 63 29 24 25 65
98
28 63 24 22 29 25 65
0
98
28
0
dx2 4
27 49
d_deviance
6
27 49 188 16 57 33 10
.4
144
119 137 162 99 172
621789 34 197 138 180 15 88 20 118 42 102 164 8260 40 5132 61 19
75 11 84 7883 4
.6
5037 442313
.8
P(Y=1)
figures 20 & 21
Les « covariate patterns » mal ajustés sont généralement représentés par des points soit en haut à droite, soit à gauche, et qui en plus se détachent sensiblement de la masse des points. Pour repérer plus facilement quels sont les sujets correspondant on a fait apparaître leur numéro d’identification. On constate sur ces figures que seuls deux « covariate patterns » sont éloignés de la masse des points et donc pas bien ajustés. Cependant, si le nombre d’observations dans les « covariate patterns » est suffisamment grand on peut montrer que ∆X 2j et ∆D j sont distribués approximativement selon une loi du Chi2 à 1 degré de liberté dont le percentile 95 est χ 02.95 (1) = 3.84 . On considérera, donc, que d’une part un « covariate pattern » est mal ajusté si la valeur de ∆X 2j ou ∆D j est plus grande que 4, d’autre part que le modèle est bien ajusté si pas plus de 5% des « covariate patterns » sont supérieurs à 4. e) Détection des points influants (effet de levier) Un résidu élevé n’implique pas forcément que le « covariate pattern » a une influence importante sur l’estimation des paramètres et l’on détermine cette influence en examinant les changements dans l’estimation des paramètres β j , j = 1,..., p (on a p covariables) lorsqu’on élimine le jème « covariate pattern ». Pour cela l’on ne va pas ré-estimer le modèle pour chaque situation mais utiliser une mesure résumant cet effet. Une mesure de l’influence d’un « covariate pattern » sur l’estimation des paramètres est fournie par :
∆βˆ j =
rs2j h j (1 − h j )
.
On remarque que cette grandeur augmente avec la magnitude du résidu standardisé de Pearson, ainsi qu’avec l’augmentation du levier. Les « covariate patterns » associés à un résidu élevé et ayant une forte influence sur l’estimation des paramètres sont a scruter attentivement.
48
Typiquement, un graphe de ∆βˆ j versus la probabilité prédite a l’allure suivante :
1
logistics regression dioagnostic
.6
133 132
98 28
.4
Pregibon dbeta
.8
188
154 77
10
0
.2
119 117 116 59 144 162 160 43138 197 22 29 65 18 20 11 99 32 202 71 56 31 187 79 172137 19 75 25 68 101 100 10262 81 8915 34 544526 17 3515546 84 88 52164 180 428260 24 47 69104 83 61 63 85 40 166 51 76 30 167 189 118 148 147 128115 111 161 163 78 45037 442313 140 67 107 170 193 192 121 130 141 210 127 206 143 9396216199 113 87 205 159 176 181 146 106 224 95 142 109 94 125 124 105 135 103 156 201 1208 50 168 97 214 209 179 177 149 91 213 86 145 211123 212 186 139 226 203 196 219 200 218 126 223 222 220 215 184 225 136 92 129 112 204 191 114 120 131 217 151 195 134 190 169 185 221 182 174 108 207 183 173 175
0
36
57 33 16 49 27
.1
.2
.3
.4
.5
.6
.7
.8
.9
1
P(Y=1)
figure 22
Dans ce graphe on constate qu’un « covariate pattern » (le no d’identification 188) a plus d’influence que les autres et est à analyser pour déterminer la raison. Néanmoins, Hosmer et Lemeshow (voir bibliographie) estiment que pour qu’un « covariate pattern » aie une influence sensible sur l’estimation des coefficients il faut que cette grandeur soit supérieure à 1, ce qui n’est pas le cas dans notre exemple et, par conséquent, nous n’en ferons pas cas. Un autre graphe résumant à la fois les « covariates patterns » mal ajustés et influants est donné par :
logistics regression dioagnostic
10
36
8
77
6
188
dx2
16 57 33 10
4
27 49 98
28 63 22 24 29 25 65
133 132
2
68 47 59 54 79 76 3045 26
154
75 8411 7883 45037
442313
0
211 200 145 218 86 196 219 213 225 136 92 203 215 184 182 174 220 159 168 185 222 221 134 190 169 131 223 217 195 151 126 204 191 226 114 120 112 173 175 129 207 183 108
81 71 35 56 119 46 137 31 69 18 43 117 116 162 99172 52 101 100 34 62 1789 15 197 138 180 88 208260 118 42 164 148 160104102 140 67189 40 5132 167 147 166 61 19 111 155 199 216 115 202 193 192 128 161 9396 187 85 87 181 146 143 163 121 130 141 113 205 170 224 142 150 208 201 95 107 156 149 103 127 176 91 179 177 124 105 135 206 109 139 125 186 94 106 209 210 214 212 97 123
144
0
.2
.4
.6
.8
P(Y=1)
figure 23
f) Evaluation du pouvoir discriminant du modèle : sensibilité, spécificité et courbe ROC On utilise le modèle Logistique pour modéliser la probabilité des attributs 0/1 de la variable dépendante y en fonction des co-variables x1, x2, …, xp. A partir des probabilités estimées on décidera en fixant un seuil, par exemple à 0.5, de classer l’individu dans la catégorie y = 1 si sa probabilité est supérieure au seuil et dans la catégorie y = 0 sinon. Il s’agit d’une règle de classement :
Si
Pˆ j = Pˆ ( y = 1 | x j ) ≥ seuil
alors
yˆ = 1 et 0 sinon.
49
Il est intéressant de déterminer la performance du classement et comment celui-ci dépend du seuil (ou de la règle) choisi. Pour cela nous allons considérer les notions de sensibilité et spécificité. La sensibilité est définie comme la probabilité de classer l’individu dans la catégorie y = 1 (on dit que le test est positif) étant donné qu’il est effectivement observé dans celle-ci : Sensibilité = P(test + | y = 1) La spécificité est définie comme la probabilité de classer l’individu dans la catégorie y = 0 (on dit que le test est négatif) étant donné qu’il est effectivement observé dans celle-ci : Spécificité = P(test - | y = 0) Voici un exemple fictif de classement obtenu pour un seuil choisi à priori où l’on a calculé la sensibilité et la spécificité :
Classement
observé
total
y=1
y=0
y=+
5
0
5
y=-
142
428
570
total
147
428
575
Sensibilité = 5/147 = 3.4% ; Spécificité = 428/428 = 100%
0.00
Sensitivity/Specificity 0.25 0.50 0.75
1.00
Lorsqu’on varie le seuil (en anglais « cutoff ») la sensibilité et la spécificité changent, puisque la règle de classement est modifiée. Afin de représenter les valeurs pour toutes les possibilités de seuils on dessine sur un graphe des courbes de sensibilités et spécificités :
0.00
0.25
0.50 Probability cutoff Sensitivity
0.75
1.00
Specificity
figure 24
50
On constate que, dans cet exemple, en fixant le seuil à 0.3 on obtient un classement avec une sensibilité et spécificité d’environ 70%.
0.00
0.25
Sensitivity 0.50
0.75
1.00
Comme indicateur de la capacité du modèle à discriminer on utilisera la courbe ROC :
0.00
0.25
0.50 1 - Specificity
0.75
1.00
Area under ROC curve = 0.7462
figure 25
La surface sous cette courbe nous permet d’évaluer la précision du modèle pour discriminer les outcomes positifs y = 1 des outcomes négatifs y = 0. On retiendra comme règle du pouce : Si
aire ROC = 0.5
il n’y a pas de discrimination
Si
aire 0.7 ≤ ROC < 0.8 la discrimination est acceptable
Si
aire 0.8 ≤ ROC < 0.9 la discrimination est excellente
Si
aire ROC ≥ 0.9
la discrimination est exceptionnelle
Remarque Un modèle mal ajusté, i.e. mal calibré, peut très bien fournir une bonne discrimination. Pour s’en convaincre, il suffit de penser à la situation où l’on ajoute 0.15 à toutes les probabilités estimées. Le modèle sera alors mal calibré, mais en déplaçant le seuil de 0.15 on obtiendra la même discrimination. Un bon modèle doit être bien calibré et permettre une bonne discrimination. g) La validation du modèle Les paramètres de notre modèle de régression logistique sont estimés à partir de l’échantillon de données dont on dispose. On peut se poser la question de la validité de la transposition de nos résultats à l’ensemble de la population d’intérêt. Pour que l’on puisse inférer des caractéristiques de la population à partir de notre échantillon il faut, bien entendu, qu’il soit représentatif de celle-ci. Cette représentativité s’obtient en particulier lorsqu’il y a eu tirage aléatoire.
51
Maintenant, nos estimations sont sensibles à la taille de notre échantillon et il est légitime de se poser la question de la reproductibilité de nos résultats. En particulier, est-ce qu’on obtiendra a peu près les mêmes estimations si l’on analyse un autre échantillon de même taille et issu de la même population ? Comme en général l’on ne dispose pas d’un deuxième échantillon de la population, pour tenter de répondre à cette dernière question des techniques de validation basées sur le reéchantillonnage ou la partition de l’échantillon ont été développées. Pour ne pas trop allonger le texte, nous n’allons pas rentrer dans la description de ces techniques, ici, mais l’on retiendra que l’évaluation de l’ajustement de notre modèle sur un autre échantillon provenant de la même population fournira très probablement des résultats moins optimistes. En définitive, on retiendra que plus l’échantillon de données disponibles sera grand plus les estimations seront fiables et reproductibles. Règle du pouce que nous avions énoncée auparavant dans le texte de « dix outcomes positifs par co-variable introduite dans le modèle » va dans le sens d’une meilleure reproductibilité des résultats.
7.4) Limitations et biais (*) Nous allons ci-dessous passer en revue quelques problèmes courrant que l’on rencontre lorsqu’on fait de la régression logistique. a) Le problème de la séparabilité ou quasi-séparabilité (*) Un problème que l’on rencontre parfois lorsqu’on estime un modèle de régression logistique est celui de la séparabilité ou quasi-séparabilité, dans quel cas les estimateurs ne convergent pas et, en principe, le logiciel affiche un message d’erreur. Dans ce cas les résultats ne sont pas fiables et il faut rechercher la cause de la séparabilité. STATA possède un algorithme permettant de détecter le problème de séparabilité avant de lancer la routine de maximisation de la vraisemblance. Pour que le modèle de régression logistique soit estimable il faut qu’il y ait une superposition des valeurs des co-variables pour les différents outcomes 0/1. Autrement dit, il ne faut pas que les co-variables discriminent parfaitement. La situation la plus simple de séparabilité est celle où pour une covariable qualitative il y a une fréquence nulle dans une table de contingence :
Outcome
1
2
3
total
1
7
12
20
39
0
13
8
0
21
total
20
20
20
60
Odds ratio
1
2.79
inf
52
Dans ce cas, tous les individus dans la catégorie 3 de la variable explicative ont l’outcome 1 et l’Odds ratio associé est infini. Dans ce cas, ces 20 individus sont éliminés de l’estimation et l’on rapportera simplement que les individus de cette catégories ont probabilié 1 d’avoir l’oucome positif. La situation impliquant une variable continue peut être représentée dans les graphes suivants :
Maladie coronarienne
oui
Maladie coronarienne
oui
Sans ce point on a séparabilité complète
non
Ici on a séparabilité Quasi-complète
non age
50
age
50
figures 26 & 27
Si l’on n’observait dans notre échantillon pas de patient sain âgé de plus de 50, on aurait un problème de séparabilité complète. Autrement dit, les patients de plus de 50 ans auraient une probabilité 1 d’avoir un outcome positif et la discrimination en fonction de l’âge serait alors parfaite. On a séparabilité est quasi-complète si le seul chevauchement observé est exactement pour l’âge de 50 ans. La situation de séparabilité ou quasi-séparabilité impliquant plusieurs co-variables catégorielles ou continue est plus complexe et nous renvoyons le lecteur intéressé à la littérature (15). b) Le problème de « l’overfitting » Le problème de « l’overfitting » intervient lorsqu’on a sélectionné beaucoup de co-variables dans le modèle ou lorsqu’on n’a pas assez d’outcomes positifs par rapport au nombre de covariables retenues. Dans ce cas, les résultats ne sont pas reproductibles et les associations fortuites. Le risque « d’overfitting » augmente, en particulier, lorsqu’on introduit plusieurs termes d’interaction dans le modèle. Le problème est d’autant plus délicat que les variables explicatives sont catégorielles polytomiques, dans quel cas il y a foison de termes d’interaction à estimer. Dans ce cas, il vaut mieux se limiter aux interaction faisant intervenir un assez grand nombre de patients. c) Le biais de sélection Un biais de sélection se produit lorsque l’échantillon de données n’est pas tiré aléatoirement ou qu’il n’y a pas eu randomisation. Dans ce cas, il se peut que les patients exposés au
53
traitement diffèrent systématiquement de ceux n’ayant pas été exposés et l’analyse de l’effet de l’exposition est biaisé par la sélection. Pour corriger le biais de sélection on peut recourir aux modèles de sélection et la littérature sur le thème est vaste et assez complexe. On distingue aussi d’autres biais potentiels, comme le biais de classement qui résulte des erreurs de mesures et le biais de confusion provenant de la non-inclusion dans le modèle d’un facteur confondant. d) Le problème de surdispersion « overdispersion » Dans la régression logistique on fait l’hypothèse que le modèle binomial est adéquat pour décrire la variabilité aléatoire de l’outcome. Or, il se peut que ce ne soit pas le cas et que la dispersion soit en réalité supérieure à celle prédite par le modèle. On parle alors de surdispersion ou, en anglais, « d’overdispersion ». Dans ce cas, les p-value et intervalles de confiances estimés par la régression logistique ordinaire sont biaisés. e) Extensions e.1) Le cas de données répétées Lorsque les données sont répétées, comme lorsqu’on suit une cohorte de patients au fil du temps, elles ne sont pas indépendantes et il faut utiliser des techniques statistiques appropriées comme des modèles marginaux (e.g. GEE) ou mixte (à coefficients aléatoires). e.2) Le cas de données agrégées « cluster » Comme pour les données répétées, les données agrégées ne sont pas indépendante et il faut utiliser des techniques statistiques qui prennent en compte la dépendance des données, e.g. les modèles multi-niveaux.
Exercice 7 Dans cet exercice nous vous proposons de considérer les données « Low birth weight » et d’appliquer les notions que nous venons de voir pour ajuster un modèle logistique à ces données.
54
8) Le logiciel statistique STATA Nous proposons ci-dessous une liste de commandes STATA utiles pour la régression logistique. Pour se fixer les idées, nous prendrons comme exemple de variables celles utilisées par Hosmer et Lemeshow dans le fichier de données « Low birth weight data » disponible sur le web (http://www.stata-press.com/data/r8/lbw). La commande pour invoquer la régression logistique est « logistic » ou « logit », la première produisant les résultats sous forme d’Odds ratio et la deuxième fournissant les coefficients estimés du modèle : * Chargement des données Low birth weight ***************************************** cd "C:\Mes documents\Cours régression logistique\" use "C:\Mes documents\Cours régression logistique\lbw.dta", clear xi: logistic low age lwt i.race smoke ptl ht ui
Pour représenter les données, ainsi que la forme de la relation « non ajustée » entre la probabilité de « low birth weight » et l’âge on utilisera la commande : twoway scatter low age, ylabel(0(0.2)1) ytitle(P(low)) title(Relation entre low birth weight et l'âge) || lowess low age, sort
Pour déterminer et représenter la relation fonctionnelle entre l’âge et « low birth weight » on peut utiliser les polynômes fractionnels : xi: fracpoly logistic low age lwt i.race smoke ptl ht ui fracplot age
On obtient toutes sortes de mesures d’adéquation du modèle aux données avec les commandes: fitstat lfit
Pour le test de Hosmer et Lemeshow : lfit, group(10) table
Pour déterminer la sensibilité et la spécificité: Lstat lsens
Pour la courbe ROC: lroc
Pour calculer les probabilités prédites par le modèle:
55
predict p
Pour calculer les résidus standardisés de Pearson: predict rstd, rs label var rstd "Standardized Residual"
Pour représenter les résidus standardisés de Pearson en fonction de la probabilité sur un graphe, où les points sont en bleu et étiquetés par le numéro d’identification et avec une ligne horizontale passant par 0 : graph twoway (scatter rstd p, msymbol(smcircle) msize(small) mcolor(red) mlabel(id) /// mlabsize(vsmall)), ytitle(Standardized Pearson Residual, size(small)) yscale(titlegap(2) /// range(-5 +5) outergap(3)) ylabel(-5(1)5, labgap(2) labsize(vsmall)) xtitle(P(Y=1), size(small)) /// xscale(titlegap(2) range(0 1) outergap(3)) xlabel(0(0.1)1, labgap(2) labsize(vsmall)) /// yline(0) title(Graph des résidus de Pearson standardisés, size(medsmall))
Pour calculer le résidu de déviance : predict deviance, deviance label var deviance "Deviance Residual"
Pour calculer la réduction de la statistique du Chi2 de Pearson, après avoir enlevé les sujets de covariate pattern xj : predict dx2, dx2
et représenter ceci dans un graphe: graph twoway (scatter dx2 p, msymbol(smcircle) msize(small) mcolor(red) mlabel(id) mlabsize(vsmall)), /// ytitle(dx2, size(medsmall)) yscale(titlegap(2) outergap(3)) xtitle(P(Y=1), size(small)) /// xscale(titlegap(2) outergap(3)) xlabel(, labgap(2) labsize(vsmall)) /// title(logistics regression dioagnostic, size(small)) xline(0.5)
Pour calculer la réduction de la déviance, après avoir enlevé les sujets de covariate pattern xj : predict dd, ddeviance
et représenter ceci dans un graphe: graph twoway (scatter dd p, msymbol(smcircle) msize(small) mcolor(green) mlabel(id) mlabsize(vsmall)), /// ytitle(d_deviance, size(medsmall)) yscale(titlegap(2) outergap(3)) xtitle(P(Y=1), size(small)) /// xscale(titlegap(2) outergap(3)) xlabel(, labgap(2) labsize(vsmall)) /// title(logistics regression dioagnostic, size(small)) xline(0.5)
Pour calculer la statistique d’influence d’un covariate pattern sur l’estimation des parameters: predict cook, dbeta label var cook "Pregibon dbeta"
et représenter ceci dans un graphe: 56
graph twoway (scatter cook p, msymbol(smcircle) msize(small) mcolor("0 0 255") mlabel(id) mlabsize(vsmall) mlabposition(12)), /// ytitle(Pregibon dbeta, size(medsmall)) yscale(titlegap(2) outergap(3)) xtitle(P(Y=1), size(small)) /// xscale(titlegap(2) range(0 1) outergap(3)) xlabel(0(0.1)1, labgap(2) labsize(vsmall)) /// title(logistics regression dioagnostic, size(small)) xline(0.5)
Le graphe préféré de Hosmer et Lemeshow : graph twoway (scatter dx2 p [weight=cook], msymbol(Oh)) || (scatter dx2 p, mlabel(id) mlabsize(vsmall) mlabposition(12)), /// ytitle(dx2, size(medsmall)) yscale(titlegap(2) outergap(3)) xtitle(P(Y=1), size(small)) /// xscale(titlegap(2) outergap(3)) xlabel(, labgap(2) labsize(vsmall)) /// title(logistics regression dioagnostic, size(small)) xline(0.5) legend(off)
57
58
Bibliographie Livres : •
Hosmer DW and Lemeshow S. Applied logistic regression. John Wiley & Son 1989.
•
Amemiya T. Advanced Econometrics. Harvard University Press 1985.
•
Maddala GS. Limited dependent and qualitative variables in econometrics. Cambrige University Press 1983.
Articles: 1. Chunrong A. and Norton EC. Interaction terms in logit and probit models. Economic letters 2003;80:123-129. 2. Zhang J and Yu KF. What’s the relative risk. JAMA 1998;280:1690-1691. 3. Peduzzi P, Concato J, Kemper E et al. A simulation study of the number of events per variable in logistic regression analysis. Journal of Clinical Epidemiology 1996;49:1373-1379. 4. Greenland Sr. Modeling and variable selection in epidemiologic analysis. American Journal of Public Health 1989;79:340-349. 5. Concato J, Feinstein AR and Holford T. The risk of determining risk with multivariable models. Annals of Internal Medicine 1993;118:201-210. 6. Greenland S and Neutra R. Control of confounding in the assessment of medical technology. International Journal of Epidemiology 1980;9:361-367. 7. Katz MH. Multivariable analysis: a primer for readers of medical research. Annals of Internal Medicine 2003;138:644-650. 8. Bagley SC, White H and Golomb BA. Logistic regression in the medical literature: standards for use and reporting, with particular attention to one medical domain. Journal of Clinical Epidemiology 2001;54:979-985. 9. Traissac P, Martin-Prével Y, Delpeuch F et al. Régression logistique vs autres modèles linéaires généralisés pour l’estimation de rapports de prévalences. Revue d’Epidémiologie et de Santé publique 1999 ;47 :593-604. 10. Hernan MA, Hernadez-Diaz S, Werler MM et al. Causal knowledge as a prerequisite for confounding evaluation : an application to birth defects epidemiology. American Journal of Epidemiology 2002;155;176-184. 11. Greenland S, Pearl J and Robins J. Causal diagrams for epidemiologic research. Epidemiology 1999;10:37-48. 12. Greenland S. Dose-response and trend analysis in epidemiology: alternatives to categorical analysis. Epidemiology 1995;6:356-365. 13. Royston P, Altman DG: Regression using fractional polynomials of continuous covariates: parsimonious parametric modelling. Appl.Stat. 1994;43:429-467.
59
14. Royston P, Ambler G, Sauerbrei W: The use of fractional polynomials to model continuous risk variables in epidemiology. Int.J.Epidemiol. 1999;28:964-974. 15. Albert A and Anderson JA. On the existence of maximum likelihood estimates in logistic regression models. Biometrika 1984;71:1-10.
Pour l’utilisation de STATA se référer aux manuels suivants : •
STATA Reference Manual de A-Z
•
Getting Started with Stata for Windows
•
Stata User’s Guide
•
Stata Graphics
On consultera aussi l’ouvrage de Rabe-Hesketh et Everitt : •
A Handbook of Statistical Analysis using Stata
60
