Как применить p-значение для каждой группы данных в R, используя facet_wrap в ggpubr

Question

У меня есть данные, которые выглядят так:

melted.df <- structure(list(Time = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 
4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L), .Label = c("24", 
"36", "48", "72"), class = "factor"), id = c(1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 20L, 21L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 
12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 1L, 2L, 3L, 
4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 
18L, 19L, 20L, 21L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L), Samples = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L), .Label = c("WT_Ago2_800", "WT_Ago2_400", "WT_Ago2_200", 
"WT_Ago4_800"), class = "factor"), Size = c(0, 0, 0, 0, 0, 0, 
0.3, 0, 0, 0.1, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.8, 0.5, 0, 0, 
0, 0, 0, 0, 0.1, 0.65, 0.2, 0.85, 0.725, 0.575, 0.1, 1.1, 0.9, 
1.325, 1, 0.8, 0.5, 2.2, 1.65, 0, 0, 0, 0, 0, 0, 0.825, 1.175, 
0.1, 0.55, 0.85, 0.85, 1.1, 1.4, 0.6, 0.95, 1.15, 0.975, 2.35, 
1.15, 2.1, 0, 0, 0, 0, 0, 0, 0.65, 1.4, 0.55, 0.1, 0.7, 1.1, 
0.95, 1.85, 0.85, 0.1, 1.5, 1.25, 1.8, 1.75, 2.15)), row.names = c(NA, 
-84L), class = "data.frame")

Эти данные состоят из 4 временных периодов (24, 36, 48 и 72 часа). Я хочу использовать приведенный ниже код, чтобы вставить значения p, рассчитанные как stat.test для каждого time.levels, и применить его к каждому facet_wrap. Если вы проверите i = 1, p-значения не будет, поэтому вы ничего не захотите применить к фигуре, а если вы сделаете i = 2, вы получите p-значения, примененные к фигуре. Проблема в том, что мне не удалось применить значение p к соответствующим фасетам. Он просто применяет одно и то же значение p во всех аспектах. Как я могу решить эту проблему?

код:

library(devtools)
# install_github("https://github.com/kassambara/rstatix")
library(rstatix)  # https://github.com/kassambara/rstatix
library(stringi)
library(ggpubr)
time.levels <- levels(melted.df$Time)
stat.test <- NULL
for (i in 1:length(time.levels)){
  stat.test <- aov(Size ~ Samples, data = melted.df[melted.df$Time == time.levels[i],]) %>%
  tukey_hsd()
# stat.test <- rbind(stat.test, tmp.stat)
 bp <-  ggboxplot(melted.df, x = "Samples", y = "Size") +
    facet_wrap(vars(Time))+
    stat_pvalue_manual(
      stat.test, label = "p.adj",
      y.position = c(2, 2.5, 3, 3.5, 3.8, 4)
    )
 bp

}

Answer 1

Вам не нужно использовать gridExtra::grid.arrange.

Вот чистое решение.

library(rstatix) # latest version
library(ggpubr)  # latest version
stat.test <- melted.df %>%
  group_by(Time) %>%
  tukey_hsd(Size ~ Samples)

ggboxplot(melted.df, x = "Samples", y = "Size", facet.by = "Time") +
  stat_pvalue_manual(
    stat.test, label = "p.adj",
    y.position = c(2, 2.5, 3, 3.5, 3.8, 4)
  )

Answer 2

Примечание.Все ваши значения в размере для Time == 24L равны нулю:

> filter(melted.df, Time == 24L) %>% select(Size) %>% summary
      Size  
 Min.   :0  
 1st Qu.:0  
 Median :0  
 Mean   :0  
 3rd Qu.:0  
 Max.   :0  

Если вы все равно хотите продолжить, вы должны составить графики по отдельности и затем использовать gridExtra::grid.arrange:

library(gridExtra)
bp <- vector("list", length = length(time.levels))
for (i in seq_along(time.levels)) {
  sdf <- melted.df[melted.df$Time == time.levels[i],]
  stat.test <- aov(Size ~ Samples, data = sdf) %>%
      tukey_hsd()

 bp[[i]] <-  ggboxplot(sdf, x = "Samples", y = "Size") +
    facet_wrap(vars(Time))+
    stat_pvalue_manual(
      stat.test, label = "p.adj",
      y.position = c(2, 2.5, 3, 3.5, 3.8, 4)
    )
}
do.call(grid.arrange, bp)

Обратите внимание , что вы должны использовать подмножество data.frame sdf в качестве входа для ggboxplot.