Il est très facile de faire une ANCOVA en langage R, mais, avant de pouvoir tirer des conclusions, il faut vérifier que les assomptions (hypothèses de base de l'ANCOVA) sont vérifiées. Ces assomptions sont :
Indépendance des observations
Normalité des résidus du modèle complet
Homogénéité des variances entre niveaux du facteur
Linéarité de la relation covariable(s) – Y (variable dépendante)
Homogénéité des pentes (parallélisme)
En plus de ces assomptions académiques, il faut aussi vérifier l'absence d'effet de levier liées à des valeurs aberrantes :
Vérifier la présence d'outliers (valeurs aberrantes, potentiellement extrêmes)
Ce faisant, si une assomption n'est pas respectée, que peut-on faire : la dernière partie de cette page présente les alternatives :
Alternatives à l'ANCOVA
En fin d'étude, comme pour l'ANOVA, on doit envisager de faire des tests post-hocs :
Tests post-hocs post-ANCOVA
Remarque extrêmement utile : la fonction m.test() de mon package {KefiR} réalise automatiquement tout ce que présente cette page : elle établit le modèle ANCOVA, testes les assomptions, bascule sur les alternatives robustes si nécessaires et va jusqu'à faire les tests post-hocs idoines. Je vous invite à l'essayer.
# Création des données
groupe <- c(rep("A", 10), rep("B", 10))
age <- c(20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28)
performance <- c(80, 85, 90, 95, 110, 105, 110, 111, 120, 125,
70, 75, 80, 85, 72, 95, 100, 105, 110, 115)
donnees <- data.frame(groupe = groupe,
age = age,
performance = performance)
# Réalisation de l'ANCOVA
resultat <- aov(performance ~ age + groupe, # ATTENTION de mettre age (covariable avant groupe)
data = donnees)
# Remarque : si on reste sur l'aov() par défaut, il faut toujours mettre la covariable en premier pour supprimer son effet.
# Si en revanche, on réalise une ANCOVA de type 3, nous n'avons plus à nous soucier de l'ordre des facteurs et covariables.
# Affichage des résultats
summary(resultat)
Il suffit d'installer {KefiR}
Puis...
m.test(performance ~ age + groupe, data=donnees)
Et tous les contrôles d'assomptions et alternatives éventuelles seront appliquées.
Dans le cadre de la construction d'un modèle ANCOVA, il faut vérifier qu'il n'y a pas de biais dans la collecte des données.
chaque ligne correspond à une observation distincte,
aucun sujet n’apparaît plusieurs fois (pas de mesures répétées non modélisées),
les observations ne sont pas structurées en groupes, lots ou clusters non pris en compte dans la formule (p. ex. classes, parcelles, familles),
les unités expérimentales sont assignées au traitement sans influence mutuelle directe.
Si ces conditions ne sont pas respectées, les tests d’ANCOVA (estimations, F, p-values) ne sont plus valides car l’hypothèse d’erreurs indépendantes est violée.
Il n'existe pas de test R pour faire ce contrôle, cela repose sur l'expérimentateur qui, par exemple si mesure l'effets de groupes en intégrant l'âge, sait que, par exemple, il n'a pas biaisé ses groupes en prenant des groupes qui présentent intrinsèquement des âges différents.
L’ANCOVA compare des groupes en ajustant l’effet d’une covariable : si les résidus (les écarts entre les données réelles et le modèle) ne sont pas normaux, les différences pourraient être mal ajustées d'un groupe à l'autre.
mes_residus <- residuals(resultat)
Méthode graphique : histogramme et QQ-plot
hist(mes_residus)
qqnorm(mes_residus) ; qqline(mes_residus)
De toute évidence, mon modèle ici n'est pas top, cela ne dessine ni une cloche ni un bel alignement en QQ-plot.
Méthode statistique : test de Shapiro ou autres tests selon le nombre de résidus.
shapiro.test(mes_residus)
La p-value = 0.007037 confirme que cela ne sent pas la normalité.
Remarque importante : Le test de Shapiro est surtout adapté aux petits échantillons. Pour des échantillons plus grands, il est préférable d’utiliser des tests plus souples (comme celui de Jarque-Bera) ou, au-delà d’une certaine taille, de se contenter de vérifier l’asymétrie (skewness, seuil absolu de 1) et l’aplatissement (kurtosis, seuil absolu de 1,5).
Pour chaque groupe de facteurs (oublions la covariable qui est continue) : il faut vérifier si les résidus ont aussi à peu près la même variance.
Contrairement à une ANOVA, on ne fait pas la variance des groupes, mais la variance des résidus des différents groupes.
Ce faisant, si les résidus sont normaux, on peut faire un test exigeant de Bartlett. Si en revanche, les résidus ne sont pas normaux, on peut passer à un test de Levene, voire de Brown-Forsyth.
bartlett.test(mes_residus, donnees$groupe)
Pour faire de même en Levene :
library(lawstat)
# Effectuer le test de Levene qui bascule automatiqueet sur Brown-Forsyth si nécessaire
levene.test(mes_residus, donnees$groupe)
Graphiquement, on peut aussi vérifier que les résidus par rapport aux valeurs ajustées sont globalement les mêmes :
res <- residuals(resultat)
fit <- fitted(resultat)
plot(fit, res,
xlab = "Valeurs ajustées",
ylab = "Résidus",
main = "Contrôle visuel de l'homogénéité des variances")
abline(h = 0)
Remarque très utile : si la variance se révèle homogène et que les résidus sont considérés comme non normaux, on peut retester cette normalité car l'ANCOVA comme l'ANOVA, tolère une violation légère de la normalité en situation d'homoscédasticité.
Il suffit de contrôler la normalité avec skweness et kurtosis et avec des seuils plus souples, ce que fait automatiquement m.test() de {KefiR}.
La relation entre covariable explicative et variable Y dépendante doit être linéaire.
Pour confirmer ce point, on ajoute un terme quadratique, la covariable au carré I(cov^2) à un modèle linéaire avec facteurs + covariable. Si l'ajout de ce terme améliore significativement le modèle, c'est que la relation n'est pas linéaire.
resultat <- aov(performance ~ age + groupe, # ATTENTION de mettre age (covariable avant groupe)
data = donnees)
resultat2 <- aov(performance ~ age + I(age^2) + groupe,
data = donnees)
anova(resultat, resultat2)
La fonction anova() montre ici que resultat2 est significativement meilleur que resultat.
Ici, la linéarité covariable - variable dépendante n'est pas respectée. On devrait passer à une alternative de l'ANCOVA...
Remarque : on notera ici que comme je n'ai pas fait du type 3, l'ordre des temes est très important : ne pas mettre le facteur avant la covariable.
Source : Huitema, B. E. (2011). The Analysis of Covariance and Alternatives: Statistical Methods for Experiments, Quasi-Experiments, and Single-Case Studies (2ᵉ éd.). John Wiley & Sons.
plot(donnees$age, donnees$performance, col=as.factor(donnees$groupe),pch=16,cex=2)
for (i in unique(donnees$groupe) ) {
abline(lm(performance~age, data=donnees[donnees$groupe==i,]), lwd=2)
}
A l'œil, si les données suivent les régressions linéaires de chaque groupe, alors c'est bon !
Dans une ANCOVA classique, les pentes doivent être les même quel que soit le groupe. En d'autres termes, il ne doit pas y avoir d'interaction entre le facteur et la covariable ce qu'il impliquerait un effet différent d'une modalité à l'autre du facteur.
Le contrôle est simple, il suffit de faire un modèle avec interaction et vérifier que cette interaction n'a pas d'effet significatif.
resultat <- aov(performance ~ age + groupe, # ATTENTION de mettre age (covariable avant groupe)
data = donnees)
resultat_interaction <- aov(performance ~ age * groupe,
data = donnees)
anova(resultat, resultat_interaction)
# Ici resultat_interaction n'est pas meilleur que resultat, donc pas de problème sur cette assomption.
Toutefois, qu'est ce que cela implique si je fais mon ANCOVA de façon non classique, avec une interaction ?
Tout simplement, comme l'ANCOVA ici a pour but d'éliminer l'effet la covariable (ex : effet de l'âge), s'il y a interaction, cela n'est plus possible. Les moyennes ajustées perdent leur sens, car l’ajustement de la covariable n’est pas le même pour chaque groupe.
De fait, on ne doit plus parler d'ANCOVA, on doit alors parler de modèle linéaire général à pentes hétérogènes (cf. 7. Les alternatives à l'ANCOVA classique, ci-après).
Le contrôle visuel de l’homogénéité des pentes utilise exactement le même graphique que celui mobilisé pour l’assomption 4. Linéarité de la relation covariable – variable dépendante.
Il s’agit de représenter la variable dépendante en fonction de la covariable, en colorant les points selon les groupes (ou les combinaisons de groupes lorsqu’il y a plusieurs facteurs), puis d’ajouter une droite de régression pour chaque groupe.
Les pentes doivent être visuellement parallèles pour que l’assomption soit considérée comme respectée.
Des divergences nettes, des pentes de signes différents, ou des croisements de lignes indiquent une violation de l’homogénéité des pentes et impliquent de passer à un modèle linéaire général avec interaction facteur × covariable.
Les outliers sont à regarder en réalité attentivement avant même le contrôle des assomptions car : des valeurs aberrantes, voire extrêmes fausses tout le modèle ANCOVA et conduisent à des rejets à tort du modèle ou, parfois inversement et plus rarement, peuvent conduirent à conserver un modèle ANCOVA qui devrait être rejeté. Les outliers :
Faussent la normalité des résidus
Faussent la variance des groupes
Induisent de mauvaises pentes entre les groupes.
Déforment les ajustement.
Inversent les conclusions.
La présence d'outliers doit donc conduire, ou à leur élimination manuelle ou à un changement de test statistique.
Avec identify_outliers() de {rstatix}, il est facile d'identifier les outliers par le code suivant :
by(donnees$performance, donnees$groupe,function(x){identify_outliers(data.frame(x))})
Plusieurs cas s'offre à nous :
Une augmentation accrue du nombre de facteurs et/ou de covariables
L'introduction d'interaction : modèle linéaire général à pentes hétérogènes
Le nom respect des assomptions ou la présence d'outliers/données aberrantes.
======================
Consigne ci-dessus :
En ANCOVA : écarter les données aberrantes avec rstatix::identify_outlier comme cela a été fait en ANOVA pour .multi_factor_analysis.
Contrôle de l'homogénéité de la variance est bien fait sur les résidus en ANCOVA : cela n'apparait pas dans les commentaires, on pourrait croire que cela est fait pour les groupes à partir de la variable dépendante.
Vérifier que en cas de non-normalité des résidus, et homogénéité des variances, on tests bien un retour à une situation paramétrique en recontrôlant la normalité en mode extrem - comme pour l'ANOVA avec .multi_factor_analysis. Sinon, l'implémenter.
Je m'interroge si m.test() en mode ANCOVA tolère les données appariées/répétées. Il faut voir si oui, et si non, qu'il dirige l'utilisateur par un message approprié. De même pour les effets aléatoires, prévoir un routage vers une action adaptée ou un message présentant les limites de m.test() à l'utilisateur.
De la même façon, il me semble que la MANOVA ignore tout de la MANCOVA. L'entrée d'une MANVOCA sur m.test() doit donc conduire à un message claire pour l'utilisateur. A vérifier et implémenter si nécessaire.
De façon gérénale, il faudrait ajouter à bp.log puis appliquer à l'ensemble de m.test() un message qui, à chaque fois que m.test() annonce une limite non traitée par elle-même et invitant l'utilisateur à se débrouiller en allant vers d'autres tests (lme4, glm...), l'informe que cette fonction n'est pas encore implémentée sous m.test(), mais qu'il ne faut pas hésiter à contacter l'auteur à antoine.masse@u-bordeaux.fr si on souhaite voir cette fonctionnalité implémenter dans m.test() [prévoir d'envoyer la ligne de code utilisé et le jeu de données concernée afin que.... à finir]
Concernant l'introduction d'interaction covariable * facteur en m.test() ==> Le traitement est fait par les fonctions d'ANCOVA, pourtant : ce n'est plus une ANCOVA mais un "modèle linéaire général à pentes hétérogènes" : cela doit être mentionné dans le verbose, mais aussi, changer l'approche car on devrait dès le début annoncer qu'on passe en modèle linéaire général à pentes hétérogènes" et non ANCOVA du fait de l'interaction ce qui dégage de fait le contrôle de cette assomption. On doit aussi vérifier que le post-hoc est adapté car il n'y a plus aucun sens de contrôler le facteur ou la covariable seul ET, de PLUS, le mode graphique doit s'adapter (ce n'est plus de l'ANCOVA). En résumé, Si une interaction facteur × covariable est présente dans le modèle, la condition d’homogénéité des pentes cesse d’exister. Le modèle autorise explicitement des pentes différentes. Les hypothèses restantes sont celles du modèle linéaire (normalité, indépendance, homoscédasticité, linéarité dans les groupes). L’interprétation ne peut plus porter sur les effets principaux, mais exclusivement sur la façon dont la relation covariable–réponse diffère selon les groupes (visualisée par les droites de régression par groupe).
Autre point à vérifier améliorer si justifié d'un point de vue académique concernant l'ANCOVA - et surtout le cas où le nombre de faceurs augmente et/ou le nombre de covariables explicatives :
Absence de test formel du modèle complet vs. modèle réduit.
Les interactions facteur×covariable ne sont jamais comparées globalement via anova(mod_reduced, mod_full).
→ À ajouter impérativement pour décider objectivement entre pentes communes et pentes spécifiques.
Absence totale d’évaluation de la multicolinéarité.
Aucun calcul de VIF ou de matrice de corrélation entre covariables.
→ Intégrer un bloc car::vif() après ajustement du modèle.
Absence d’évaluation de l’indépendance entre covariables et facteurs.
Aucune ANOVA du type covariate ~ factor.
→ Ajouter une boucle systématique testant chaque covariable contre chaque facteur.
=========================
Ainsi tu dois, reprendre l'ensemble de mes consignes et interrogations ci-dessus pour améliorer la consigne ci-dessous avant de l'exécuter :
==========================
Consigne ci-dessous :
Auditer l’ensemble du pipeline ANCOVA de m.test()
(.ancova_analysis(), .posthoc_ANCOVA(), .variance(), .normality())
en appliquant strictement les standards statistiques académiques.
Ajouter :
un module de sélection automatique du modèle (interaction ou non) basé sur un test hiérarchique mod_reduced vs mod_full ;
un module de détection de multicolinéarité (VIF > 5 ou > 10 selon standard retenu) ;
un module de détection de dépendance covariable–facteur (tests univariés systématiques).
Documenter clairement chaque décision (pentes homogènes ou non, colinéarité détectée, dépendance détectée) dans la structure assumptions_checked.
Conserver et respecter la structure actuelle de m.test(), ses objets, ses messages bilingues, et son fonctionnement en mode robuste.
Liens externes :
Un site intéressant pour se donner un autre regard sur le sujet : le site de Datanovia.