Prerequisites

install.packages("ggpubr")

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

library(ggpubr)

data("ToothGrowth")
head(ToothGrowth)
   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

Methods for comparing means

The standard methods to compare the means of two or more groups in R, have been largely described at: comparing means in R.

The most common methods for comparing means include:

A practical guide to compute and interpret the results of each of these methods are provided at the following links:

R functions to add p-values

Here we present two new R functions in the ggpubr package:

• compare_means(): easy to use solution to performs one and multiple mean comparisons.
• stat_compare_means(): easy to use solution to automatically add p-values and significance levels to a ggplot.

compare_means(formula, data, method = "wilcox.test", paired = FALSE,
  group.by = NULL, ref.group = NULL, ...)

stat_compare_means(mapping = NULL, comparisons = NULL hide.ns = FALSE,
                   label = NULL,  label.x = NULL, label.y = NULL,  ...)

Compare two independent groups

Perform the test:

compare_means(len ~ supp, data = ToothGrowth)
# A tibble: 1 x 8
    .y. group1 group2          p      p.adj p.format p.signif   method
  <chr>  <chr>  <chr>      <dbl>      <dbl>    <chr>    <chr>    <chr>
1   len     OJ     VC 0.06449067 0.06449067    0.064       ns Wilcoxon

Returned value is a data frame with the following columns:

• .y.: the y variable used in the test.
• p: the p-value
• p.adj: the adjusted p-value. Default value for p.adjust.method = “holm”
• p.format: the formatted p-value
• p.signif: the significance level.
• method: the statistical test used to compare groups.

Create a box plot with p-values:

p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
          color = "supp", palette = "jco",
          add = "jitter")
#  Add p-value
p + stat_compare_means()

# Change method
p + stat_compare_means(method = "t.test")

Add p-values and significance levels to ggplots

Note that, the p-value label position can be adjusted using the arguments: label.x, label.y, hjust and vjust.

The default p-value label displayed is obtained by concatenating the method and the p columns of the returned data frame by the function compare_means(). You can specify other combinations using the aes() function.

For example,

• aes(label = ..p.format..) or aes(label = paste0(“p =”, ..p.format..)): display only the formatted p-value (without the method name)
• aes(label = ..p.signif..): display only the significance level.
• aes(label = paste0(..method.., “\n”, “p =”, ..p.format..)): Use line break (“\n”) between the method name and the p-value.

As an illustration, type this:

p + stat_compare_means( aes(label = ..p.signif..), 
                        label.x = 1.5, label.y = 40)

Add p-values and significance levels to ggplots

If you prefer, it’s also possible to specify the argument label as a character vector:

p + stat_compare_means( label = "p.signif", label.x = 1.5, label.y = 40)

Compare two paired samples

Perform the test:

compare_means(len ~ supp, data = ToothGrowth, paired = TRUE)
# A tibble: 1 x 8
    .y. group1 group2           p       p.adj p.format p.signif   method
  <chr>  <chr>  <chr>       <dbl>       <dbl>    <chr>    <chr>    <chr>
1   len     OJ     VC 0.004312554 0.004312554   0.0043       ** Wilcoxon

Visualize paired data using the ggpaired() function:

ggpaired(ToothGrowth, x = "supp", y = "len",
         color = "supp", line.color = "gray", line.size = 0.4,
         palette = "jco")+
  stat_compare_means(paired = TRUE)

Add p-values and significance levels to ggplots

Compare more than two groups

• Global test:

# Global test
compare_means(len ~ dose,  data = ToothGrowth, method = "anova")
# A tibble: 1 x 6
    .y.            p        p.adj p.format p.signif method
  <chr>        <dbl>        <dbl>    <chr>    <chr>  <chr>
1   len 9.532727e-16 9.532727e-16  9.5e-16     ****  Anova

Plot with global p-value:

# Default method = "kruskal.test" for multiple groups
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means()

# Change method to anova
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means(method = "anova")

Add p-values and significance levels to ggplots

• Pairwise comparisons. If the grouping variable contains more than two levels, then pairwise tests will be performed automatically. The default method is “wilcox.test”. You can change this to “t.test”.

# Perorm pairwise comparisons
compare_means(len ~ dose,  data = ToothGrowth)
# A tibble: 3 x 8
    .y. group1 group2            p        p.adj p.format p.signif   method
  <chr>  <chr>  <chr>        <dbl>        <dbl>    <chr>    <chr>    <chr>
1   len    0.5      1 7.020855e-06 1.404171e-05  7.0e-06     **** Wilcoxon
2   len    0.5      2 8.406447e-08 2.521934e-07  8.4e-08     **** Wilcoxon
3   len      1      2 1.772382e-04 1.772382e-04  0.00018      *** Wilcoxon
# Visualize: Specify the comparisons you want
my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)     # Add global p-value

Add p-values and significance levels to ggplots

If you want to specify the precise y location of bars, use the argument label.y:

ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons, label.y = c(29, 35, 40))+
  stat_compare_means(label.y = 45)

Add p-values and significance levels to ggplots

(Adding bars, connecting compared groups, has been facilitated by the ggsignif R package )

• Multiple pairwise tests against a reference group:

# Pairwise comparison against reference
compare_means(len ~ dose,  data = ToothGrowth, ref.group = "0.5",
              method = "t.test")
# A tibble: 2 x 8
    .y. group1 group2            p        p.adj p.format p.signif method
  <chr>  <chr>  <chr>        <dbl>        <dbl>    <chr>    <chr>  <chr>
1   len    0.5      1 6.697250e-09 6.697250e-09  6.7e-09     **** T-test
2   len    0.5      2 1.469534e-16 2.939068e-16  < 2e-16     **** T-test
# Visualize
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means(method = "anova", label.y = 40)+      # Add global p-value
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = "0.5")                    # Pairwise comparison against reference

Add p-values and significance levels to ggplots

• Multiple pairwise tests against all (base-mean):

# Comparison of each group against base-mean
compare_means(len ~ dose,  data = ToothGrowth, ref.group = ".all.",
              method = "t.test")
# A tibble: 3 x 8
    .y. group1 group2            p        p.adj p.format p.signif method
  <chr>  <chr>  <chr>        <dbl>        <dbl>    <chr>    <chr>  <chr>
1   len  .all.    0.5 1.244626e-06 3.733877e-06  1.2e-06     **** T-test
2   len  .all.      1 5.667266e-01 5.667266e-01     0.57       ns T-test
3   len  .all.      2 1.371704e-05 2.743408e-05  1.4e-05     **** T-test
# Visualize
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means(method = "anova", label.y = 40)+      # Add global p-value
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = ".all.")                  # Pairwise comparison against all

Add p-values and significance levels to ggplots

A typical situation, where pairwise comparisons against “all” can be useful, is illustrated here using the myeloma data set from the survminer package.

We’ll plot the expression profile of the DEPDC1 gene according to the patients’ molecular groups. We want to know if there is any difference between groups. If yes, where the difference is?

To answer to this question, you can perform a pairwise comparison between all the 7 groups. This will lead to a lot of comparisons between all possible combinations. If you have many groups, as here, it might be difficult to interpret.

Another easy solution is to compare each of the seven groups against “all” (i.e. base-mean). When the test is significant, then you can conclude that DEPDC1 is significantly overexpressed or downexpressed in a group xxx compared to all.

# Load myeloma data from survminer package
if(!require(survminer)) install.packages("survminer")
data("myeloma", package = "survminer")

# Perform the test
compare_means(DEPDC1 ~ molecular_group,  data = myeloma,
              ref.group = ".all.", method = "t.test")
# A tibble: 7 x 8
     .y. group1           group2            p        p.adj p.format p.signif method
   <chr>  <chr>            <chr>        <dbl>        <dbl>    <chr>    <chr>  <chr>
1 DEPDC1  .all.       Cyclin D-1 1.496896e-01 4.490687e-01  0.14969       ns T-test
2 DEPDC1  .all.       Cyclin D-2 5.231428e-01 1.000000e+00  0.52314       ns T-test
3 DEPDC1  .all.     Hyperdiploid 2.815333e-04 1.689200e-03  0.00028      *** T-test
4 DEPDC1  .all. Low bone disease 5.083985e-03 2.541992e-02  0.00508       ** T-test
5 DEPDC1  .all.              MAF 8.610664e-02 3.444265e-01  0.08611       ns T-test
6 DEPDC1  .all.            MMSET 5.762908e-01 1.000000e+00  0.57629       ns T-test
7 DEPDC1  .all.    Proliferation 1.241416e-09 8.689910e-09  1.2e-09     **** T-test
# Visualize the expression profile
ggboxplot(myeloma, x = "molecular_group", y = "DEPDC1", color = "molecular_group", 
          add = "jitter", legend = "none") +
  rotate_x_text(angle = 45)+
  geom_hline(yintercept = mean(myeloma$DEPDC1), linetype = 2)+ # Add horizontal line at base mean
  stat_compare_means(method = "anova", label.y = 1600)+        # Add global annova p-value
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = ".all.")                      # Pairwise comparison against all

Add p-values and significance levels to ggplots

# Visualize the expression profile
ggboxplot(myeloma, x = "molecular_group", y = "DEPDC1", color = "molecular_group", 
          add = "jitter", legend = "none") +
  rotate_x_text(angle = 45)+
  geom_hline(yintercept = mean(myeloma$DEPDC1), linetype = 2)+ # Add horizontal line at base mean
  stat_compare_means(method = "anova", label.y = 1600)+        # Add global annova p-value
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = ".all.", hide.ns = TRUE)      # Pairwise comparison against all

Add p-values and significance levels to ggplots

Multiple grouping variables

• Two independent sample comparisons after grouping the data by another variable:

Perform the test:

compare_means(len ~ supp, data = ToothGrowth, 
              group.by = "dose")
# A tibble: 3 x 9
   dose   .y. group1 group2           p      p.adj p.format p.signif   method
  <dbl> <chr>  <chr>  <chr>       <dbl>      <dbl>    <chr>    <chr>    <chr>
1   0.5   len     OJ     VC 0.023186427 0.04637285    0.023        * Wilcoxon
2   1.0   len     OJ     VC 0.004030367 0.01209110    0.004       ** Wilcoxon
3   2.0   len     OJ     VC 1.000000000 1.00000000    1.000       ns Wilcoxon

Visualize (1/2). Create a multi-panel box plots facetted by group (here, “dose”):

# Box plot facetted by "dose"
p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
          color = "supp", palette = "jco",
          add = "jitter",
          facet.by = "dose", short.panel.labs = FALSE)
# Use only p.format as label. Remove method name.
p + stat_compare_means(label = "p.format")

Add p-values and significance levels to ggplots

# Or use significance symbol as label
p + stat_compare_means(label =  "p.signif", label.x = 1.5)

Add p-values and significance levels to ggplots

Visualize (2/2). Create one single panel with all box plots. Plot y = “len” by x = “dose” and color by “supp”:

p <- ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "supp", palette = "jco",
          add = "jitter")
p + stat_compare_means(aes(group = supp))

Add p-values and significance levels to ggplots

# Show only p-value
p + stat_compare_means(aes(group = supp), label = "p.format")

Add p-values and significance levels to ggplots

# Use significance symbol as label
p + stat_compare_means(aes(group = supp), label = "p.signif")

Add p-values and significance levels to ggplots

• Paired sample comparisons after grouping the data by another variable:

Perform the test:

compare_means(len ~ supp, data = ToothGrowth, 
              group.by = "dose", paired = TRUE)
# A tibble: 3 x 9
   dose   .y. group1 group2          p      p.adj p.format p.signif   method
  <dbl> <chr>  <chr>  <chr>      <dbl>      <dbl>    <chr>    <chr>    <chr>
1   0.5   len     OJ     VC 0.03296938 0.06593876    0.033        * Wilcoxon
2   1.0   len     OJ     VC 0.01905889 0.05717667    0.019        * Wilcoxon
3   2.0   len     OJ     VC 1.00000000 1.00000000    1.000       ns Wilcoxon

Visualize. Create a multi-panel box plots facetted by group (here, “dose”):

# Box plot facetted by "dose"
p <- ggpaired(ToothGrowth, x = "supp", y = "len",
          color = "supp", palette = "jco", 
          line.color = "gray", line.size = 0.4,
          facet.by = "dose", short.panel.labs = FALSE)
# Use only p.format as label. Remove method name.
p + stat_compare_means(label = "p.format", paired = TRUE)

Add p-values and significance levels to ggplots

Other plot types

• Bar and line plots (one grouping variable):

# Bar plot of mean +/-se
ggbarplot(ToothGrowth, x = "dose", y = "len", add = "mean_se")+
  stat_compare_means() +                                         # Global p-value
  stat_compare_means(ref.group = "0.5", label = "p.signif",
                     label.y = c(22, 29))                   # compare to ref.group

# Line plot of mean +/-se
ggline(ToothGrowth, x = "dose", y = "len", add = "mean_se")+
  stat_compare_means() +                                         # Global p-value
  stat_compare_means(ref.group = "0.5", label = "p.signif",
                     label.y = c(22, 29))     

Add p-values and significance levels to ggplots

• Bar and line plots (two grouping variables):

ggbarplot(ToothGrowth, x = "dose", y = "len", add = "mean_se",
          color = "supp", palette = "jco", 
          position = position_dodge(0.8))+
  stat_compare_means(aes(group = supp), label = "p.signif", label.y = 29)

ggline(ToothGrowth, x = "dose", y = "len", add = "mean_se",
          color = "supp", palette = "jco")+
  stat_compare_means(aes(group = supp), label = "p.signif", 
                     label.y = c(16, 25, 29))

Add p-values and significance levels to ggplots

Infos

This analysis has been performed using R software (ver. 3.3.2) and ggpubr (ver. 0.1.3).