R Graphics Essentials for Great Data Visualization: 200 Practical Examples You Want to Know for Data Science
NEW!!
In this article, we’ll describe how to easily i)
compare means
of two or multiple groups; ii) and to automatically add
p-values
and
significance levels
to a ggplot (such as box plots, dot plots, bar plots and line plots …).
Contents
Prerequisites
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
Load ggpubr:
library(ggpubr)
Official documentation of ggpubr is available at:
http://www.sthda.com/english/rpkgs/ggpubr
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:
Methods
R function
Description
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()
As we’ll show in the next sections, it has multiple useful options compared to the standard R functions.
The simplified format is as follow:
compare_means(formula, data, method = "wilcox.test", paired = FALSE,
group.by = NULL, ref.group = NULL, ...)
formula
: a formula of the form
x ~ group
, where x is a numeric variable and group is a factor with one or multiple levels. For example,
formula = TP53 ~ cancer_group
. It’s also possible to perform the test for multiple response variables at the same time. For example,
formula = c(TP53, PTEN) ~ cancer_group
.
data
: a data.frame containing the variables in the formula.
method
: the type of test. Default is
“wilcox.test”
. Allowed values include:
“t.test”
(parametric) and
“wilcox.test”
" (non-parametric). Perform comparison between two groups of samples. If the grouping variable contains more than two levels, then a pairwise comparison is performed.
“anova”
(parametric) and
“kruskal.test”
(non-parametric). Perform one-way ANOVA test comparing multiple groups.
ref.group
: a character string specifying the reference group. If specified, for a given grouping variable, each of the group levels will be compared to the reference group (i.e. control group). ref.group can be also
“.all.”
. In this case, each of the grouping variable levels is compared to all (i.e. base-mean).
stat_compare_means()
This function extends ggplot2 for adding mean comparison p-values to a ggplot, such as box blots, dot plots, bar plots and line plots.
The simplified format is as follow:
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
## 1 len OJ VC 0.0645 0.0645 0.064 ns Wilcoxon
By default
method = “wilcox.test”
(non-parametric test). You can also specify
method = “t.test”
for a parametric t-test.
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:
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)
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
## 1 len OJ VC 0.00431 0.00431 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)
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
## 1 len 9.53e-16 9.53e-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")
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
## 1 len 0.5 1 7.02e-06 1.40e-05 7.0e-06 **** Wilcoxon
## 2 len 0.5 2 8.41e-08 2.52e-07 8.4e-08 **** Wilcoxon
## 3 len 1 2 1.77e-04 1.77e-04 0.00018 *** Wilcoxon
# Visualize: Specify the comparisons you want
my_comparisons
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)
(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
## 1 len 0.5 1 6.70e-09 6.70e-09 6.7e-09 **** T-test
## 2 len 0.5 2 1.47e-16 2.94e-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
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
## 1 len .all. 0.5 1.24e-06 3.73e-06 1.2e-06 **** T-test
## 2 len .all. 1 5.67e-01 5.67e-01 0.57 ns T-test
## 3 len .all. 2 1.37e-05 2.74e-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
A typical situation, where pairwise comparisons against “all” can be useful, is illustrated here using the
myeloma
data set available on Github.
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 GitHub
myeloma
## # A tibble: 7 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## 1 DEPDC1 .all. Cyclin D-1 0.149690 0.44907 0.14969 ns T-test
## 2 DEPDC1 .all. Cyclin D-2 0.523143 1.00000 0.52314 ns T-test
## 3 DEPDC1 .all. Hyperdiploid 0.000282 0.00169 0.00028 *** T-test
## 4 DEPDC1 .all. Low bone disease 0.005084 0.02542 0.00508 ** T-test
## 5 DEPDC1 .all. MAF 0.086107 0.34443 0.08611 ns T-test
## 6 DEPDC1 .all. MMSET 0.576291 1.00000 0.57629 ns T-test
## # ... with 1 more rows
# 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
From the plot above, we can conclude that DEPDC1 is significantly overexpressed in proliferation group and, it’s significantly downexpressed in Hyperdiploid and Low bone disease compared to all.
# 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
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
## 1 0.5 len OJ VC 0.02319 0.0464 0.023 * Wilcoxon
## 2 1.0 len OJ VC 0.00403 0.0121 0.004 ** Wilcoxon
## 3 2.0 len OJ VC 1.00000 1.0000 1.000 ns Wilcoxon
In the example above, for each level of the variable “dose”, we compare the means of the variable “len” in the different groups formed by the grouping variable “supp”.
Visualize (1/2). Create a multi-panel box plots facetted by group (here, “dose”):
# Box plot facetted by "dose"
# Or use significance symbol as label
p + stat_compare_means(label = "p.signif", label.x = 1.5)
To hide the ‘ns’ symbol, use the argument hide.ns = TRUE.
Visualize (2/2). Create one single panel with all box plots. Plot y = “len” by x = “dose” and color by “supp”:
# Show only p-value
p + stat_compare_means(aes(group = supp), label = "p.format")
# Use significance symbol as label
p + stat_compare_means(aes(group = supp), label = "p.signif")
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
## 1 0.5 len OJ VC 0.0330 0.0659 0.033 * Wilcoxon
## 2 1.0 len OJ VC 0.0191 0.0572 0.019 * Wilcoxon
## 3 2.0 len OJ VC 1.0000 1.0000 1.000 ns Wilcoxon
Visualize. Create a multi-panel box plots facetted by group (here, “dose”):
# Box plot facetted by "dose"
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))
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))
Enjoyed this article? Give us 5 stars
(just above this text block)! Reader needs to be STHDA member for voting. I’d be very grateful if you’d help it spread by emailing it to a friend, or sharing it on Twitter, Facebook or Linked In.
Show me some love with the like buttons below... Thank you and please don't forget to share and comment below!!
Avez vous aimé cet article? Donnez nous 5 étoiles
(juste au dessus de ce block)! Vous devez être membre pour voter. Je vous serais très reconnaissant si vous aidiez à sa diffusion en l'envoyant par courriel à un ami ou en le partageant sur Twitter, Facebook ou Linked In.
Montrez-moi un peu d'amour avec les like ci-dessous ... Merci et n'oubliez pas, s'il vous plaît, de partager et de commenter ci-dessous!
Recommended for you
This section contains best data science and self-development resources to help you on your path.
Coursera - Online Courses and Specialization
Data science
Popular Courses Launched in 2020
Trending Courses
Books - Data Science
Our Books
Others
Comments
You are not authorized to post a comment
Comment
Can someone please help me how to get p-values of the Tukey test instead of anova above the different box plots?
thanks
Hello,
Is there a way to perform the ANOVA or Kruskal-Wallis and then follow-up with a Tukey HSD or Dunns' pairwise analysis instead of t-test or Wilcoxon pairwise test?
I know this option exists: https://rdrr.io/github/kassambara/ggpubr/man/stat_pvalue_manual.html but it would be so much more elegant to perform the test and have the values plotted as part of stat_compare.
Thanks
[img]https://imgur.com/nMW3iVS[/img] Is it possible to set p values separately in each group?
I have a data set with four separate treatment groups and I would like to have them all on the same graph with each group compared to treatment timepoint zero (there are 5 timepoints). Attached is an image of my data. I would like to have p.signif for each line (group) for timepoint (1,2,3,4,5) compared to my control timepoint (0). Is there a way to break up the stat_ into the individual groups (blue,black,red,green)? Thanks!
Is it possible to use the package ANOM (analysis of means) instead of t.test while comparing means against all? For the comparison of means against ".all.", this may be a more robust option than only using t.test without any correction, to reduce the chance of false positives.
Here is the line I am referring to:
stat_compare_means(label = "p.signif", method = "t.test",
ref.group = ".all.", hide.ns = TRUE)
Thanks!
When explicitly specifying comparisons
[e.g. stat_compare_means(comparisons = list( c("0.5","1"), c("1","2"), c("2","3")),label.y = c(15, 15.5, 16),paired = T) ]
- is it possible to:
1. define the number of digits after the decimal point
2. add "p=" before the p value (aes(label = paste0(“p =”, ..p.format..)) doesn't seem to work)
3. define that any value below 0.001 will appear as '>0.001' but still use the p.format notation
Solution:
I installed the most current version of ggpubr from github
devtools::install_github("kassambara/ggpubr")
Thean added to the stat_compare_means
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))
Thanks !
Hello,
thanks a lot for this amazing guide.
Is there also a way to get a compact letter display? So instead of the method you used to Compare more than two groups, could u also write the letters a, b, ab as significance indicators?
Best Felix
Hi Kassambara,
I'm having the same issue as Hao (post #594). Every time I try to create a plot specifying grouping by a color (color = "factor"), if, for example, I'm plotting a line plot, I get separate lines by groups, but if I then use add = "mean_se", it doesn't recognize the grouping variable (i.e. it will collapse across the factors and plot a single mean and SE).
This happened after I upgraded to R 3.5.1 on Mac OSX.
Hope you can help, this is a very useful package,
Is it also possible to compare both the significance of the supp OJ vs. VC for each of the doses (as above), but also for the significance between the plots per dose in the same graph?
So, significance for OJ vs. VS for doses 0.5, 1 and 2 AND significance for dose 0.5 vs. 1 vs 2 for OJ and VS?
Hi Kassambara,
Thank you very much for the codes ! For some reasons, the code in "Bar and line plots (two grouping variables)" section only produces one error bar in between the bars and lines instead of two error bars on both bard and lines. Could you please help? My ggplot2 and ggpubr are up-to-date.
Thanks a lot!
I get the following error on the first " p + stat_compare_means()"
Don't know how to automatically pick scale for object of type quosure/formula. Defaulting to continuous.
Don't know how to automatically pick scale for object of type quosure/formula. Defaulting to continuous.
Error in validDetails.text(x) :
'pairlist' object cannot be coerced to type 'double'
Is there a way to resolve this. I am using R x64 3.5.0
Hi everybody,
I tried the the stat_compare_means from the myeloma data. But got the following error:
Don't know how to automatically pick scale for object of type quosure/formula. Defaulting to continuous.
Don't know how to automatically pick scale for object of type quosure/formula. Defaulting to continuous.
Error: length(rows) == 1 ist nicht TRUE
I want to use the funtion for my own data but I get the same error.
Who can help me?
Thanks
