使用stat_signif和自定义测试ggplot2

问题描述

我的问题与ggsignif package error stat_signif requires the following missing aesthetics: y这里发布的问题非常相似,除了我有多个小组,因此除了代表总体卡方检验之外,我还需要跟进每个小组之间的卡方比较。据我所知,stat_signif没有用于后验卡方的参数,但是它确实接受自定义公式,只要它具有一个包含名为p值的条目的列表即可。

我正在使用不推荐使用的fifer软件包中的代码进行事后卡方检验,我还编辑了公式,因此输出中带有一列标记为p.value的列(请参阅本文底部的已编辑公式)。

下面是我的数据示例(虚构)

require(tidyverse)
require(ggplot2)
require(ggsignif)
require(ggpubr)


condition <- c("a","a","b","c","c")
binary_1 <- c(1,1,0)

df <- data.frame(condition,binary_1)
df

gg_df <- df %>%
  mutate(binary_1 = as.factor(binary_1))

gg_melt <- melt(gg_df)


gg_1.1 <- gg_melt %>%
  group_by(condition) %>%
  count(binary_1) %>%
  mutate(Freq = n) %>%
  ggplot() +
  aes(condition,Freq,fill = binary_1) +
  geom_col(position = "stack") +
  geom_signif(
    comparisons = list(c("a","b"),c("a","c"),c("b","c")),test = "chisq.post.hoc",map_signif_level = TRUE
  ) +
  scale_fill_manual(values = c("#FDAE61","#9E0142"),name = "Binary_1")

gg_1.1

但是,这会产生以下错误

Warning message:
computation Failed in `stat_signif()`:
'c(4,2)' is not a function,character or symbol

如果我使用相同的语法但未指定测试,则它可以工作,但是会自动执行wilcox测试,并且出现以下错误

Warning message:
In wilcox.test.default(c(4,1),c(1,4)) :
  cannot compute exact p-value with ties
gg_1.2 <- gg_melt %>%
  group_by(condition) %>%
  count(binary_1) %>%
  mutate(Freq = n) %>%
  ggplot() +
  aes(condition,name = "Binary_1")

gg_1.2

下面是我用来计算事后卡方检验的函数

chisq.post.hoc <- function(tbl,test=c("fisher.test"),popsInRows=TRUE,control=c("fdr","BH","BY","bonferroni","holm","hochberg","hommel"),digits=4,...) {
  #### extract correction method
  control <- match.arg(control)
  
  #### extract which test (fisher or chi square) 
  test = match.fun(test)
  
  #### test rows or columns
  if (!popsInRows) tbl <- t(tbl)
  popsNames <- rownames(tbl)
  
  #### come up with all possible comparisons
  prs <- combn(1:nrow(tbl),2)
  
  #### preallocate  
  tests <- ncol(prs)
  pvals <- numeric(tests)
  lbls <- character(tests)
  for (i in 1:tests) {
    pvals[i] <- test(tbl[prs[,i],],...)$p.value
    lbls[i] <- paste(popsNames[prs[,i]],collapse=" vs. ")
  }
  adj.pvals <- p.adjust(pvals,method=control)
  cat("Adjusted p-values used the",control,"method.\n\n")
  data.frame(comparison=lbls,raw.p=round(pvals,digits),p.value=round(adj.pvals,digits))
}

这是获取p值的输出

chisq_matrix <- as.matrix(table(df$condition,df$binary_1))
chisq.test(chisq_matrix)

chisq.post.hoc(chisq_matrix,test='chisq.test')

我尝试将输出保存为列表并重新运行代码,但出现错误

Error in get(as.character(FUN),mode = "function",envir = envir) : 
  object 'chisq_list' of mode 'function' was not found

chisq_list <- chisq.post.hoc(chisq_matrix,test='chisq.test')

gg_1.3 <- gg_melt %>%
  group_by(condition) %>%
  count(binary_1) %>%
  mutate(Freq = n) %>%
  ggplot() +
  aes(condition,test = "chisq_list",name = "Binary_1")

如何保存事后测试中的p值,以使其显示在图形上?

解决方法

暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!

如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@)