© 2025 Tous droits réservés
Il m’arrive très régulièrement d’animer des formations de remise à niveau en biostatistiques auprès de chercheurs, de techniciens, ou de chargés de projet, qui évoluent dans les domaines de la biologie, de l’agronomie, ou encore de la nutrition animale. Et lorsqu’on aborde les ANOVA à un facteur (analyse de variance à 1 facteur), neuf fois sur dix, au moins un stagiaire me demande comment identifier les groupes homogènes après les tests post-hoc, et comment les afficher sur un graphique en utilisant des lettres.
J’avais déjà écrit des articles sur le sujet, mais cela fait maintenant plusieurs années, et depuis, j’ai amélioré mes procédures (j’emploie notamment le package emmeans
pour les comparaisons multiples).
Alors, j’ai préféré réécrire un tutoriel complet pour vous montrer, en pas à pas, comment réaliser une ANOVA à un facteur (paramétrique et non paramétrique) avec R, comment vérifier les conditions d’application, comment réaliser les tests post-hoc si l’ANOVA est significative, et finalement, comment synthétiser les résultats en identifiant les groupes homogènes et comment les représenter sur un graphique, à l’aide de lettres.
Nous allons employer les données chickwts
qui sont présentes dans le package dataset
.
Dans ce jeu de données, des poussins nouvellement éclos ont été répartis aléatoirement en six groupes, et chaque groupe a reçu un complément alimentaire différent. Leurs poids en grammes après six semaines sont donnés ainsi que les types de compléments alimentaires.
data("chickwts")
head(chickwts)
weight feed
1 179 horsebean
2 160 horsebean
3 136 horsebean
4 227 horsebean
5 217 horsebean
6 168 horsebean
Nous pouvons commencer par visualiser les données :
library(tidyverse)
ggplot(chickwts, aes(x=feed, y=weight,colour=feed ))+
geom_jitter(height=0, width=0.15)+
geom_boxplot(alpha=0.25, outlier.shape = NA)+
theme_minimal()+
theme(legend.position="none")
Cette étape n’est pas indispensable, mais comme le type de complément alimentaire n’est pas une variable ordinale, on peut se permettre de modifier l’ordre des modalités pour obtenir un graphique ordonné :
chickwts$feed <- chickwts$feed %>%
fct_relevel("horsebean",
"linseed",
"soybean",
"meatmeal",
"sunflower",
"casein")
ggplot(chickwts, aes(x=feed, y=weight,colour=feed ))+
geom_jitter(height=0, width=0.15)+
geom_boxplot(alpha=0.25, outlier.shape = NA)+
theme_minimal()+
theme(legend.position="none")
La première étape, consiste à réaliser l’ANOVA. Ici, j’emploie la fonction aov()
:
chickwts.aov1 <- aov(weight~feed, data=chickwts)
Visualisation des résultats :
Une fois l’ANOVA réalisée, j’aime bien afficher les résultats pour vérifier que l’ajustement s’est bien déroulé. Il s’agit seulement de vérifier que toutes les quantités ont été calculées, et que les degrés de libertés correspondent bien à ce que nous pensons.
summary(chickwts.aov1)
Df Sum Sq Mean Sq F value Pr(>F)
feed 5 231129 46226 15.37 5.94e-10 ***
Residuals 65 195556 3009
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Le nombre de degrés de liberté pour la variable feed
correspond au nombre de modalités (complément alimentaire) – 1. Il y a 6 compléments différents donc 5 degrés de liberté. Il n’y a donc pas d’erreur. Ce qui peut arriver c’est qu’une modalité soit enregistrée avec deux orthographes différentes et crée artificiellement un groupe supplémentaire.
Le nombre de degrés de liberté pour les résidus correspond au nombre total de données – le nombre de groupe ; ici 71-6 = 65. Là encore, c’est conforme.
À ce stade, je n’interprète pas les résultats, car pour avoir confiance dans les résultats, les conditions d’application de l’ANOVA doivent être satisfaites.
L’ANOVA à un facteur a 3 conditions d’application :
Indépendance des données :
L’hypothèse d’indépendance des données se juge sur le design expérimental. Ici, les poussins ont été répartis aléatoirement en six groupes, et chaque groupe a reçu un complément alimentaire différent ; les données sont donc indépendantes. Lorsque les données ne sont pas indépendantes, les unités expérimentales reçoivent successivement les différents traitements.
Normalité des résidus :
Cette hypothèse peut être évaluée visuellement par un qqplot :
library(car)
qqPlot(chickwts.aov1)
Ici les points sont bien alignés le long de la droite et sont tous inclus dans l’enveloppe de confiance. Ainsi, rien ne permet de penser qu’il existe un défaut de normalité.
L’hypothèse de normalité peut également être évaluée par un test statistique ; le plus courant étant le test de Shapiro-Wilk .
shapiro.test(residuals(chickwts.aov1))
Shapiro-Wilk normality test
data: residuals(chickwts.aov1)
W = 0.98616, p-value = 0.6272
Pour rappel, l’hypothèse nulle du test de Shapiro-Wilk suppose la normalité des données, alors que l’hypothèse alternative suppose un défaut de normalité. La pvalue étant > 0.05, là encore rien ne permet d’affirmer qu’il existe un défaut de normalité. Nous pouvons finalement considérer que l’hypothèse de normalité des résidus est satisfaite.
Égalité des variances
Cette hypothèse peut être explorée en comparant la hauteur des boites des boxplots. Sur le graphique réalisé précédemment, on peut voir que la variabilité des données est plus faible pour le régime “sunflower” mais il est difficile de savoir si cela est “suffisamment” important pour rejeter l’hypothèse d’égalité des variances.
Cette hypothèse d’égalité des variances peut être évaluée par un test statistique, classiquement le test de Bartlett ou le test de Levene. Pour rappel, l’hypothèse nulle de ces tests suppose l’égalité des variances, alors que l’hypothèse alternative suppose un défaut d’égalité des variances.
Voici plusieurs façons de les réaliser :
library(performance)
performance :: check_homogeneity(chickwts.aov1)
OK: There is not clear evidence for different variances across groups (Bartlett Test, p = 0.660).
stats::bartlett.test(chickwts$weight~chickwts$feed)
Bartlett test of homogeneity of variances
data: chickwts$weight by chickwts$feed
Bartlett's K-squared = 3.2597, df = 5, p-value = 0.66
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 0.7493 0.5896
65
Les pvalues étant > 0.05, rien ne nous permet de conclure qu’il existe un défaut d’égalité des variances et nous pouvons finalement considérer que l’hypothèse d’égalité des variances est satisfaite.
Les conditions d’application étant satisfaites, nous pouvons avoir confiance dans les résultats de l’ANOVA, et donc les interpréter.
Pour rappel, l’hypothèse nulle de l’ANOVA suppose que les moyennes des différents groupes sont égales, alors que l’hypothèse alternative suppose qu’elles sont globalement différentes (ce qui veut dire qu’il existe au moins un couple de moyennes différentes).
summary(chickwts.aov1)
Df Sum Sq Mean Sq F value Pr(>F)
feed 5 231129 46226 15.37 5.94e-10 ***
Residuals 65 195556 3009
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ici, la pvalue associée à la variable feed
est <0.05. Nous pouvons donc rejeter l’hypothèse nulle et conclure que les moyennes des poids des poussins sont globalement différentes. Nous pouvons aussi formuler le rejet de l’hypothèse nulle comme “il existe un lien entre le type de complément alimentaire et le poids des poussins”
Remarque : Si les conditions n’étaient pas satisfaites, il faudrait réaliser une ANOVA non paramétrique (nommé également test de Kruskal Wallis (voir plus bas))
Quand l’effet de la variable étudié est significatif (pval <0.05) comme ici, on cherche ensuite à explorer où se situent les différences. Pour cela, on réalise des tests post-hoc. Il s’agit de comparaisons multiples qui consistent à faire des tests de Student (lorsque l’ANOVA paramétrique a été employée) deux à deux.
Les comparaisons multiples entraînent une augmentation du risque alpha global, c’est-à-dire le risque global de se tromper au moins une fois en concluant à la différence entre deux moyennes si celle-ci n’existe pas en réalité. Pour ne pas voir le risque alpha augmenter, les p_values des comparaisons multiples doivent être ajustées (ou parle aussi de correction). Il existe différentes procédures : Tukey, Holm Bonferoni, etc… L’approche qui me semble la plus employée est celle de “holm ». C’est celle que je vais employer ici.
Pour réaliser les tests post hoc, nous pouvons employer la fonction emmeans
du package emmeans
, comme cela :
library(emmeans)
emmeans(chickwts.aov1, specs=pairwise~feed, adjust="holm")$contrasts
contrast estimate SE df t.ratio p.value
horsebean - linseed -58.55 23.5 65 -2.493 0.0944
horsebean - soybean -86.23 22.7 65 -3.797 0.0030
horsebean - meatmeal -116.71 24.0 65 -4.870 0.0001
horsebean - sunflower -168.72 23.5 65 -7.184 <.0001
horsebean - casein -163.38 23.5 65 -6.957 <.0001
linseed - soybean -27.68 21.6 65 -1.283 0.5177
linseed - meatmeal -58.16 22.9 65 -2.540 0.0944
linseed - sunflower -110.17 22.4 65 -4.920 0.0001
linseed - casein -104.83 22.4 65 -4.682 0.0002
soybean - meatmeal -30.48 22.1 65 -1.379 0.5177
soybean - sunflower -82.49 21.6 65 -3.823 0.0030
soybean - casein -77.15 21.6 65 -3.576 0.0053
meatmeal - sunflower -52.01 22.9 65 -2.271 0.1322
meatmeal - casein -46.67 22.9 65 -2.039 0.1823
sunflower - casein 5.33 22.4 65 0.238 0.8125
P value adjustment: holm method for 15 tests
Pour évaluer la significativité statistiques des différences, nous devons comparer les pvalues qui ont été ajustées, à la valeur 0,05 du risque alpha. Ici, les différences significatives mises en évidence sont :
Il est possible de visualiser les différences comme ceci :
mc <- emmeans(chickwts.aov1, specs=pairwise~feed, adjust="holm")
plot(mc, comparisons=TRUE)
Les barres bleues indiquent les intervalles de confiances à 95% des moyennes et les flèches rouges indiquent les différences significatives entre les moyennes. D’après ma compréhension, les flèches rouges dépassent les zones bleues, car les pvalues ont été ajustées pour tenir compte des comparaisons multiples. Si deux flèches se superposent alors cela correspond à une pvalue ajustée > 0,05. Inversement, lorsque 2 flèches ne se superposent pas, alors cela correspond à une pvalue ajustée < 0,05.
Maintenant que nous savons où se situent les différences significatives, nous souhaitons synthétiser ces résultats en identifiant des groupes homogènes. Pour cela, nous pouvons employer la fonction cld
du package multcomp
. Cette fonction permet d’identifier les sous groupes homogènes à partir des résultats des tests post-hoc de l’ANOVA, en utilisant des lettres pour indiquer les différences significatives entre les moyennes ou les traitements.
Lorsque deux groupes partagent la même lettre, cela signifie que rien ne nous permet d’affirmer qu’ils sont significativement différents.
En revanche, si deux groupes ne partagent pas la même lettre, alors cela signifie qu’ils sont significativement différents au risque de 5%.
library(multcomp)
LetterResults <-
emmeans(chickwts.aov1, specs = pairwise ~ feed, adjust = "holm") %>%
multcomp::cld(Letters = letters)
LetterResults
feed emmean SE df lower.CL upper.CL .group
horsebean 160 17.3 65 126 195 a
linseed 219 15.8 65 187 250 ab
soybean 246 14.7 65 217 276 b
meatmeal 277 16.5 65 244 310 bc
casein 324 15.8 65 292 355 c
sunflower 329 15.8 65 297 361 c
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 6 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Pour présenter les résultats, j’aime bien employer la fonction ggbetweenstats
du package ggstatsplot
. Et il est très facile de rajouter l’information sur les sous-groupes homogènes que nous venons de définir, en utilisant une couche geom_text()
. L’argument y dans la fonction aes()
permet de spécifier la position de la lettre sur l’axe des y ; ici à la valeur 440. `
library(ggstatsplot)
ggbetweenstats(data=chickwts,
x=feed,
y=weight,
type="parametric",
var.equal = TRUE,
p.adjust.method = "holm",
title = "Efficacité de 6 suppléments alimentaires sur le taux de croissance des poulets")+
geom_text(data=LetterResults, aes(label=.group, y=440))
Si vous souhaitez supprimer les pvalues, je n’ai rien trouver de mieux que de :
results.subtitle = FALSE)
dans la fonction ggbetweenstats()
library(ggstatsplot)
ggbetweenstats(
data = chickwts,
x = feed,
y = weight,
type = "parametric",
var.equal = TRUE,
p.adjust.method = "holm",
title = "Efficacité de 6 suppléments alimentaires sur le taux de croissance des poulets",
results.subtitle = FALSE) + # Masque le résultat global de l'ANOVA et les p-values)
scale_y_continuous(
limits = c(100, 445),
breaks = seq(from = 100, to = 444, by = 100)
)+
geom_text(data=LetterResults, aes(label=.group, y=440))
La même démarche peut être réalisé en cas d’utilisation d’une ANOVA non paramétrique?
L’ANOVA non paramétrique, qui s’appelle aussi le test de Kruskal- Wallis peut être réalisée à l’aide de la fonction kruskal.wallis() :
chickwts.kw <- kruskal.test(weight~feed, data=chickwts)
chickwts.kw
Kruskal-Wallis rank sum test
data: weight by feed
Kruskal-Wallis chi-squared = 37.343, df = 5, p-value = 5.113e-07
L’ANOVA non paramétrique étant significative (pval <0.05), l’étape suivante consiste à réaliser les tests post hoc, puisqu’aucune condition d’application n’est à évaluer (à part l’indépendance des données qui est vérifiée par le design expérimental).
Lorsque l’approche est non paramétrique, les comparaisons multiples peuvent être réalisées à l’aide du test de Dunn. Les résultats sont fournis sous la forme d’une matrice de p-values :
library(PMCMRplus)
dunn_results<-PMCMRplus::kwAllPairsDunnTest(weight~feed, data=chickwts,p.adjust.method = "holm")
Warning in kwAllPairsDunnTest.default(c(179, 160, 136, 227, 217, 168, 108, :
Ties are present. z-quantiles were corrected for ties.
dunn_results
Pairwise comparisons using Dunn's all-pairs test
data: weight by feed
horsebean linseed soybean meatmeal sunflower
linseed 0.5830 - - - -
soybean 0.0834 0.9900 - - -
meatmeal 0.0092 0.4815 0.9900 - -
sunflower 9.2e-06 0.0062 0.0715 0.5830 -
casein 2.1e-05 0.0103 0.0994 0.6274 0.9900
P value adjustment method: holm
alternative hypothesis: two.sided
Ici, les différences significatives mises en évidence sont :
Certaines différences mises en évidence par l’ANOVA paramétrique ne sont pas retrouvées par l’ANOVA non paramétrique.
Pour synthétiser les résultats et identifier les groupes homogènes, nous pouvons employer la fonction multcompLetters()
du package multcompView
:
# Créer un tableau avec les résultats et obtenir les sous-groupes
library(multcompView)
p_values <- dunn_results$p.value
group_letters <- multcompLetters(p_values, compare = "<", threshold = 0.05)
group_letters
linseed soybean meatmeal sunflower casein
"a" "ab" "bc" "c" "c"
# Extraire les lettres associées aux groupes
letters_df <- data.frame(feed = names(group_letters$Letters),
Letters = group_letters$Letters)
letters_df
feed Letters
linseed linseed a
soybean soybean ab
meatmeal meatmeal bc
sunflower sunflower c
casein casein c
# Créer le plot et ajouter les lettres sur le graphique
ggbetweenstats(data = chickwts,
x = feed,
y = weight,
type = "nonparametric",
var.equal = TRUE,
p.adjust.method = "holm",
title = "Efficacité de 6 suppléments alimentaires sur le taux de croissance des poulets") +
geom_text(data = letters_df, aes(x = feed, y = 440, label = Letters))
Et pour supprimer les pvalues du graphique :
ggbetweenstats(data = chickwts,
x = feed,
y = weight,
type = "nonparametric",
var.equal = TRUE,
p.adjust.method = "holm",
title = "Efficacité de 6 suppléments alimentaires sur le taux de croissance des poulets",
results.subtitle = FALSE)+
scale_y_continuous(
limits = c(100, 445),
breaks = seq(from = 100, to = 444, by = 100)
)+
geom_text(data = letters_df, aes(x = feed, y = 440, label = Letters))
Dans ce tutoriel, nous avons vu comment réaliser des ANOVA paramétrique et non paramétrique à un facteur avec R, comment vérifier les conditions de validité, puis comment réaliser les tests post-hoc pour identifier où se situent les différences significatives entre les moyennes (si le test global est significatif), et enfin, comment synthétiser les résultats de ces comparaisons en identifiant les groupes homogènes et en les affichant sur un graphique, en employant des lettres.
J’espère que ces explications vous aideront à résoudre un problème concret dans vos analyses de données.
Si vous avez des questions, des difficultés ou des retours d’expérience, N’hésitez pas à laisser un commentaire ci-dessous pour partager vos réflexions !
Retrouvez le programme de ma formation d’initiation à Quarto ainsi que le planning des prochaines session en cliquant ici
C’est possible en faisant un don sur la page Tipeee du blog
© 2025 Tous droits réservés
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.
12 réponses
Bonjour,
Vraiment très très intéressant
Merci grandement !
C est génial vraiment merci beaucoup
Merci pour cet article enrichissant, mais je racontre un problème avec cette partie du script, pourriez m’aider ?
« `{r}
library(multcomp)
LetterResults %
multcomp::cld(Letters = letters)
LetterResults
« `
Erreur dans UseMethod(« cld ») :
pas de méthode pour ‘cld’ applicable pour un objet de classe « emmGrid »
Bonjour,
Est-ce que vous obtenez cette erreur en reproduisant l’exemple ou sur d’autres données ?
Merci pour le billet, il est très clair dans l’ensemble.
Quelques remarques/questions toutefois :
– Il me semble que les tests de vérifications des hypothèses doivent être réalisés sur les résidus du modèle de régression et non sur les valeurs brutes. Le Shapiro est réalisé sur les résidus mais pas les Bartlett et Levene. J’ai déjà vu un modèle dont les résidus ne répondaient pas au critère d’homogénéité alors qu’un Levene sur les variances des données brutes était OK.
– Je ne connaissais pas le test d’homoscédasticité sur les résidus du modèle, je faisais toujours une inspection visuelle mais c’est vrai que lorsque j’avais un doute, j’étais bien en peine pour tester avec une procédure formelle (i.e. un indicateur statistique et une p-val).
– A quoi correspondent les contrastes dans les modèles linéaires ?
– J’utilise classiquement Tukey, Bonferroni ou Holm mais sans jamais trop savoir lequel utiliser et pourquoi. Un petit billet sur ces sujets là pour compléter celui-ci ?
Bonjour,
Pour la normalité, oui, c’est sur les résidus.
Par contre, pour l’homogénéité, d’après ma compréhension, cela revient au même de le faire sur les résidus ou sur les données brutes (dans le cadre de l’ANOVA à un facteur)). J’ai regardé dans 2 livres (Stats facile avec R et Introductory Statistics with R) et la condition est évaluée sur les données brutes. Si vous avez l’exemple sous la maain dont vous parlez, je veux bien y jeter un coup d’œil.
Les contrastes sont des différences en quelque sorte.
Concernant les approches pour les comparaisons multiples, j’avais écrit un article dessus (il ya trèèès longtemps) : https://delladata.fr/comparaisons-multiples-et-ajustement-des-pvalues-avec-le-logiciel-r/
Bonne continuation
Merci Claire pour cette révision de l’Anova à un facteur. Moi j’avais un souci sur les tests post-hoc. Il m’arrive parfois d’utiliser différents tests pour la comparaison multiples afin être sûr des résultats obtenus (SNK test, tukey et autres). Parfois je trouve des résultats différents. La question est de savoir quel test le plus robuste pour ces analyses? Et t’il aussi possible de réaliser volontairement le test student paramétrique pour la comparaison multiples ?
Bonjour,
je ne sais pas si une approche est plus robuste que les autres dans toutes les situations. Je ne pense pas.
Je vous conseille d’en choisir une et de vous y tenir (vous pouvez employer celle qui est le plus généralement employée dans votre domaine d’étude).
Je ne suis pas sûre de comprendre votre question sur les tests de stuent pour les comparaisons multiples. Vous pouvez regarder du coté de la fonction pairwise.t.test() qui permet d’utiliser différentes approches d’ajustement.
Bonne continuation
Bonjour, merci beaucoup pour ce tutoriel, ça aide beaucoup!
Je me demandais cependant, dans l’ANOVA non-paramétrique, pourquoi les lettres ne sont pas attribuées à tous les groupes? Dans votre cas, « horsebean » n’a pas de lettre et n’est pas comparé aux autres groupes. J’avais une analyse à peu près similaire à faire et la même chose s’est produite: mon premier groupe n’est pas considéré. Comment corriger cela?
Merci beaucoup!
Bonjour,
Merci pour votre question, elle est très pertinente.
Sur les résultats des pvalues, on voit que horsbean est significativement différent des autres groupes.
Je ne sais pas pourquoi la fonction multcompLetters() ne le prend pas en compte.
Voici un code qui devrait vous aider :
# Réaliser le test de Dunn
dunn_results <- PMCMRplus::kwAllPairsDunnTest(weight ~ feed, data = chickwts, p.adjust.method = "holm")# Obtenir les p-values p_values <- dunn_results$p.value# Vérifier que tous les groupes sont bien présents dans le jeu de données original all_groups <- levels(chickwts$feed)# Appliquer multcompLetters pour obtenir les lettres des groupes group_letters <- multcompLetters(p_values, compare = "<", threshold = 0.05)# Créer un tableau de lettres pour les groupes présents dans les p-values letters_df <- data.frame(feed = names(group_letters$Letters), Letters = group_letters$Letters)# Ajouter les groupes manquants (s'ils existent) avec une lettre distincte (par exemple "d") missing_groups <- setdiff(all_groups, letters_df$feed)if(length(missing_groups) > 0) {
# Ajouter les groupes manquants avec une lettre distincte (ou ajustée selon vos besoins)
for(group in missing_groups) {
letters_df <- rbind(letters_df, data.frame(feed = group, Letters = "d")) } }# Résultat final letters_df #plot ggbetweenstats(data = chickwts, x = feed, y = weight, type = "nonparametric", var.equal = TRUE, p.adjust.method = "holm", title = "Efficacité de 6 suppléments alimentaires sur le taux de croissance des poulets", results.subtitle = FALSE)+ scale_y_continuous( limits = c(100, 445), breaks = seq(from = 100, to = 444, by = 100) )+ geom_text(data = letters_df, aes(x = feed, y = 440, label = Letters))
Bonne continuation
Astuce transmise par Antoine Soetewey :
Si vous désirez supprimer les p-valeurs et accolades des tests post-hoc tout en gardant les résultats de l’ANOVA/Kruskal-Wallis dans le sous-titre du graphique, il existe l’argument pairwise.display = « none » (dans la fonction ggbetweenstats et ggwithinstats).
Merci beaucoup madame !
Pouvez vous nous faire un article sur les modèles linéaires mixtes ?