Introduction à l'analyse de covariance (ANCOVA)

Introduction à l'analyse de covariance avec R

Table des matières

Introduction

L’analyse de covariance ou ANCOVA est un modèle de régression linéaire dans lequel une réponse quantitative est expliquée en fonction de deux variables :

  • l’une quantitative ( que l’on nomme parfois co-variable)
  • l’autre qualitative, (à plusieurs modalités).

Un modèle ANCOVA pourrait, par exemple, être employé pour expliquer la fréquence cardiaque maximale de sujets (réponse), en fonction de leur âge (variable quantitative) et de leur statut par rapport à une pathologie cardiaque (variable qualitative : malade ou non malade)

L’analyse de covariance (ANCOVA) peut être envisagée sous deux angles. Le premier est celui d’une régression linéaire simple par sous-groupe. L’ANCOVA peut alors nous permettre de répondre à la question “est ce que la relation linéaire entre la réponse et la variable quantitative est différente entre les sous-groupes (modalités de la variable qualitative) ? 
Par exemple, ”est-ce que la relation linéaire entre la fréquence cardiaque maximale et l’âge, est différente en présence ou en absence de pathologie cardiaque ?“. Exprimée autrement :”est ce que la relation linéaire entre la fréquence cardiaque maximale et l’âge dépend du statut vis-à-vis de la pathologie (malade ou non malade)”

Lorsque l’ANCOVA est envisagée sous cet angle de la régression linéaire, on représente, assez intuitivement, les données sous formes de droites, comme ceci.

analyse de covariance : plusieurs droites différentes

Et comme on le verra plus loin, différentes droites peuvent être envisagées.

 

Le deuxième angle, sous lequel on peut envisager l’ANCOVA, c’est celui de la comparaison de moyennes prédites de la réponse, entre les sous-groupes (modalités de la variable qualitative). Les moyennes étant prédites par le modèle, elles sont ajustées sur la variable quantitative. Dans cette situation, l’ANCOVA va nous permettre de répondre à la question : “est ce que les fréquences cardiaques maximales des sujets malades et non malades sont différentes, une fois l’effet de l’âge des sujets pris en compte” ? Formulé autrement : “est ce que le niveau moyen de la fréquence cardiaque maximale, dépend du statut pathologique (malade / pas malade), une fois que l’effet de l’âge a été retiré.

Dans cette situation, une des représentations graphique qui me semble le plus adéquat est la suivante : 

Analyse de covariance : Représentation conjointe des données observées et des moyennes prédites par le modèle ancova

Les croix représentent les fréquences cardiaques observées des sujets (en rouge pour les sujets non malades, en bleu pour les malades). Les ronds noirs représentent les fréquences cardiaques maximales moyennes, prédites par le modèle. On les appelle aussi parfois “moyennes conditionnelles” ou “moyennes marginales”. Ce sont les moyennes prédites par le modèle ANCOVA, une fois que l’effet de l’âge a été pris en compte. Les enveloppes grisées représentent l’intervalle de confiance à 95% de ces moyennes prédites.

L’ ANCOVA comme une régression par sous-groupe

Lorsque l’ANCOVA est employée pour répondre à la question “est ce que la relation linéaire entre la réponse et la covariable est différente entre les sous-groupes”, plusieurs sous-questions peuvent être posées et explorées en ajustant des modèles de régressions linéaires différents.

La première sous-question que l’on va se poser est “est ce que l’évolution de la réponse en fonction de la co-variable est différente selon les sous-groupes (modalités de la variable qualitative)”. Derrière le terme “évolution”, il y a la notion de pente. On cherche donc à répondre à la question “est ce que les pentes des droites des sous-groupes sont différentes ?

Pour cela, on ajuste un modèle spécifiant k droites différentes, avec k pentes non contraintes (elles peuvent être différentes). Pour permettre l’ajustement de ces k pentes non contraintes, il est nécessaire d’inclure un terme d’interaction entre la co-variable quantitative et la variable qualitative. On parle de modèle complet :

# data : dataset heart_disease du package funModeling, limité aux hommes
library(funModeling) 
mydata <- heart_disease %>% 
filter(gender=="male")

# ajustement du model complet
lm_full <- lm(max_heart_rate~age + 
                has_heart_disease + 
                age +
                has_heart_disease:age,
                data = mydata) 
# représentation graphique
ggplot(mydata, aes(y=max_heart_rate, x=age, colour=has_heart_disease))+
    geom_point()+
    geom_smooth(method="lm") 
analyse de covariance : plusieurs droites différentes

Les estimations des pentes (slopes)  et des ordonnées à l’origine  (intercepts) peuvent être obtenues avec la fonction summary(), comme ceci :

library(broom)

tidy(summary(lm_full)) 
## # A tibble: 4 x 5
##   term                     estimate std.error statistic  p.value
##   <chr>                       <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)               221.       12.6       17.6  1.12e-42
## 2 age                        -1.17      0.243     -4.80 3.07e- 6
## 3 has_heart_diseaseyes      -56.2      17.9       -3.15 1.90e- 3
## 4 age:has_heart_diseaseyes    0.691     0.330      2.09 3.78e- 2 

L’estimate de la ligne “intercept” correspond à l’ordonnée à l’origine du sous-groupe de référence (ici pas de pathologie – no). Ainsi, la fréquence cardiaque maximale moyenne prédite pour un âge = 0, chez les patients sans pathologie cardiaque est de 221 bpm

L’estimate de la ligne “age” correspond à la pente de l’évolution de la fréquence cardiaque maximale en fonction de l’âge, chez les patients sans pathologie cardiaque (groupe de référence). Ainsi, la fréquence cardiaque maximale diminue de 1.17 bpm par an.

L’estimate de la ligne “has_heart_diseaseyes” correspond à la différence des ordonnées à l’origine des patients avec et sans pathologie. Ainsi la fréquence cardiaque maximale moyenne prédites, à 0 ans, des sujets avec pathologie cardiaque est de 221-56.24 = 164.76 bpm.

L’estimate de la ligne “age:has_heart_diseaseyes” correspond à la différence de pente entre des sujets avec et sans pathologie. Ainsi la pente des sujets avec pathologie cardiaque est : -1.17+0.69 = -0.48 bmp/an. Environ deux fois moins que les sujets sans pathologie.

 

La p-value de l’effet global de l’interaction peut être obtenue comme ceci :


Anova(lm_full)
## Anova Table (Type III tests)
## 
## Response: max_heart_rate
##                       Sum Sq  Df  F value    Pr(>F)    
## (Intercept)           123684   1 310.1207 < 2.2e-16 ***
## age                     9192   1  23.0466 3.073e-06 ***
## has_heart_disease       3950   1   9.9050  0.001898 ** 
## age:has_heart_disease   1743   1   4.3707  0.037813 *  
## Residuals              80563 202                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Ici, l’interaction est significative (p=0.037), l’évolution de la fréquence cardiaque maximale (la pente) est différente chez les sujets malades et non malades.

Lorsque aucune différence d’évolution n’est mise en évidence, autrement dit que rien ne nous permet d’affirmer que les pentes sont différentes, dans un second temps, nous pouvons poursuivre l’analyse et explorer si les ordonnées à l’origine sont différentes.

En pratique, il s’agit de réaliser deux ajustements : l’un avec k droites de pentes identiques, et un autre avec une seule droite (pentes et ordonnées identiques pour les k droite). Puis de comparer ces ajustements avec un test F.

Si le test est significatif, cela signifie que les ordonnées à l’origine sont différentes et donc que les relations entre la réponse et la covariable sont différentes, mais avec des pentes identiques.

Si le test n’est pas significatif, cela signifie que rien ne nous permet d’affirmer que les ordonnées à l’origine sont différentes (en plus des pentes que rien ne nous permettait d’affirmer qu’elles étaient différentes). Ainsi,  rien ne nous permet d’affirmer que les relations entre la réponse et la covariable sont différentes dans les sous-groupes. On supposera alors qu’il s’agit d’une relation indentique (une seule droite)

Pour ajuster un modèle avec des pentes identiques, il suffit de retirer le terme d’interaction :

lm_same_slope <- lm(max_heart_rate ~ age + has_heart_disease, data=mydata) 

Pour obtenir l’estimation des paramètres :

tidy(summary(lm_same_slope))
## # A tibble: 3 x 5
##   term                 estimate std.error statistic  p.value
##   <chr>                   <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)           202.        8.73      23.2  6.60e-59
## 2 age                    -0.792     0.166     -4.77 3.46e- 6
## 3 has_heart_diseaseyes  -19.4       2.94      -6.58 3.83e-1 

Voici une visualisation de cette modélisation :

library(moderndive) # pour geom_parallel_slopes()
ggplot(mydata, aes(x=age, y=max_heart_rate, colour=has_heart_disease))+
    geom_point()+
    geom_parallel_slopes() 
Analyse de covariance : plusieurs droites avec une pente commune

L’ajustement du modèle avec une droite de régression linéaire unique se réalise en retirant la variable qualitative du modèle : 

# ajustement du modèle
lm_one_line <-  lm(max_heart_rate~has_heart_disease,data=mydata)

# affichage de l'estimation des paramètres
tidy(summary(lm_one_line ))
## # A tibble: 2 x 5
##   term                 estimate std.error statistic   p.value
##   <chr>                   <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)             162.       2.21     73.3  1.82e-148
## 2 has_heart_diseaseyes    -23.4      2.97     -7.87 1.99e- 13 

La visualisation de ce modèle : 

ggplot(mydata, aes(x=age, y=max_heart_rate))+
    geom_point(aes(colour=has_heart_disease))+
    geom_smooth(method="lm") 
analyse de covariance : modélisation d'une seule droite
# comparaison des ajustements
anova(lm_same_slope, lm_one_line)
## Analysis of Variance Table
## 
## Model 1: max_heart_rate ~ age + has_heart_disease
## Model 2: max_heart_rate ~ has_heart_disease
##   Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
## 1    203 82306                                  
## 2    204 91545 -1   -9239.2 22.788 3.457e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Ici, le test est significatif (pval < 0.001) , cela signifie qu’une différence entre les ordonnées à l’origine est mise en évidence.

Si le test n’est pas significatif alors cela signifie que rien ne nous permet d’affirmer que la relation entre la réponse et la covariable est différente selon les sous-groupes. On supposera alors une seule et même relation (une droite unique)

L'ANCOVA comme une comparaison de moyennes ajustées

Dans cette situation, on commence, là aussi, par ajuster le modèle ANCOVA complet :

# ajustement du model complet
lm_full <- lm(max_heart_rate~age + 
                has_heart_disease + 
                age +
                has_heart_disease:age,
                data = mydata)

Anova(lm_full)
## Anova Table (Type III tests)
## 
## Response: max_heart_rate
##                       Sum Sq  Df  F value    Pr(>F)    
## (Intercept)           123684   1 310.1207 < 2.2e-16 ***
## age                     9192   1  23.0466 3.073e-06 ***
## has_heart_disease       3950   1   9.9050  0.001898 ** 
## age:has_heart_disease   1743   1   4.3707  0.037813 *  
## Residuals              80563 202                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Nous pouvons obtenir une estimation des moyennes prédites ajustées (ou marginales) à l’aide de la fonction emmeans() du package emmeans, en utilisant l’argument specs=pairwise~has_heart_disease

library(emmeans)
emmeans(lm_full1, specs=pairwise~has_heart_disease)
## $emmeans
##  has_heart_disease emmean   SE  df lower.CL upper.CL
##  no                   159 2.19 202      154      163
##  yes                  139 1.94 202      136      143
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE  df t.ratio p.value
##  no - yes     19.1 2.92 202   6.518  <.0001 

La moyenne de la fréquence cardiaque maximale marginale (une fois l’effet de l’âge retiré), prédite par le modèle est de 159 bpm pour les sujets sans pathologie et 139 bpm pour les sujets ayant une pathologie cardiaque.

Le test statistique met en évidence une différence entre ces deux moyennes dans le sens d’une moyenne prédite plus faible pour les sujets avec pathologie cardiaque (pval < 0.0001).

Une visualisation peut être obtenue avec la fonction emmip() :

emmip(lm_full1,  ~ has_heart_disease, CIs=TRUE, linearg=NULL) + 
    geom_jitter(aes(x = has_heart_disease, 
                    y = max_heart_rate, 
                    colour = has_heart_disease),
              data = mydata, pch = 4, width = 0.1)+
     labs(y = "Estimated marginal mean (max_heart_rate)", y = "Heart_disease")+
    theme_bw()+
    theme(legend.position = "none") 
Analyse de covariance : Représentation conjointe des données observées et des moyennes prédites par le modèle ancova

Si l’interaction est non significative, on peut la retirer du modèle pour ajuster le modèle avec des droites parallèles, puis estimer et réaliser les comparaisons de moyennes à partir de ce modèle plus restreint.

Conclusion

Dans cet article, nous avons vu les deux grandes situations dans lesquelles un modèle ANCOVA peut être employé : pour évaluer si la relation linéaire entre une réponse et une covariable est différente dans des sous-sous-groupes (cela peut être très utile pour explorer les données), ou encore pour comparer des moyennes après prise en compte d’une covariable.

Je me suis ici focalisée sur le principe de l’ANCOVA, sans jamais vérifier les hypothèses du modèle linéaire. Je réaliserai prochainement un tutoriel plus complet dans lequel je vous montrerai la procédure que j’utilise, étape par étape, pour réaliser une ANCOVA, en vérifiant les hypothèses, en linéarisant les relations si nécessaires, en vous parlant aussi des différentes paramétrisations.

D’ici là, si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee.

Commentaires

4 réponses

    1. Bonjour Gabriel,

      non il n’est pas possible d’avoir l’article en pdf, il est uniquement consultables sur le site.
      Bonne continuation

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée. Les champs obligatoires sont indiqués avec *

Bonjour !

vous venez souvent ?

Identifiez-vous pour avoir accès à toutes les fontionnalités !

Aide mémoire off'R ;)

Enregistrez vous pour recevoir gratuitement mes fiches “aide mémoire” (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.