p-values ajustées : affichez les facilement sur un graph

Contexte

Il y a quelque temps, j’ai écrit un article expliquant comment, après une ANOVA, indiquer les  p-values ou les différences significatives sur un graph , en employant le package `ggpubr`.

 

Par exemple, pour une ANOVA non paramétrique (test de Kruskal Wallis) :

Je disais qu’une des limites du package était que les p-values des comparaisons multiples, affichées sur le graph, étaient les p-values non ajustéesOr quand on réalise une ANOVA (paramétrique ou non paramétrique), les p-values des comparaisons multiples doivent être
ajustées. Et donc que se sont ces p-values ajustées qu’il faudrait pouvoir afficher.

Pour plus d’informations sur « pourquoi il faut ajuster les p-values des comparaisons multiples », consultez l’ article : « Comparaisons multiples et ajustement des pvalues avec le logiciel R« 

Et puis la semaine dernière, un lecteur du blog m’a demandé si j’avais finalement trouvé une solution pour afficher ces p-values ajustées.

Je ne m’y étais plus intéressée alors en en première approche je lui ai dit de faire lui même les traits sous ggplot2 avec la fonction `geom_segment()`, puis ajouter les p-values avec la fonction `annotate()`

Et puis, en cherchant un peu plus, j’ai trouvé une solution qui avait été postée par un utilisateur du package `ggpubr` dans l’onglet « Issue » du répertoire GitHub du packageC’est cette solution que je vous explique ici en pas à pas.

Le Principe

La solution consiste à :

1. Indiquer les comparaisons deux à deux souhaitées sous forme d’une list

2. Créer un data frame qui contiendra, les p-values ajustées en utilisant la fonction compare_means(), et en précisant :

  • la méthode d’ajustement souhaitée, grâce à l’argument `p.adj`
  • les tests à employer pour les comparaisons multiples, via l’argument `method=`. On utilise `method=’t.test’` pour des tests t et `method=’wilcox.test’` pour des tests de Wilcoxon.

3. Ajouter, à ce data frame, une variable y.position dans laquelle on précise le futur emplacement, sur le graph, au niveau de l’axe des ordonnées, des p-values des comparaisons multiples

4. Faire le graph

P-values ajustées après une ANOVA paramétrique

Réalisation

J’utilise ici le jeu de données ToothGrowth (contenu dans le package dataset, chargé par défaut !) :

data("ToothGrowth")
TG <- ToothGrowth

head(TG)

    ##    len supp dose
    ## 1  4.2   VC  0.5
    ## 2 11.5   VC  0.5
    ## 3  7.3   VC  0.5
    ## 4  5.8   VC  0.5
    ## 5  6.4   VC  0.5
    ## 6 10.0   VC  0.5

str(TG)

    ## 'data.frame':    60 obs. of  3 variables:
    ##  $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
    ##  $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
    ##  $ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ... 

Je commence par passer la variable dose en facteur (par défaut, R la considère comme une variable numérique) :

TG$dose <- as.factor(TG$dose)

str(TG)

    ## 'data.frame':    60 obs. of  3 variables:
    ##  $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
    ##  $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
    ##  $ dose: Factor w/ 3 levels "0.5","1","2": 1 1 1 1 1 1 1 1 1 1 ... 

On réalise les 3 étapes décrites précédemment :

Définir les comparaisons deux à deux souhaitées

# on indique les comparaisons deux à deux
my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") ) 

Créer un data frame qui contiendra les pvalues

Pour cela, nous utilisons la fonction compare_means(),  en précisant la méthode d’ajustement souhaitée, et le type de tests.

Ici, j’appelle mc_t_test le data frame qui contiendra les informations des comparaisons et notamment les p-values ajustées. Je choisis d’utiliser la méthode d’ajustement d’Holm. D’autres choix sont possibles (cf. paragraphe 5). Et comme il s’agit d’une ANOVA paramétrique, je précise , dans l’argument method que les tests à employer sont des tests t.

mc_t_test <- compare_means(len ~ dose, comparisons = my_comparisons, p.adj= "holm", method='t.test', data = ToothGrowth)

    mc_t_test

    ## # A tibble: 3 x 8
    ##   .y.   group1 group2        p    p.adj p.format p.signif method
    ##                         
    ## 1 len   0.5    1      1.27e- 7 2.50e- 7 1.3e-07  ****     T-test
    ## 2 len   0.5    2      4.40e-14 1.30e-13 4.4e-14  ****     T-test
    ## 3 len   1      2      1.91e- 5 1.90e- 5 1.9e-05  ****     T-test 

Ajouter une variable y.position dans ce data frame

Cette variable  `y.position` permet de préciser le futur emplacement, sur le graph, au niveau de l’axe des ordonnées, des p-values des comparaisons multiples:

Je choisis les valeurs y=29 ,y=35, et y=39 en m’aidant du premier graph réalisé dans l’introduction.

# on ajoute la hauteur ou indiquer les p-values sur le graphs
mc_t_test <- mc_t_test %>%
        mutate(y.position = c(29, 35, 39))

mc_t_test

    ## # A tibble: 3 x 9
    ##   .y.   group1 group2        p    p.adj p.format p.signif method y.position
    ##                               
    ## 1 len   0.5    1      1.27e- 7 2.50e- 7 1.3e-07  ****     T-test         29
    ## 2 len   0.5    2      4.40e-14 1.30e-13 4.4e-14  ****     T-test         35
    ## 3 len   1      2      1.91e- 5 1.90e- 5 1.9e-05  ****     T-test         39 

Pour plus d’infos sur la création d’une variable à l’aide de la fonction mutate() et de l’utilisation du symbole pipe %>%, consultez l’article « Initiation à la manipulation de données avec le package dplyr « 

Faire le graph

On indique method = "anova" dans la fonction stat_compare_means() et on indique le data frame mc_t_test dans la fonction stat_pvalue_manual() :

ggboxplot(TG, x = "dose", y = "len", color = "dose",
        palette =c("#00AFBB", "#E7B800", "#FC4E07"),
        add = "jitter", shape = "dose")+
    
    # ajout de la p-value globale
    stat_compare_means(method = "anova", label.y=45)+
    
     # ajout des p-values des comparaisons deux à deux grace à la fonction
     # stat_pvalue_manual et l'objet mc_t_test
     stat_pvalue_manual(mc_t_test, label = "p.adj")+
     theme(legend.position="none") 

Vérifications

Réalisation de l'ANOVA


my.anova <- lm(len~dose, data=TG)
 
    library(car)
    Anova(my.anova)
 
    ## Anova Table (Type II tests)
    ## 
    ## Response: len
    ##           Sum Sq Df F value    Pr(>F)    
    ## dose      2426.4  2  67.416 9.533e-16 ***
    ## Residuals 1025.8 57                      
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

La p-value est bien identique à celle du graph.

Réalisation des comparaisons multiples

Lorsque l’argument method='t.test' est utilisé dans la fonction compare_mean() du package ggpubr , les tests t employés sont des tests pour variances inégales. Ces tests t sont aussi appelés tests de Welch.
Il est possible de réaliser tous les tests t deux à deux, en utilisant la fonction pairwise.t.test. Pour obtenir des tests de Welch, on utilise l’argument pool.sd =FALSE.

pairwise.t.test(TG$len, TG$dose, p.adj="holm", pool.sd =FALSE)

    ## 
    ##  Pairwise comparisons using t tests with non-pooled SD 
    ## 
    ## data:  TG$len and TG$dose 
    ## 
    ##   0.5     1      
    ## 1 2.5e-07 -      
    ## 2 1.3e-13 1.9e-05
    ## 
    ## P value adjustment method: holm 

Lorsque l’argument method='t.test' est utilisé dans la fonction compare_mean() du package ggpubr , les tests t employés sont des tests pour variances inégales. Ces tests t sont aussi appelés tests de Welch.
Il est possible de réaliser tous les tests t deux à deux, en utilisant la fonction pairwise.t.test. Pour obtenir des tests de Welch, on utilise l’argument pool.sd =FALSE.

Ces p-values correspondent bien à celle du graph.

P-values ajustées après une ANOVA non paramétrique

L’ANOVA non paramétrique est appelée test de Kruskal Wallis.

Réalisation

Définir les comparaisons deux à deux

On conserve les comparaisons deux à deux précédentes:

 # on indique les comparaisons deux à deux
    my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )
    
mc_wilcox_test <- compare_means(len ~ dose, comparisons = my_comparisons, p.adj= "holm", method='wilcox.test', data = ToothGrowth)     

Créer un data frame qui contiendra les p-values ajustées

Pour cela, nous utilisons la fonction compare_means()`, en  précisant la méthode d’ajustement souhaitée, et le type de tests.

Ici j’appelle `mc_wilcox_test` le data frame qui contiendra les informations des comparaisons et notamment les p-values ajustées. Comme précédemment, je choisis d’utiliser la méthode d’ajustement d’Holm. D’autres choix sont possibles (paragraphe 5 ). Et comme il s’agit d’une ANOVA non paramétrique, je précise que les tests à employer sont des tests de wilcoxon.

 mc_wilcox_test

    ## # A tibble: 3 x 8
    ##   .y.   group1 group2            p      p.adj p.format p.signif method  
    ##                                 
    ## 1 len   0.5    1      0.00000702   0.000014   7.0e-06  ****     Wilcoxon
    ## 2 len   0.5    2      0.0000000841 0.00000025 8.4e-08  ****     Wilcoxon
    ## 3 len   1      2      0.000177     0.00018    0.00018  ***      Wilcoxon 

Ajouter une variable y.position dans ce data frame

Je garde les même valeurs que précédemment : y=29 ,y=35, et y=39

# on ajoute la hauteur ou indiquer les p-values sur le graphs
    mc_wilcox_test <- mc_wilcox_test %>%
        mutate(y.position = c(29, 35, 39))

    mc_wilcox_test

    ## # A tibble: 3 x 9
    ##   .y.   group1 group2        p    p.adj p.format p.signif method y.position
    ##                               
    ## 1 len   0.5    1       7.02e-6  1.40e-5 7.0e-06  ****     Wilco~         29
    ## 2 len   0.5    2       8.41e-8  2.50e-7 8.4e-08  ****     Wilco~         35
    ## 3 len   1      2       1.77e-4  1.80e-4 0.00018  ***      Wilco~         39 

Faire le graph

On indique method = "kruskal.test" dans la fonction stat_compare_means() et on indique le data frame mc_wilcox_test dans la fonction stat_pvalue_manual() 

ggboxplot(TG, x = "dose", y = "len", color = "dose", 
              palette =c("#00AFBB", "#E7B800", "#FC4E07"),
              add = "jitter", shape = "dose")+
              # ajout de la p-value globale
              stat_compare_means(method = "kruskal.test", label.y=45)+
              # ajout des p-values des comparaisons deux à deux grace à la fonction
              # stat_pvalue_manual et l'objet mc_wilcox_test
             stat_pvalue_manual(mc_wilcox_test, label = "p.adj")+
             theme(legend.position="none") 

Vérification

Réalisation du test de Kruskal Wallis

my.ks <- kruskal.test(len~dose, data=TG)
    my.ks

    ## 
    ##  Kruskal-Wallis rank sum test
    ## 
    ## data:  len by dose
    ## Kruskal-Wallis chi-squared = 40.669, df = 2, p-value = 1.475e-09 

La p-value est bien identique à celle du graph.

Réalisation des comparaisons multiples

Il est possible de réaliser tous les tests de Wilcoxon deux à deux, en utilisant la fonction pairwise.wilcox.test :

pairwise.wilcox.test(TG$len, TG$dose, p.adj="holm")

    ## 
    ##  Pairwise comparisons using Wilcoxon rank sum test 
    ## 
    ## data:  TG$len and TG$dose 
    ## 
    ##   0.5     1      
    ## 1 1.4e-05 -      
    ## 2 2.5e-07 0.00018
    ## 
    ## P value adjustment method: holm 

Ces p-values correspondent bien à celle du graph.

Remarques : les messages de warning indiquent simplement qu’il existe des valeurs de `len` identiques dans plusieurs groupes et que cela rend impossible le calcule d’une p-value exacte. Une approximation sera alors utilisée (ce n’est pas grave !).

Remarques

Il existe plusieurs méthodes d’ajustement des p-values qui peuvent être employées dans les fonctions compare_means()pairwise.t.test() et pairwise.wilcox.test(), via l’arugment p.adj.

Pour connaitre comment préciser ces méthodes, utiliser la commande suivante :

p.adjust.methods

    ## [1] "holm"       "hochberg"   "hommel"     "bonferroni" "BH"        
    ## [6] "BY"         "fdr"        "none" 

Pour plus d’information, consulter la page d’aide, en utilisant la commande :

?p.adjust 

Pour aller plus loin

Je vous conseille vivement d’aller faire un tour sur la page de présentation du package `ggpubr`car ce package propose de nombreuses options de graphes très intéressantes.Et si vous rencontrez des difficultés, consultez l’onglet Issue, c’est là que les utilisateurs posent leur questions, remontent les bugs rencontrés, ou encore postent des réponses !Si vous souhaitez plus d’informations, théoriques ou pratiques, sur l’ANOVA à un facteur, vous pouvez consulter mes articles : Et vous, avez vous d’autres astuces, ou d’autres routines pour présenter vos résultats d’ANOVA à un facteur ? 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 🙏

Crédits photos : geralt

13 Responses

  1. Je trouve cela particulièrement intéressant, car, personnellement, je ne trouvais une pas grande utilité à comparer des boxplots. Mais avec les p-values indiqués comme ceci ça a une vraie valeur ajoutée !

    Merci Claire

  2. Bonjour, tout d’abord merci pour l’article.
    J’aimerai savoir s’il est possible de spécifier le nombre de comparaison que l’on souhaite faire dans le p.adj. Car spécifier que la correction se fait avec Holms est une chose mais j’aimerai être sur que le nombre de comparaisons est exactement celui que je souhaite.
    Bine à vous,
    MP

    1. Bonjour,

      On ne précise pas directement le nombre de comparaisons mais on spécifie les comparaisons individuellement via
      my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )et ensuite my_comparaisons passe en argument de stat_compare_means(). Bonne continuation

    1. Et ben, je pense qu’il faut faire les segments horizontaux sous ggplot2 avec geom_segment() et ajouter les p-values ajustées avec annotate(), mais c’est fastidieux. J’essaierai de le montrer dans un prochain article.
      Bonne continuation

  3. Bonjour, Je vous remercie pour cet article vraiment intéressant. J’avais une question pour le test de kruskal Wallis, comment est-ce qu’on pourrait procéder quand on se retrouve dans un cas avec un rang, au lieu de plusieurs.

Laisser un commentaire

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

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.