R ggplot2: boxplots with significance level (more than 2 groups: kruskal.test and wilcox.test pairwise) and multiple facets

Question

Following up on this question, I am trying to make boxplots and pairwise comparisons to show levels of significance (only for the significant pairs) again, but this time I have more than 2 groups to compare and more complicated facets.

I am going to use the iris dataset here for illustration purposes. Check the MWE below where I add an additional "treatment" variable.

library(reshape2)
library(ggplot2)
data(iris)
iris$treatment <- rep(c("A","B"), length(iris$Species)/2)
mydf <- melt(iris, measure.vars=names(iris)[1:4])
ggplot(mydf, aes(x=variable, y=value, fill=Species)) + geom_boxplot() +
stat_summary(fun.y=mean, geom="point", shape=5, size=4) +
facet_grid(treatment~Species, scales="free", space="free_x") +
theme(axis.text.x = element_text(angle=45, hjust=1)) 

This produces the following plot:

The idea would be to perform a Kruskal-Wallis test across the "variable" groups (Sepal.Length, Sepal.Width, Petal.Length, Petal.Width), and pairwise Wilcoxon tests between them, PER FACET defined by "Species" and "treatment".

It would most likely involve updating the annotation like in my previous question.

In other words, I want to do the same as in this other question I posted, but PER FACET.

I am getting horribly confused and stuck, though the solution should be quite similar... Any help would be appreciated!! Thanks!!

Godsend

Answer 1

You can try

library(ggsignif)
ggplot(mydf,aes(x=variable, y=value)) + 
  geom_boxplot(aes(fill=Species)) + # define the fill argument here
  facet_grid(treatment~Species) +
  ylim(0,15)+
  theme(axis.text.x = element_text(angle=45, hjust=1)) +
  geom_signif(test="wilcox.test", comparisons = combn(levels(mydf$variable),2, simplify = F)[-4],
              step_increase = 0.2)

Kruskal.wallis can be included by adding

library(ggpubr)
stat_compare_means(test="kruskal.test")

Urbas

Hi Jimobu, many thanks for your time. I was reaching something like that, but I have a few concerns with it... First thing is I wonder why I don't get the same p-values with geom_signif() than with pv <- tidy(with(mydf, pairwise.wilcox.test(value, interaction(treatment,Species,variable), p.adjust.method = "BH"))) – Godsend

The other thing would be to put p-values in the form of * and ** not showing those non-significant comparisons. Exactly the same we did in #45477450 but per facet... I actually need the same p-values, cause the plots we made in the previous question were subsets from the same data – Godsend

p-values: you have to compare non-adjuasted pvalues. thus use p.adjust.method = "none". Use this solution for the other problem. – Urbas

Hi @Jimbou thanks for your comment! the question is the opposite, as in how to show BH-adjusted p.values with ggsignif – Godsend

In #45553215 I am confused with the step to go from "pv" to "pv.final", since I'n not familiar with dplyr and it shows errors with this data... – Godsend

You don't have to use dplyr. Simply filter and order the pv according the head(myplot2$data[[2]]) data.frame. – Urbas

That's exactly the step I cannot figure out, cause myplot2$data[[2]] does not have "treatment" or "Species" information, so I cannot do any sort of merge or anything... – Godsend

I've posted the following question, I would need to get a list in the same order as myplot2$data[[2]]$annotation, and then merge pv to that and get pv.final to replace the annotation. #46467690 – Godsend

Guys, what if I don't have the treatment variable? How could I still get the comparisons showed? – Grochow

@Grochow what do you mean with.... "I don't have the treatment variable"? As your vairable also has stored the x-values you would not be able to create the plot. – Urbas

How can I do the same thing when I also have groups? For example, can you do the same comparisons when A and B are groups, instead of vertical facets? My code does not work for that. – Russom

@Russom this depends how the group treatment is specified in aes(). Is it provided in color or fill, then my answer is: no, not easy and very complicated since you have to hack in the underlying grids. When the group is added using i.e., interaction(treatment, variable) than it is easily possible. – Urbas

@Urbas are you saying that I need to use groups and not fill? Because grouping with 'fill' hasn't worked so far. – Russom

@Russom No, i want to say that fill = treatment works for grouping along the x-axis, but works not or only in a very difficult way for the p-values. – Urbas