logo标签列表

R ggplot2箱线图 - ggpubr的stat_compare_means函数不能正常工作

rggplot2boxplotp-valueggpubr
7
我正在尝试使用ggplot2和ggpubr包，以星号的形式向我的箱线图添加显著性水平，但我有许多比较，只想显示显著的那些。
我尝试在stat_compare_means中使用hide.ns=TRUE选项，但它明显不起作用，这可能是ggpubr包中的一个错误。
此外，您可以看到我从成对的wilcox.test比较中省略了组“PGMC4”；我如何在kruskal.test中也省略该组？
我最后一个问题是显著性水平是如何工作的？例如，*表示小于0.05的显着性，**表示小于0.025的显着性，***表示小于0.01的显着性？ ggpubr使用的是什么约定？它显示p值还是调整后的p值？如果是后者，调整方法是什么？BH？
请查看下面的MWE以及参考此链接和这个链接。
##############################
##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()
##############################

这个 MWE 将会生成以下的箱线图：

test

问题如下：

1- 如何使hide.ns=TRUE起作用？

2- 如何增加*的大小？

3- 如何从kruskal.test比较中排除一组？

4- ggpubr使用的*约定是什么，显示的p值是否已经调整？

非常感谢！

编辑

另外，在进行时

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

我得到的p值与进行时不同

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

其中，mydf.sub是给定两组比较的mydf的子集()。

ggpubr在这里做什么？它如何计算p值？

编辑2

请帮助我，解决方案不一定要使用ggpubr（但必须使用 ggplot2 ），我只需要能够隐藏NS并使星号的大小更大，并且进行与wilcox.test（）+ p.adjust（method"BH"）相同的p值计算。

谢谢！

  1. 是的确实。看起来是一个bug。
  2. 没有头绪。
  3. 使用stat_compare_means(data=mydf[ mydf$Group != "PGMC4", ],aes(x=Group, y=Value, fill=Group), size=5)。
  4. 将结果与pairwise.wilcox.test(mydf$Value, mydf$Group, p.adjust.method = "none")进行比较。
 - Roman
自己绘制所有内容，例如使用 library(ggsignif);geom_signif() 和注释。请参见此处的最后一个答案：https://dev59.com/Y2Qn5IYBdhLWcg3wIUFn#27073333 - Roman
你能否将这个问题进一步发展成一个答案？使用ggsignif的geom_signif()函数，我似乎无法弄清如何去除NS比较、如何改变**的大小以及如何放置括号，使它们不会重叠（就像ggpubr一样）...所以我还是停留在同一个点上。 - DaniCee
4- ggpubr使用的*约定是什么：人们通常喜欢根据p值的数量级来解释显著性水平，例如<0.05 = *，<0.001 = **，<0.0001 ***。这是一个完全错误的观点。如果p值小于alpha值（通常为0.05），则p值是显著的，否则不显著。所有这些“显著性水平”都是无意义的。（NHST本身就是一种巫术，但这是另一个讨论话题） - Scransom
如果有人能找到一种去除NS并增加*大小的方法，我将感到非常满意...谢谢！ - DaniCee
显示剩余3条评论
1个回答
11

编辑：自从我发现了rstatix包之后，我会这样做：

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

你可以尝试以下方法。思路是使用pairwise.wilcox.test自己计算统计数据。然后使用ggsignif函数中的geom_signif来添加预先计算好的p值。使用y_position可以放置括号，以避免重叠。
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

参数vjust、textsize和size无法正常工作。最新版本ggsignif_0.3.0似乎存在一个错误。

编辑：当您只想显示重要比较时，可以轻松地对数据集CN进行子集处理。自从我更新到ggsignif_0.4.0和R version 3.4.1以来，vjust和textsize现在都按预期工作。您可以尝试使用step_increase而不是y_position。
# 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)

你也可以使用ggpubr。添加：

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

enter image description here

我正在使用ggsignif_0.4.0中的geom_signif，而不是ggpubr...无法更改*大小... - DaniCee
哦，我明白了！即使您在实际绘图中没有使用它，加载ggpubr也会使其出现问题...让我消化一下这个答案并接受它。 - DaniCee
@DaniCee 目前我在 ggsignif 中没有找到解决方案。你可以尝试使用 geom_text(x=2, y=6, label=paste("Kruskal-Wallis, p=", round(with(mydf[ mydf$Group != "PGMC4", ], kruskal.test(Value~ Group)$p.value),2)))。 - Roman
你能否看一下我刚刚发布的这个问题（https://dev59.com/a6Tia4cB1Zd3GeqP8yHw）？我感觉它肯定和某些问题非常相似，但是我无法解决它... 非常感谢！ - DaniCee
嗨 @Jimbou，我正在尝试做类似的事情，但这次我要同时比较多个组和多个因素...如果你能看一下，我会很感激！谢谢 https://stackoverflow.com/questions/46446392/r-ggplot2-boxplots-with-significance-level-more-than-2-groups-kruskal-test-an - DaniCee
显示剩余4条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接