2

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:

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!!

DaniCee
  • 2,397
  • 6
  • 36
  • 59

1 Answers1

3

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)

enter image description here

Kruskal.wallis can be included by adding

library(ggpubr)
stat_compare_means(test="kruskal.test")
Roman
  • 17,008
  • 3
  • 36
  • 49
  • 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")))` – DaniCee Sep 28 '17 at 03:17
  • 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 https://stackoverflow.com/questions/45476950/r-ggplot2-boxplots-ggpubr-stat-compare-means-not-working-properly 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 – DaniCee Sep 28 '17 at 03:20
  • p-values: you have to compare non-adjuasted pvalues. thus use `p.adjust.method = "none"`. Use this [solution](https://stackoverflow.com/questions/45552715/r-ggplot2-perform-pairwise-tests-per-pair-in-a-facet-and-show-the-p-values-wit) for the other problem. – Roman Sep 28 '17 at 07:34
  • Hi @Jimbou thanks for your comment! the question is the opposite, as in how to show BH-adjusted p.values with ggsignif – DaniCee Sep 28 '17 at 08:42
  • In https://stackoverflow.com/questions/45552715/r-ggplot2-perform-pairwise-tests-per-pair-in-a-facet-and-show-the-p-values-wit 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... – DaniCee Sep 28 '17 at 08:59
  • You don't have to use dplyr. Simply filter and order the `pv` according the `head(myplot2$data[[2]])` data.frame. – Roman Sep 28 '17 at 09:05
  • 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... – DaniCee Sep 28 '17 at 09:55
  • 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. https://stackoverflow.com/questions/46467190/r-make-a-list-with-the-combinations-of-one-factor-levels-for-each-combination – DaniCee Sep 28 '17 at 10:47
  • Guys, what if I don't have the treatment variable? How could I still get the comparisons showed? – guidebortoli Apr 07 '18 at 16:19
  • @guidebortoli 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. – Roman Apr 09 '18 at 08:02
  • 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. – Ben Dec 19 '19 at 10:05
  • @Ben 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. – Roman Dec 19 '19 at 10:43
  • @Roman are you saying that I need to use groups and not fill? Because grouping with 'fill' hasn't worked so far. – Ben Dec 19 '19 at 10:46
  • @Ben 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. – Roman Dec 19 '19 at 10:49