ggpubr stat_compare_means displaying inconsistent labels

Question

This is my first post on Stack Overflow. I have been searching for the solution for days without luck.

Please, consider the following code which will output the boxplot below. Unfortunately I can't post inline images because I don't have enough reputation points.

library(tidyverse)
library(reshape2)
library(ggplot2)
library(ggpubr)
library(rstatix)

df <- data.frame(group = c("g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g2", "g2", "g2", "g2", "g2", "g2", "g2", "g2", "g2", "g3", "g3", "g3", "g3", "g3", "g3", "g3", "g3", "g3", "g3", "g3", "g3", "g4", "g4", "g4", "g4", "g4", "g4", "g4", "g4", "g4", "g4", "g4", "g4", "g4"),
                 sample_1 = c(5.38, 74.43, 22.62, 51.71, 22.75, 69.08, 40.43, 26.92, 38.13, 36.22, 19.38, 9.04, 36.99, 18.15, 34.9, 73.98, 54.56, 88.25, 31.45, 8.02, 4.42, 59.8, 29.23, 104.45, 156.11, 155.33, 82.27, 302.51, 512.11, 155.99, 56.89, 391.07, 154.28, 138.87, 1466.86, 213.1, 83.42, 42.57, 76.08, 34.69, 36.53, 58.89, 14.83, 46.85, 66.38, 70.22, 106.51, 29.68),
                 sample_2 = c(0.01, 0, 0.02, 0.01, 0, 0.01, 0, 0.01, 0, 0.03, 0, 0, 0.09, 0, 0.01, 0.02, 0, 0.03, 0, 0.06, 0, 0, 0, 0.06, 0.04, 0.03, 0.06, 0.07, 0.05, 0.02, 0.02, 0.06, 0.03, 0, 0.23, 0.02, 0.01, 0.01, 0.01, 0.01, 0, 0.04, 0.01, 0.01, 0.02, 0.02, 0, 0.01),
                 sample_3 = c(0.61, 1.2, 0.8, 0.72, 0.74, 1.06, 0.46, 0.66, 0.9, 0.76, 0.63, 0.44, 4.38, 0.63, 3.77, 0.61, 0.77, 0.6, 0.49, 3.46, 0.37, 0.69, 1.18, 0.46, 0.55, 0.64, 0.99, 0.93, 5.95, 0.53, 1.01, 0.71, 0.62, 0.74, 1.61, 0.34, 0.43, 0.6, 0.48, 0.38, 0.34, 0.72, 0.54, 0.4, 0.37, 0.56, 1.93, 0.31))

long_df <- melt(df, id="group", variable.name="sample")

comparisons <- list(c("g1", "g2"), c("g2", "g3"), c("g3", "g4"), c("g1", "g3"), c("g2", "g4"), c("g1", "g4"))

long_df %>%
  ggplot(aes(x=group, y=value, fill=group)) +
  geom_boxplot() +
  geom_point(position=position_jitterdodge(0.2)) +
  stat_compare_means(aes(label = paste0("p = ", ..p.format..)), comparisons=comparisons, method="wilcox.test") +
  scale_y_continuous(trans="log10") +
  facet_wrap(~sample, scales="free")

Output from the code above:

Desired result (mockup): https://i.stack.imgur.com/I99F1.png Problem using stat_compare_means:

The p-values drawn by stat_compare_means are displayed in different formats, some are in scientific notation, and others are in decimal format with different number of decimal places. I would like to have all those p-values displayed with three decimal places (ie. 0.019) or "< 0.001" when applicable.

I have tried using the rstatix package with stat_pvalue_manual but couldn't manage to get the desired result with the following code:

stat_test <- long_df %>%
  group_by(sample) %>%
  wilcox_test(value ~ group) %>%
  mutate(p.label = ifelse(p < 0.001, "< 0.001", round(p, 3))) %>%
  add_y_position()
stat_test

long_df %>%
  ggplot(aes(x=group, y=value)) +
  geom_boxplot(aes(fill=group)) +
  stat_pvalue_manual(stat_test, label="p.label") +
  facet_wrap(~sample, scales="free") +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))

Output from the code above: https://i.stack.imgur.com/vciV5.png Problems using stat_pvalue_manual:

The pairwise comparison lines are not in the order I need like in the first image;
If I add the jitterdoge code and error comes up (Error: position_jitterdodge() requires at least one aesthetic to dodge by);
Can't have log10 y axis.

user_10415749 · Answer

You can work a bit around with it by calculating the positions yourself.

stat_test <- long_df %>%
  group_by(sample) %>%
  wilcox_test(value ~ group) %>%
  mutate(p.label = ifelse(p < 0.001, "< 0.001", round(p, 3))) %>%
  add_y_position()

m <- log10(max(long_df$value))
stat_test <- stat_test %>% mutate(y.position = rep(seq(m + m / 6, 2 * m, m / 6), 3))


long_df %>%
  ggplot(aes(x=group, y=value)) +
  geom_boxplot(aes(fill=group)) +
  stat_pvalue_manual(stat_test, label="p.label") +
  facet_wrap(~sample, scales="free") +
  scale_y_continuous(trans="log10")

What I did is taking your max value and say that you want to divide the 6 lines over the same amount of space, just make a sequence there in log10 values and repeat that 3 times (once for each group) and add them to y.position in your stat_test. Of course you can decide to give the lines less space and use any other new y.positions you like. You can do this even nicer to take the max values seperately for each group and calculate the positions for each facet instead of replicating.

m <- log10(max(long_df$value))
rep(seq(m + m / 6, 2 * m, m / 6), 3)

ggpubr stat_compare_means displaying inconsistent labels

1 Answers