R - stat_compare_means 从 Kruskal-Wallis 测试返回不同的值

Posted 2023-02-16

【中文标题】R - stat_compare_means 从 Kruskal-Wallis 测试返回不同的值

【英文标题】：R - stat_compare_means return differnt value from Kruskal-Wallis test

【发布时间】：2020-04-17 00:38:55

【问题描述】：

我想使用 ggpubr 包中的 R 函数 stat_compare_means 将 Kruskal-Wallis 测试的 p 值绘制到我的 ggplot。

但是，如果我只是运行该函数，则绘制的值与该值不同：

kruskal.test(value ~ type, data = Profile_melt)

我绘制 p 值的代码是：

ggplot(Profile_melt, aes(type, value)) + 
  geom_boxplot(aes(fill = factor(type), alpha = 0.5), 
               outlier.shape = NA, show.legend = FALSE) +
  geom_jitter(width = 0.2, size = 2, show.legend = FALSE,
              aes(colour = factor(type)), alpha = 0.5) +
  theme_bw() +
  facet_grid(Case ~ Marker, scales = 'free') +
  stat_compare_means(comparison = list(c("Real", "Binomial")),method = 'kruskal.test')+
  background_grid(major = 'y', minor = "none") + # add thin horizontal lines 
  xlab('Category') +
  ylab('Cell counts (Frequencies)')+
  theme(axis.text = element_text(size = 15), 
        axis.title = element_text(size = 20), 
        legend.text = element_text(size = 38),
        legend.title = element_text(size = 30), 
        strip.background = element_rect(colour="black", fill="white"),
        strip.text = element_text(margin = margin(10, 10, 10, 10), size = 25)) +
  panel_border()

这是我的数据sample data

【参考方案1】：

有许多代码行可能与问题无关。也许，您的问题可能是：

为什么

kruskal.test(value ~ type, data = Profile_melt)

#Kruskal-Wallis chi-squared = 4.9673, df = 1, p-value = 0.02583

产生不同的 p 值

ggboxplot(Profile_melt, x="type", y = "value") + 
  stat_compare_means(comparison = list(c("Real", "Binomial")), method = 'kruskal.test')

# p-value = 0.49

您可以通过检查原始代码找出原因。 ggpubr 的开发人员可能会更好地解释这一点，如果有问题，可能会在那里修复它。要获得正确且一致的 p 值，请删除 comparison = list(c("Real", "Binomial"))：

ggboxplot(Profile_melt, x="type", y = "value") + 
  stat_compare_means(method = 'kruskal.test')

或

编辑

ggboxplot(Profile_melt, x="type", y = "value") + 
  stat_compare_means(comparison = list(c("Real", "Binomial")))

使用您的其他代码，图表如下所示：

【讨论】：

您好，志强，谢谢您的回复。删除这个比较 = ... 行确实会产生一致的结果，但是，绘图的格式也会改变，这不是我想要的。

我明白了。你想要漂亮的水平线。您可以通过删除 method = 'kruskal.test' 而不是 comparison = list(c("Real", "Binomial") 来实现。我会修改我的答案。

请运行我提供的完整代码。您的建议在方面不起作用。

我已经在我这边运行了代码。它与facet 的工作方式相同，因为您的样本很小，它会产生一些警告cannot compute exact p-value with ties。正如我所说，真正的解决办法可能是更改ggpubr。

嗯，这很奇怪...您能否提供您身边的完整代码（带方面）？谢谢

【参考方案2】：

stat_compare_means 来自 ggpubr 调用 compare_means 默认使用 wilcox.test。因此，正如@ZhiqiangWang 指出的那样，如果您删除该方法或比较，它将进入默认值，这与您首先获得的 p 值相似，因为 2 个样本的 wilcoxon 和 kruskal 非常相似：

kruskal.test(value ~ type, data = Profile_melt)
#Kruskal-Wallis chi-squared = 4.9673, df = 1, p-value = 0.02583
wilcox.test(value ~ type, data = Profile_melt)
#W = 1034939, p-value = 0.02583

现在，对于您拥有的数据，您很可能需要每个单独的案例和标记的 p 值，而不是使用 kruskal.test(value ~ type, data = Profile_melt) 的泛比较。为所有方面打印相同的 p 值是没有意义的。

我们首先检查我们需要的 p 值：

compare_means(value ~ type, Profile_melt, group.by = c("Case","Marker"),
method="kruskal")
# A tibble: 30 x 8
   Case    Marker .y.            p   p.adj p.format p.signif method        
   <fct>   <fct>  <chr>      <dbl>   <dbl> <chr>    <chr>    <chr>         
 1 Case 1A CD3    value 0.000470   0.0085  0.00047  ***      Kruskal-Wallis
 2 Case 1A CD4    value 0.00000915 0.00022 9.2e-06  ****     Kruskal-Wallis
 3 Case 1A CD8    value 0.00695    0.09    0.00695  **       Kruskal-Wallis
 4 Case 1A CD20   value 0.707      1       0.70724  ns       Kruskal-Wallis
 5 Case 1A FoxP3  value 0.00102    0.014   0.00102  **       Kruskal-Wallis
 6 Case 1B CD3    value 0.0000415  0.00091 4.1e-05  ****     Kruskal-Wallis

类似于：

Profile_melt %>% 
group_by(Case,Marker) %>% 
summarize(k_p=kruskal.test(value ~ type)$p.value)

# A tibble: 30 x 3
# Groups:   Case [6]
   Case    Marker        k_p
   <fct>   <fct>       <dbl>
 1 Case 1A CD3    0.000470  
 2 Case 1A CD4    0.00000915
 3 Case 1A CD8    0.00695   
 4 Case 1A CD20   0.707     
 5 Case 1A FoxP3  0.00102   

我们可以绘图，使用 ggpubr 包中的 ggboxplot 一定更容易：

p = ggboxplot(Profile_melt,x="type",y="value",add="jitter",
facet.by=c("Case","Marker"),scales="free_y",ggtheme=theme_pubclean())

p+stat_compare_means(
aes(label =paste("p=",scientific(as.numeric(..p.format..)))),
method="kruskal",size=2)

【讨论】：

非常感谢！这正是我想要的。还有一件事，我想要每个方面的水平比较线，我想将数字格式统一为科学记数法（1.23e-3 等）。我怎么能这样做？

我猜你的意思是网格线，因为你必须使用 ggboxplot 中的 ggtheme 选项选择一个主题，或者你使用 + theme(..) 设置它。至于科学，你可以看到我上面的内容，不幸的是 ggpubr 已经有 1 个有效数字，所以我对此无能为力

这有点超出您的问题范围，这是关于 p 值的。因此，如果您需要对情节进行更多修改，我建议您将其作为一个单独的问题发布..