如何解决R ggplot2:将 kruskal Wallis 和成对 Wilcoxon 检验添加到每个组和方面内具有多个组/子组的箱线图中
我正在尝试将 kruskal Wallis 和成对 Wilcoxon 检验添加到图中以显示哪些组有显着差异,但我在每个组和方面中有多个组/子组,这使它变得复杂。
这里是使用 iris 数据集作为示例的 R 代码,其想法是针对不同的变量(Sepal.Length、Sepal.Width、Petal.Length、 Petal.Width) 每个物种,以及它们之间成对的 wilcox.test:
rm(list=ls(all=TRUE)); cat('\014') # clear workspace
library(tidyverse)
library(ggplot2)
library(viridis)
library(rstatix)
data(iris)
iris$treatment <- rep(c("A","B","C"),length(iris$Species)/3)
mydf <- gather(iris,Variable,value,Sepal.Length:Petal.Width)
# change number to create more difference
mydf[mydf$treatment=="B",]$value <- mydf[mydf$treatment=="B",]$value*1.2
#mydf[mydf$treatment=="C",]$value <- mydf[mydf$treatment=="C",]$value+0.3
# do pairwise Wilcoxon test for pairwise comparisons between groups
df_wilcox <- mydf %>%
group_by(Species,Variable) %>%
pairwise_wilcox_test(value ~ treatment) %>%
add_y_position(step.increase = 0.02)
# do Kruskal Wallis test to see whether or not there is statistically significant difference between three or more groups
df_kw <- compare_means(value ~ treatment,mydf,group.by = c("Species","Variable"),method="kruskal")
# plot boxplot with wilcoxon and kruskal test results
P <- ggplot(data=mydf,aes(x=treatment,y=value,fill=Variable))+
stat_boxplot(geom = "errorbar")+geom_boxplot(outlier.shape = NA)+
facet_wrap(~Species,nrow=1)+
theme_bw()+
theme(axis.text=element_text(size=12),axis.title=element_text(size=16),plot.title=element_text(size=20)) +
theme(strip.text = element_text(size=14))+
scale_fill_viridis(discrete = TRUE) +
guides(fill=guide_legend(title="Variable"))+
stat_pvalue_manual(df_wilcox,color ="Variable",step.group.by="Variable",tip.length = 0,step.increase = 0.02)
#stat_pvalue_manual(df_wilcox,step.increase = 0.02,hide.ns=T) #hide non-significant
# change legend title and wilcoxon test color
ggpar(P,legend.title = "Wilcoxon test",palette = c("#440154FF","#3B528BFF","#21908CFF","#FDE725FF"))
为了改善这个数字,我想:
- 自动将来自“df_kw”的 Kruskal 测试结果作为文本添加到图中,并且只显示显着的 p 值(例如 KW(petal.length) p = 0.003)
- 使不同变量(例如花瓣/Speal 长度/宽度)的处理(例如“A”、“B”、“C”)之间的 wilcoxon 线看起来很整齐(例如全部在箱线图的顶部,线间距一致)
- 使 wilcoxon 测试线的颜色与箱线图的颜色相同(现在,如果 wilcoxon 测试变量小于实际变量,则如果我隐藏不显着的 'ggpar' 并不总是有效)
我被困在这里,想知道有人有解决方案吗?非常感谢!
解决方法
我可以回答关于如何自动将 pvalues 标签添加到图中的问题的第一部分。一种方法是将mydf
和df_kw
组合起来,以便df_kw
包含与mydf
相同的所有列。在这里,我使用 data.table
包这样做:
setDT(mydf); setDT(df_kw) # convert to data.tables by reference
df_kw <- mydf[df_kw,mult = "first",on = c("Variable","Species"),nomatch=0L] #creates data table with the same columns as mydf
df_kw <- df_kw[df_kw$p < 0.05,] #removes non-significant values
然后您可以使用 geom_text
自动添加标签。我会先生成一个值的字符向量来定位这样的标签:
y_lab_placement <- c(sort(rep(seq(max(mydf$value)*1.25,by = -0.35,length.out = length(unique(mydf$Variable))),length(unique(mydf$Species))),decreasing = T)) # creates y values of where to place the labels
y_lab_placement <- y_lab_placement[1:nrow(df_kw)] # adjusts length of placements to the length of significant values
然后我会将此行添加到您的 ggplot 以添加标签:
geom_text(data = df_kw,aes(x = 2,y = y_lab_placement,label = c(paste(Variable,"KW p ~",round(p,5)))))+ #adds labels to the plot based on your data
这是包括这些版本在内的整个代码块。
rm(list=ls(all=TRUE)); cat('\014') # clear workspace
library(tidyverse)
library(ggplot2)
library(viridis)
library(rstatix)
library(data.table) # used in creating combined data table
data(iris)
iris$treatment <- rep(c("A","B","C"),length(iris$Species)/3)
mydf <- gather(iris,Variable,value,Sepal.Length:Petal.Width)
# change number to create more difference
mydf[mydf$treatment=="B",]$value <- mydf[mydf$treatment=="B",]$value*1.2
#mydf[mydf$treatment=="C",]$value <- mydf[mydf$treatment=="C",]$value+0.3
# do pairwise Wilcoxon test for pairwise comparisons between groups
df_wilcox <- mydf %>%
group_by(Species,Variable) %>%
pairwise_wilcox_test(value ~ treatment) %>%
add_y_position(step.increase = 0.02)
# do Kruskal Wallis test to see whether or not there is statistically significant difference between three or more groups
df_kw <- compare_means(value ~ treatment,mydf,group.by = c("Species","Variable"),method="kruskal")
setDT(mydf); setDT(df_kw) # convert to data.tables by reference
df_kw <- mydf[df_kw,] #removes non-significant values
# plot boxplot with wilcoxon and kruskal test results
y_lab_placement <- c(sort(rep(seq(max(mydf$value)*1.25,decreasing = T)) # creates y values of where to place the labels
y_lab_placement <- y_lab_placement[1:nrow(df_kw)] # adjusts length of placements to the length of significant values
P <- ggplot(data=mydf,aes(x=treatment,y=value,fill=Variable))+
stat_boxplot(geom = "errorbar")+geom_boxplot(outlier.shape = NA)+
facet_wrap(~Species,nrow=1)+
theme_bw()+
theme(axis.text=element_text(size=12),axis.title=element_text(size=16),plot.title=element_text(size=20)) +
theme(strip.text = element_text(size=14))+
scale_fill_viridis(discrete = TRUE) +
guides(fill=guide_legend(title="Variable"))+
geom_text(data = df_kw,5)))))+ #adds labels to the plot based on your data
stat_pvalue_manual(df_wilcox,color ="Variable",step.group.by="Variable",tip.length = 0,step.increase = 0.02)
#stat_pvalue_manual(df_wilcox,step.increase = 0.02,hide.ns=T) #hide non-significant
# change legend title and wilcoxon test color
ggpar(P,legend.title = "Wilcoxon test",palette = c("#440154FF","#3B528BFF","#21908CFF","#FDE725FF"))
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。