R ggplot2 boxplots - ggpubr stat_compare_means not working properly

r ggplot2 boxplot p-value ggpubr
DaniCee picture DaniCee · Aug 3, 2017 · Viewed 10.2k times · Source

I am trying to add significance levels to my boxplots in the form of asterisks using ggplot2 and the ggpubr package, but I have many comparisons and I only want to show the significant ones.

I try to use the option hide.ns=TRUE in stat_compare_means, but it clearly does not work, it might be a bug in the ggpubr package.

Besides, you see that I leave out group "PGMC4" from the pairwise wilcox.test comparisons; how can I leave this group out also for the kruskal.test?

The last question I have is how the significance level works? As in * is significant below 0.05, ** below 0.025, *** below 0.01? what is the convention ggpubr uses? Is it showing p-values or adjusted p-values? If the latter, what's the adjusting method? BH?

Please check my MWE below and this link and this other one for reference

##############################
##MWE
set.seed(5)
#test df
mydf <- data.frame(ID=paste(sample(LETTERS, 163, replace=TRUE), sample(1:1000, 163, replace=FALSE), sep=''),
                   Group=c(rep('C',10),rep('FH',10),rep('I',19),rep('IF',42),rep('NA',14),rep('NF',42),rep('NI',15),rep('NS',10),rep('PGMC4',1)),
                   Value=rnorm(n=163))
#I don't want to compare PGMC4 cause I have only onw sample
groups <- as.character(unique(mydf$Group[which(mydf$Group!="PGMC4")]))
#function to make combinations of groups without repeating pairs, and avoiding self-combinations
expand.grid.unique <- function(x, y, include.equals=FALSE){
    x <- unique(x)
    y <- unique(y)
    g <- function(i){
        z <- setdiff(y, x[seq_len(i-include.equals)])
        if(length(z)) cbind(x[i], z, deparse.level=0)
    }
    do.call(rbind, lapply(seq_along(x), g))
}
#all pairs I want to compare
combs <- as.data.frame(expand.grid.unique(groups, groups), stringsAsFactors=FALSE)
head(combs)
my.comps <- as.data.frame(t(combs), stringsAsFactors=FALSE)
colnames(my.comps) <- NULL
rownames(my.comps) <- NULL
#pairs I want to compare in list format for stat_compare_means
my.comps <- as.list(my.comps)
head(my.comps)
pdf(file="test.pdf", height=20, width=25)
print(#or ggsave()
  ggplot(mydf, aes(x=Group, y=Value, fill=Group)) + geom_boxplot() +
    stat_summary(fun.y=mean, geom="point", shape=5, size=4) +
    scale_fill_manual(values=myPal) +
    ggtitle("TEST TITLE") +
    theme(plot.title = element_text(size=30),
      axis.text=element_text(size=12),
      axis.text.x = element_text(angle=45, hjust=1),
      axis.ticks = element_blank(),
      axis.title=element_text(size=20,face="bold"),
      legend.text=element_text(size=16)) +
  stat_compare_means(comparisons=my.comps, method="wilcox.test", label="p.signif", size=14) + #WHY DOES hide.ns=TRUE NOT WORK??? WHY DOES size=14 NOT WORK???
  stat_compare_means(method="kruskal.test", size=14) #GLOBAL COMPARISON ACROSS GROUPS (HOW TO LEAVE PGMC4 OUT OF THIS??)
)
dev.off()
##############################

The MWE will produce the following boxplots:

test

The questions would be:

1- How to make hide.ns=TRUE work?

2- How to increase the size of the *?

3- How to exclude a group from the kruskal.test comparison?

4- What is the * convention used by ggpubr, and are the p-values shown adjusted or not?

Many thanks!!

EDIT

Besides, when doing

stat_compare_means(comparisons=my.comps, method="wilcox.test", p.adjust.method="BH")

I do not obtain the same p-values as when doing

wilcox.test(Value ~ Group, data=mydf.sub)$p.value

where mydf.sub is a subset() of mydf for a given comparison of 2 groups.

What is ggpubr doing here? How does it calculate the p.values?

EDIT 2

Please help, the solution does not have to be with ggpubr (but it has to be with ggplot2), I just need to be able to hide the NS and make the size of the asterisks bigger, as well as a p-value calculation identical to wilcox.test() + p.adjust(method"BH").

Thanks!

Answer

Roman picture Roman · Aug 7, 2017

Edit: Since I discovered the rstatix package I would do:

set.seed(123)
#test df
mydf <- data.frame(ID=paste(sample(LETTERS, 163, replace=TRUE), sample(1:1000, 163, replace=FALSE), sep=''),
                   Group=c(rep('C',10),rep('FH',10),rep('I',19),rep('IF',42),rep('NA',14),rep('NF',42),rep('NI',15),rep('NS',10),rep('PGMC4',1)),
                   Value=c(runif(n=100), runif(63,max= 0.5)))


library(tidyverse)

stat_pvalue <- mydf %>% 
 rstatix::wilcox_test(Value ~ Group) %>%
 filter(p < 0.05) %>% 
 rstatix::add_significance("p") %>% 
 rstatix::add_y_position() %>% 
 mutate(y.position = seq(min(y.position), max(y.position),length.out = n())

ggplot(mydf, aes(x=Group, y=Value)) + geom_boxplot() +
  ggpubr::stat_pvalue_manual(stat_pvalue, label = "p.signif") +
  theme_bw(base_size = 16)

enter image description here

Old Answer:

You can try following. The idea is that you calculate the stats by your own using pairwise.wilcox.test. Then you use the ggsignif function geom_signif to add the precalculated pvalues. With y_position you can place the brackets so they don't overlap.

library(tidyverse)
library(ggsignif)
library(broom)
# your list of combinations you want to compare
CN <- combn(levels(mydf$Group)[-9], 2, simplify = FALSE)
# the pvalues. I use broom and tidy to get a nice formatted dataframe. Note, I turned off the adjustment of the pvalues. 
pv <- tidy(with(mydf[ mydf$Group != "PGMC4", ], pairwise.wilcox.test(Value, Group, p.adjust.method = "none")))
#  data preparation 
CN2 <- do.call(rbind.data.frame, CN)
colnames(CN2) <- colnames(pv)[-3]
# subset the pvalues, by merging the CN list
pv_final <- merge(CN2, pv, by.x = c("group2", "group1"), by.y = c("group1", "group2"))
# fix ordering
pv_final <- pv_final[order(pv_final$group1), ] 
# set signif level
pv_final$map_signif <- ifelse(pv_final$p.value > 0.05, "", ifelse(pv_final$p.value > 0.01,"*", "**"))  

# the plot
ggplot(mydf, aes(x=Group, y=Value, fill=Group)) + geom_boxplot() +
  stat_compare_means(data=mydf[ mydf$Group != "PGMC4", ], aes(x=Group, y=Value, fill=Group), size=5) + 
  ylim(-4,30)+
  geom_signif(comparisons=CN,
              y_position = 3:30, annotation= pv_final$map_signif) + 
  theme_bw(base_size = 16)

enter image description here

The arguments vjust, textsize, and size are not properly working. Seems to be a bug in the latest version ggsignif_0.3.0.

Edit: When you want to show only the significant comparisons, you can easily subset the dataset CN. Since I updated to ggsignif_0.4.0 and R version 3.4.1, vjust and textsize are working now as expected. Instead of y_position you can try step_increase.

# subset 
gr <- pv_final$p.value <= 0.05
CN[gr]

ggplot(mydf, aes(x=Group, y=Value, fill=Group)) + 
  geom_boxplot() +
  stat_compare_means(data=mydf[ mydf$Group != "PGMC4", ], aes(x=Group, y=Value, fill=Group), size=5) + 
  geom_signif(comparisons=CN[gr], textsize = 12, vjust = 0.7, 
             step_increase=0.12, annotation= pv_final$map_signif[gr]) + 
  theme_bw(base_size = 16)

You can use ggpubr as well. Add:

stat_compare_means(comparisons=CN[gr], method="wilcox.test", label="p.signif", color="red")

enter image description here

Related questions

Put stars on ggplot barplots and boxplots - to indicate the level of significance (p-value)

It's common to put stars on barplots or boxplots to show the level of significance (p-value) of one or between two groups, below are several examples: The number of stars are defined by p-value, for example one can put 3 stars …

Read More
Plot multiple boxplot in one graph

I saved my data in as a .csv file with 12 columns. Columns two through 11 (labeled F1, F2, ..., F11) are features. Column one contains the label of these features either good or bad. I would like to plot a boxplot of …

Read More
Transform only one axis to log10 scale with ggplot2

I have the following problem: I would like to visualize a discrete and a continuous variable on a boxplot in which the latter has a few extreme high values. This makes the boxplot meaningless (the points and even the "body" …

Read More