Introduction à l'analyse de covariance (ANCOVA)
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.
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 :
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") 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() 
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 :
La visualisation de ce modèle :
# ajustement du modèle
lm_one_line <- lm(max_heart_rate~age,data=mydata)
# affichage de l'estimation des paramètres
## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 208. 9.54 21.8 2.89e-55
## 2 age -1.10 0.175 -6.31 1.71e- 9 ggplot(mydata, aes(x=age, y=max_heart_rate))+
geom_point(aes(colour=has_heart_disease))+
geom_smooth(method="lm") 
# 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_full, 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_full, ~ 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") 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.
14 réponses
Merci pour ce travail,
Est-il possible de l’avoir en copie PDF?
Merci
Bonjour Gabriel,
non il n’est pas possible d’avoir l’article en pdf, il est uniquement consultables sur le site.
Bonne continuation
Bonjour Claire,
Avez vous un code qui permet de tracer un graphe à trois axes sur r
Merci
Omar
Bonjour Omar,
je ne vois pas bien de quoi vous parlez, est-ce que vous pouvez, s’il vous plait, compléter votre message avec un lien vers une page internet illustrant votre demande ?
Merci
Excellente
Bonjour madame je suis très ravie de ce travail je voulais avoir des codes sur r
Bonsoir madame merci beaucoup pour le travail, je suis ravi de l’avoir reçu. S’il vous plaît est-ce possible d’avoir le pdf
Bonjour Odilon,
aucune version pdf n’est disponible, les articles sont seulement consultables sur le site.
Bonne continuation
Claire Della vedova
Excellent article. Nous espérons connaître encore plus avec vous. 🙏🙏
Dans l’écriture du modèle : lm_full <- lm(…) , la variable "age" apparaît deux fois. C"est volontaire ?
Bonjour,
Vous écrivez “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” mais dans le modèle qui suit vous mettez la fréquence cardiaque max et fonction du groupe (pathologie oui ou non) et non pas en fonction de l’age. Donc vous retirez finalement la variable quantitative.
Quelle est donc la bonne formule/formulation ?
Merci,
Clara
Bonjour,
Vous avez raison, c’est une coquille !
Voici la correction
lm_one_line <- lm(max_heart_rate~age,data=mydata) tidy(summary(lm_one_line )) ## # A tibble: 2 × 5 ## term estimate std.error statistic p.value ## ## 1 (Intercept) 208. 9.54 21.8 2.89e-55
## 2 age -1.10 0.175 -6.31 1.71e- 9
Merci.
Bonjour,
Merci pour le partage de ce travail.
Dans le code vous utilisez le modele .._full1, moi j’utilise le modele .._full et j’obtiens les memes résultats. S’agit-il d’une erreur?
Par ailleurs vous aviez promis de communiquer un travail plus détaillé d’ANCOVA (avec la verification des hypothèses notamment)
On espere le voir bientôt
Bonjour,
oui le full1 était une coquille, c’est full tout court.
Merci pour la piqure de rappel concernant le tuto, je le fait remonter dans ma litse, mais le temps me manque….
Bonen continuation