Comment afficher les p-values sur un graph ?
Lorsqu’on souhaite synthétiser les résultats d’une comparaison de moyennes, on réalise généralement une représentation graphique sur laquelle on reporte les comparaisons réalisées, et les résultats des tests statistiques employés.
Les résultats peuvent être exprimés en termes de p-values, par exemple, comme ceci :

Ou encore, en employant des lettres, comme ci dessous. La convention veut que si deux moyennes partagent une même lettre alors elles ne sont pas significativement différentes, et au contraire, si deux moyennes ne partagent pas une même lettre alors elles sont significativement différentes

La réalisation de ces deux types de figures, comportant une annotation des résultats des comparaisons n’est pas très aisée avec le logiciel R.
Néanmoins, deux outils permettent de les construire plus facilement. Il s’agit d’une part du package ggpubr développé par Alboukadel Kassambara, qui permet d’indiquer les comparaisons réalisées, et les p-values des tests employés. Et d’autre part, d’une fonction (nomméetri.to.squ
) publiée sur le blog de Fabio Marroni qui, couplée à la fonction multcompLetters
du package multcompView
, permet facilement d’employer des lettres pour indiquer les résultats des comparaisons.
Dans ce post, je vais donc vous montrer, en pas à pas, comment construire les deux figures précédentes.
Afficher les résultats des comparaisons à l'aide des p-values
Représenter les données
Le package ggpubr
propose des fonctions “toutes prêtes” pour représenter les données à comparer. Ces fonctions sont basées sur une utilisation indirecte du package ggplot2
. On parle de fonctions “wrapper” de ggplot2
, puisqu’elles simplifient son utilisation.
A mon sens les fonctions les plus utiles pour illustrer une comparaison de moyennes sont ggboxplot
et ggviolin
.
library(ggpubr)
data("ToothGrowth")
df <- ToothGrowth
ggviolin(df, x = "dose", y = "len", fill = "dose", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
add = "boxplot", add.params = list(fill = "white"))+
# permet de ne pas afficher la légende
theme(legend.position="none")

Ou encore, par exemple :
ggboxplot(df, x = "dose", y = "len",
color = "dose",
palette =c("#00AFBB", "#E7B800", "#FC4E07"),
add = "jitter",
shape= "dose") +
theme(legend.position="none")

Ajouter les p-values
Les comparaisons statistiques des moyennes, et l’affichage des p-values correspondantes, sont réalisés par la fonction stat_compare_means
.
Lorsque plus de deux moyennes sont à comparer, la fonction permet de réaliser deux niveaux de test : un test de comparaison gobal, et des comparaisons multiples deux à deux.
La p-value du test global
Le premier niveau est le test global. Il permet d’évaluer si au moins deux moyennes sont différentes. Il s’agit d’une analyse de variance (ANOVA). En pratique, la fonction stat_compare_means
permet de réaliser:
- une ANOVA paramétrique, à l’aide de l’argument
method="aov"
oumethod="anova"
- une ANOVA non paramétrique, à l’aide de l’argument
method="kruskal.test"
. Il s’agit de l’option par défaut.
La fonctionstat_compare_means
s’ajoute simplement à la fonction permettant de faire le plot, comme ceci :
ggviolin(df, x = "dose", y = "len",
fill = "dose",
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
add = "boxplot",
add.params = list(fill = "white")) +
# ajout de la p-value du test global
stat_compare_means(method = "kruskal.test", label.y=45)+
theme(legend.position="none")

Remarque : l’argument label.y =
permet de spécifier la position de la p-value, sur l’axe de y.
Ci dessous avec une ANOVA paramétrique :
ggviolin(df, x = "dose", y = "len",
fill = "dose",
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
add = "boxplot", add.params = list(fill = "white")) +
# ajout de la p-value du test global
stat_compare_means(method = "anova", label.y=30)+
theme(legend.position="none")

Les p-values des comparaisons deux à deux
Le deuxième niveau de test correspond aux comparaisons deux à deux entre les moyennes. Pour cette étape, la fonction stat_compare_means
permet de réaliser :
- des tests t de Student (qui sont des tests paramétriques) par l’intermédiaire de l’argument :
method = "t.test"
. - des tests non paramétriques de Wilcoxon, en utilisant l’argument :
method = "wilcox.test"
.
Pour afficher les pvalues des comparaisons deux à deux, celles-ci doivent être spécifiées en amont, puis passées à la fonction stat_compare_means
par l’intermédiaire de l’argument comparisons
, comme ceci :
my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )
ggviolin(df, x = "dose", y = "len",
fill = "dose",
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
add = "boxplot",
add.params = list(fill = "white")) +
# ajout de la p-value globale
stat_compare_means(method = "kruskal.test", label.y=45)+
# ajout des p-values des comparaisons deux à deux
stat_compare_means(comparisons = my_comparisons, method="wilcox.test",label.y = c(29, 35, 40))+
theme(legend.position="none")

Lorsque seulement deux moyennes sont comparées, un seul niveau de p-value est évidemment affiché :
p <- ggboxplot(ToothGrowth, x = "supp", y = "len", color = "supp", palette = "jco",
add = "jitter") +
theme(legend.position="none")
# Add p-value
p + stat_compare_means()

Afficher les seuils de significativité :
Plutôt que les p-values, il est possible d’afficher dees seuils de significativité, sous la forme de nombre d’étoiles, par l’intermédiaire de l’argument label = "p.signif"
de la fonction stat_compare_means
.
ggboxplot(df, x = "dose", y = "len",
color = "dose",
palette =c("#00AFBB", "#E7B800", "#FC4E07"),
add = "jitter",
shape = "dose")+
## add global test p-value
stat_compare_means(method = "kruskal.test", label.y=45)+
## add pairwise comparison p values
stat_compare_means(comparisons = my_comparisons, method="wilcox.test",label.y = c(30, 35, 40), label="p.signif")+
theme(legend.position="none")

Les limites du package
Lorsqu’on réalise des comparaisons multiples de moyennes, deux à deux, il est usuel d’ajuster les p-values pour conserver un risque alpha global de 5%.
Alors que le package `ggpubr` propose une fonction, pour réaliser des comparaisons de moyennes, avec des p-values ajustées (fonction `compare_means()`), il ne semble pas possible d’afficher ces p-values ajustées sur les plots.
Mise à jour : Vous trouverez un article décrivant comment afficher les pvalues ajustées ici :
En pratique, la fonction `compare_means()`permet de réaliser des ANOVA paramétriques et non paramétriques (test de Kruskal Wallis ), ainsi que les comparaisons multiples avec différentes méthodes d’ajustement des p-values:


D’autres méthodes d’ajustement des pvalues, sont utilisables. Pour plus de détails, utilisez ?p.adjust.methods
.
Pour plus d’infos sur le package wilcox.t.test
, vous pouvez consulter [cette page, ou encore celle-ci.
Afficher les résultats des différences à l'aide de lettres
Principe
Pour réaliser ce type de représentation graphique, le principe est de :
- Faire les comparaisons de moyennes deux à deux avec les fonctions
pairwise.t.test
ouwilcox.t.test
. Et en ajustant les pvalues si on le souhaite à l’aide des méthodes disponibles (voir les méthodes disponibles de la fonctionp.adjust.methods
). - Récupérer la matrice triangulaire de p-values issues de ces fonctions, et la transformer en matrice rectangulaire par la fonction
tri.to.squ
. Le code de cette fonction est donnée plus bas. - Entrer la matrice rectangulaire de p-values en argument de la fonction
multcompLetters
du packagemultcompView
afin d’obtenir les lettres correspondant à chaque groupe (ou moyenne) comparé. - Conserver les lettres dans un data.frame.
- Faire le plot principal avec le package
ggplot2
- Ajouter les lettres sur ce plot avec la fonction
geom_text
La fonction tri.to.squ
Le code de la fonctiontri.to.squ
provient du blog de Fabio Marroni.
tri.to.squ<-function(x)
{
rn <- row.names(x)
cn <- colnames(x)
an <- unique(c(cn,rn))
myval <- x[!is.na(x)]
mymat <- matrix(1,nrow=length(an),ncol=length(an),dimnames=list(an,an))
for(ext in 1:length(cn))
{
for(int in 1:length(rn))
{
if(is.na(x[row.names(x)==rn[int],colnames(x)==cn[ext]])) next
mymat[row.names(mymat)==rn[int],colnames(mymat)==cn[ext]]<-x[row.names(x)==rn[int],colnames(x)==cn[ext]]
mymat[row.names(mymat)==cn[ext],colnames(mymat)==rn[int]]<-x[row.names(x)==rn[int],colnames(x)==cn[ext]]
}
}
return(mymat)
}
Exemple en pas à pas
# chargement des données
myeloma <- read.delim("https://raw.githubusercontent.com/kassambara/data/master/myeloma.txt")
# réalisation du test global par une anova non paramétrique (test de Kruskal Walis)
kruskal.test(DEPDC1 ~group, data=myeloma)
##
## Kruskal-Wallis rank sum test
##
## data: DEPDC1 by group
## Kruskal-Wallis chi-squared = 57.611, df = 6, p-value = 1.374e-10
En principe, les comparaisons multiples ne peuvent être réalisées que si le test global est significatif.
# réalisation des comparaisons deux à deux des moyennes par des tests de Wilcoxon, avec ajustement des pvalues par la méthode de Holm
# et récupération de la matrice triangulaire des p-values
pp <- pairwise.wilcox.test(myeloma$DEPDC1, myeloma$group,p.adjust.method ="holm" )
pp
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: myeloma$DEPDC1 and myeloma$group
##
## cyclin.d.1 cyclin.d.2 hyperdiploid low.bone.disease
## cyclin.d.2 1.0000 - - -
## hyperdiploid 0.0534 0.0351 - -
## low.bone.disease 0.0650 0.0650 1.0000 -
## maf 1.0000 1.0000 0.1306 0.1503
## mmset 1.0000 1.0000 0.0271 0.0650
## proliferation 0.0650 1.7e-05 1.9e-09 5.5e-09
## maf mmset
## cyclin.d.2 - -
## hyperdiploid - -
## low.bone.disease - -
## maf - -
## mmset 1.0000 -
## proliferation 0.1306 0.0014
##
## P value adjustment method: holm
# transformation de la matrice triangulaire de p-values en une matrice rectangulaire
mymat <-tri.to.squ(pp$p.value)
mymat
## cyclin.d.1 cyclin.d.2 hyperdiploid low.bone.disease
## cyclin.d.1 1.00000000 1.000000e+00 5.335428e-02 6.503880e-02
## cyclin.d.2 1.00000000 1.000000e+00 3.508640e-02 6.503880e-02
## hyperdiploid 0.05335428 3.508640e-02 1.000000e+00 1.000000e+00
## low.bone.disease 0.06503880 6.503880e-02 1.000000e+00 1.000000e+00
## maf 1.00000000 1.000000e+00 1.305730e-01 1.503143e-01
## mmset 1.00000000 1.000000e+00 2.705327e-02 6.503880e-02
## proliferation 0.06503880 1.699148e-05 1.859846e-09 5.518264e-09
## maf mmset proliferation
## cyclin.d.1 1.0000000 1.00000000 6.503880e-02
## cyclin.d.2 1.0000000 1.00000000 1.699148e-05
## hyperdiploid 0.1305730 0.02705327 1.859846e-09
## low.bone.disease 0.1503143 0.06503880 5.518264e-09
## maf 1.0000000 1.00000000 1.305730e-01
## mmset 1.0000000 1.00000000 1.399310e-03
## proliferation 0.1305730 0.00139931 1.000000e+00
library(multcompView)
# création des lettres correspondant à chaque moyenne
myletters <- multcompLetters(mymat,compare="<=",threshold=0.05,Letters=letters)
myletters
## cyclin.d.1 cyclin.d.2 hyperdiploid low.bone.disease
## "abc" "a" "b" "ab"
## maf mmset proliferation
## "abc" "a" "c"
# conserver les lettre dans un data.frame
myletters_df <- data.frame(group=names(myletters$Letters),letter = myletters$Letters )
myletters_df
## group letter
## cyclin.d.1 cyclin.d.1 abc
## cyclin.d.2 cyclin.d.2 a
## hyperdiploid hyperdiploid b
## low.bone.disease low.bone.disease ab
## maf maf abc
## mmset mmset a
## proliferation proliferation c
# plot
ggplot(myeloma, aes(x=group, y=DEPDC1, colour=group, fill=group))+
geom_boxplot(outlier.alpha = 0, alpha=0.25)+
geom_jitter(width=0.25)+
stat_summary(fun.y=mean, colour="black", geom="point",
shape=18, size=3) +
theme_classic()+
theme(legend.position="none")+
theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))+
geom_text(data = myletters_df, aes(label = letter, y = 1500 ), colour="black", size=5)

J’espère que cet article vous permettra de réaliser plus facilement vos figures illustrant des comparaisons de moyennes, notamment pour vos publications. Si vous avez d’autres astuces pour réaliser ce type de représentation graphique, indiquez les moi en commentaire.
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
Poursuivez votre lecture
- Comment retrouver sous R une couleur employée avec Excel ?
- Choisissez et modifiez très facilement les couleurs de vos graphs R avec Colour Picker
- Versionnage de vos scripts avec RStudio + Git
- Comment insérer des références bibliographiques dans un document Rmarkdown ?
- Estimer l’intervalle de confiance d’une moyenne par bootstrap en une seule ligne de commande
bonjour
Impossible de faire passer la fonction tri.to.squ
Une idée de la cause de ce disfonctionnement
Cordialement
Merci pour vos explications et pour l’aide fournie.
Je débute sous R et un certain nombre d’aspects restent abscons, qui m’empêchent de suivre le pas à pas (car il manque des pas pour moi).
J’essaye de faire une représentation graphique d’un test de wilcoxon où je compare 2 médianes.
Je cherche à reproduire le graphique avec les boxplot, les points, et l’affichage de la pvalue du test.
Le code que vous utilisez est le suivant :
p <- ggboxplot(ToothGrowth, x = "supp", y = "len", color = "supp", palette = "jco",
add = "jitter")
theme(legend.position="none")
# Add p-value
p stat_compare_means()
Mais je n'arrive pas à le mettre en œuvre : "objet p introuvable".
J'ai installé les package ggplot2 et ggpubr
j'ai du mal ensuite à m'y retrouver avec les significations des x et y que vous indiquez dans le code : je ne sais pas s'il s'agit de la formule à part entière ou du nom de vos variables (mais à ce moment là où indique t-on les variables à utiliser?).
Dans mon cas j'ai un nom de data qui contient l'ensemble de mes variables…. et les 2 variables que je croise dans mon test de wilcoxon. Mettons respectivement "data", "var1", "var2"…
Pourriez-vous svp m'éclairer pour cet exemple ? cela me permettrait ensuite de mettre en œuvre les autres. Encore merci pour votre patience
Bonjour donnez nous souvent les scripts en pdf pour qu’on puisse les télécharger et les utiliser
Vraiment pour les pistes données.
Cependant si on fait les stats avec le package PMCMR, devenu récemment PMCMRplus, il n’est plus possible de l’interfacer avec multcompLetters. PMCMRplus possède bien une fonction summaryGroup qui produit les lettres pour les groupes significatifs, mais sous forme de text. Je n’arrive pas à récupérer ces lettres pour les injecter dans ggplot. Une piste serait la fonction R capture.output…
Cordialement
Bonjour,
je ne connais pas du tout ce package…
Bonne continuation.
la fonction tri.to.squ ne veut pas fonctionner , est-ce que c’est pareil pour tout le monde ??
Salut.
Je crois que vous devez installer d’abord le package multcompView, ensuite exécutez la fonction tri.to.squ. Merci.
Bonjour,
Je vous remercie pour les explications et les codes. ça fonctionne parfaitement.
Je dois afficher les lettres pour plusieurs graphes générés en utilisant la fonction facet_wrap. Sauriez-vous m’aider ?
J’ai essayé plusieurs codes mais ça ne fonctionne pas…
Merci d’avance.
Bonjour Cécile,
non je n’ai pas de code tout prêt que je pourrais vous partager.
Bonne continuation.
Bonjour
Merci pour ses explications, j’essaye de représenter la moyenne, or à chaque fois c’est la médiane qui s’affiche, et comme elle est à zéro mes résultats ne sont pas joli…
Une idée de cet affichage?
Merci beaucoup
Bonne journée
Bonjour,
si vous représentez vos données avec la fonction ggboxplot() c’est la médiane qui est représentée. Et, à priori, il n’est pas possible d’ajouter la moyenne.
Mais vous pouvez regarder sur la page github du package : https://github.com/kassambara/ggpubr/issues.
Sinon, vous devez faire le boxplot avec ggplot2 et afficher la moyenne, comme avec cet exemple :
ggplot(data=PlantGrowth, aes(x=group, y=weight, fill=group))
geom_boxplot()
stat_summary(fun.y=mean, colour=”darkred”, geom=”point”,
shape=18, size=3,show_guide = FALSE)
Puis afficher les significativité à la main, , en suivant cet article : https://statistique-et-logiciel-r.com/ajouter-pvalues-sur-ggplot/
Bonne continuation.
Bonjour,
Merci pour ces explications très interessantes !
Je cherche à afficher l’asterisk de significativité uniquement au dessus du box plot de traitement lors d’une comparaison de moyennes entre deux groupes (groupe contrôle et groupe traitement).
Auriez vous une idée pour réaliser cela ?
Merci d’avance
Bonjour,
Je n’arrive pas à charger le package multcompView, quel miroir CRAN faut-il choisir ? J’ai essayé de télécharger aussi le package PMCMR mais impossible de le trouver dans les miroir CRAN français…
Merci pour vos réponses
Bonjour,
je viens de les télécharger sans problème depuis R Studio. Ma version de R est la 4.0.2.
Merci beaucoup, pour l’article, après plusieurs jours j’ai enfin réussi à afficher les lettres. Merci
Bonjour,
Je n’arrive pas à afficher le graphe final. En effet j’ai ce message d’erreur qui s’affiche : Error in FUN(X[[i]], …) : objet ‘Classe’ introuvable (où “Classe” est ma variable qui correspond à votre “group”). J’ai vérifié, il n’y a pas d’erreur de syntaxe, la variable “Classe” existe bel et bien dans ma base. Merci
Bonjour,
Est-ce que vous voyez bien votre variable Classe quand vous faites str(dataset) ?
Si oui, est-elle bien codée en facteur ? Si ce n’est pas le cas, il faudrait la passer en facteur avec : dataset$Classe <- as.factor(datset$Classe). J'espère que cela vous aide...
Bonjour,
j’ai le même problème, tout fonctionne sauf au moment de faire le ggplot où l’on m’affiche “Error in `geom_text()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 4th layer.
Caused by error in `FUN()`:
! object ‘nom’ not found
Run `rlang::last_error()` to see where the error occurred.”
J’ai fait ce que vous indiquez à la personne au dessus de moi mais j’ai toujours le message d’erreur…
Si quelqu’un peut m’aider…
Bonjour,
pourriez-vous s’il vous plait recopier vos lignes de commandes ?
Merci.
Bonjour,
J’ai un soucis pour afficher les lettres, je reçois cette erreur :
Erreur : Aesthetics must be either length 1 or the same as the data (4): x and fill
J’ai bien suivi les commandes à partir la création de l’anova, j’obtiens bien la table myletters_df, mais j’ai l’impression qu’il y a un conflit entre le ggplot principal qui utilise mon tableau de données principal qui a 21 lignes avec le second tableau de données myletters_df dans la fonction geomtext() qui lui possède 4 lignes. Je ne sais vraiment pas quoi faire.
Bonjour,
Merci pour ce tutoriel bien mené.
J’ai une question concernant la position des lettres. Ici vous fixez leur position par rapport à une valeur sur l’axe des y (1500 dans ce cas).
Est-ce qu’il y aurait un moyen d’afficher les lettres à des coordonnées y différentes pour chaque facteurs ? Et si oui, est-ce qu’il y a un moyen pour que R adapte la position de la lettre à chaque facteur tout seul (exemple : à n distance du quatrième quartile) ou les coordonnées devront être spécifiées pour chaque facteur ?
Merci d’avance,
Martin
Bonjour Martin,
il me semble qu’on peut appliquer des coordonnées différentes manuellement.
On peut sans doute la faire automatiquement, en calculant les positions en amont, et en les stockant dans des variables, qu’on utilise ensuite dans le graph
J’espère que cela vous aide un peu…
Bonne continuation.
Bonjour, Tous d’abord merci beaucoup pour vos articles très intéressant, ils m’ont aidé à de nombreuses reprise (étant étudiante). Cependant je rencontre un problème avec la fonction geom_text qui m’affiche cet erreur :
Error in `geom_bar()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 1st layer.
Caused by error in `FUN()`:
! object ‘group’ not found
Run `rlang::last_error()` to see where the error occurred.
toutes les lignes précédentes fonctionnent correctement avec votre script et j’obtient le même type de tableau pour les lettres en fonction des groupes mais une fois que je doit faire mon graphique cela ne fonctionne plus, merci d’avance si vous prendrez du temps à répondre à mon problème
Bonjour Manon,
je crois comprendre que dans vos données, il n’y a pas de variable “group”, ou alors elle se nomme autrement.
Tenez-moi informée.
Bonjour,
Désolé pour ma réponse tardive, j’ai réussi à résoudre mon problème en faisant :
p <- ggplot(histo, aes(x=group, y=test, fill = group)) +
geom_bar(stat="identity", position=position_dodge())+
geom_errorbar(aes(ymin=test-sd, ymax=test+sd), width=.2,
position=position_dodge(.9))
p+geom_text(aes(label=letter),position = position_dodge(0.9),
vjust = -15)
Car en effet j'essayais de rajouter les lettres mais sur un histogramme et pas un boxplot, je l'ai fais à partir d'un tableau histo contenant les moyennes, écart-type et lettre
Merci beaucoup pour votre réponse je comprend mieux pourquoi le script ne marchais pas
Bonjour, le boxplot avec les lettres de significativité fonctionne à merveille, j’aimerai maintenant affichier ce même graphique sans les points dans le boxplot. Savez vous comment faire ? De plus serait-il possible d’afficher les intervalle de confiance ?
Voici mes ligne de codes :
kruskal.test(CIRajc~Espece, data= fichierTravail2)
source(“./fonctions/fonction_tri_to_squ.R”)
pp <- pairwise.wilcox.test(fichierTravail2$CIRajc, fichierTravail2$Espece,p.adjust.method ="holm" )
pp
mymat <-tri.to.squ(pp$p.value)
mymat
library(multcompView)
myletters <- multcompLetters(mymat,compare="<=",threshold=0.05,Letters=letters)
myletters
myletters_df <- data.frame(Espece=names(myletters$Letters),letter = myletters$Letters )
myletters_df
ggplot(fichierTravail2, aes(x =Espece, y=CIRajc, colour=Espece, fill=Espece)) +
geom_boxplot(outlier.alpha = 0, alpha = 0) +
geom_jitter(width = 0.25) +
stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3) +
theme_classic() +
theme(legend.position = "none") +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust= 1)) +
geom_text(data = myletters_df, aes(x = Espece, label = letter, y = 1600), colour = "black", size = 5)
Bonjour,
pour ne pas afficher les points il faut supprimer la couche geom_jitter (vous pouvez mettre un diese davant).
Pour ajouter l’IC, vous pouvez peut être vous inspirer du code de cette page : http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/
Bonne continuation