R ggplot2：具有wilcoxon显着性水平和方面的箱子图。只显示与星号的重要比较 [英] R ggplot2: boxplots with wilcoxon significance levels, and facets. Show only significant comparisons with asterisks
本文介绍了R ggplot2：具有wilcoxon显着性水平和方面的箱子图。只显示与星号的重要比较的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
为了完整起见，我修改了接受的答案并定制了结果图，但我仍然面临一些重要问题。
$总而言之，我正在做箱形图反映Kruskal-Wallis的意义和成对Wilcoxon测试比较。
我想替换p值数字用星号表示，只显示重要的比较结果，将垂直间距减小到最大值。
基本上我想要做
现在颜色问题变得更加突出，饰面高度不均匀，也可以使用多余的小平面文字。
我被困在这一点，所以希望有任何帮助。很抱歉，这个问题很长，但我认为它已经差不多了！谢谢!!
解决方案
您可以尝试以下操作。由于你的代码真的很忙，对我来说太复杂了，我建议采用不同的方法。我试图避免循环，并尽可能地使用 tidyverse 。因此，首先我创建了你的数据。然后计算kruskal wallis测试，因为这在 ggsignif 中是不可能的。之后，我将使用 geom_signif 来绘制所有p.values。最后，微不足道的将被删除，并增加一个步骤。
1-使着色工作完成
2-显示星号而不是数字完成
...以及获胜：
<3>制作一个常见的图例完成
<4>将Kruskal-Wallis线放在最上面完成后，我将值放在底部
5-更改标题和y轴文本的大小（和对齐方式）完成
library（tidyverse）
library（ggsignif）
＃1。您的数据
set.seed（2）
df< - as.tbl（iris）％>％
mutate（treatment = rep（c（A，B）） （key，value，-Species，-treatment）％>％
mutate（value = rnorm（n（））） ％>％
mutate（key = factor（key，levels = unique（key）））％>％
mutate（both = interaction（treatment，key，sep =））
＃2.克鲁斯卡尔测试
KW < - df％>％
group_by（种类）％>％
汇总（p = round（kruskal.test（value〜both）$ p.value，2），
y = $（$）= min（value），
x = 1）％>％
mutate（y = min（y））
＃3.地块
P < df％>％
ggplot（aes（x = both，y = value））+
geom_boxplot（aes（fill =物种））+
facet_grid（〜物种）+
ylim（-3,7）+
theme（axis.text.x = element_text（angle = 45，hjust = 1））+
geom_signif（comparisons = combn（levels（df $ both）， 2，简化= F），
map_signif_level = T）+
stat_summary（fun.y = mean，geom =point，shape = 5，size = 4）+
xlab ）+
geom_text（data = KW，aes（x，y = y，label = paste0（KW p =，p）），hjust = 0）+
ggtitle（Plot） + ylab（这是我自己的y实验室）
＃4.删除不重要的值并添加步骤增加
P_new < - ggplot_build（P）
P_new $数据[[2]] < - P_new $ data [[2]]％>％
过滤器（注释！=NS。）％>％
group_by（PANEL）％> ％
mutate（ index =（as.numeric（group [drop = T]） - 1）* 0.5）％>％
mutate（y = y + index，
yend = yend + index）％>％
select（-index）％>％
as.data.frame（）
＃最终情节
plot（ggplot_gtable（P_new））
和类似的方法使用两个方面
＃--------------------
＃5. Kruskal
KW <-df％>％
group_by（物种，处理）％>％
汇总（p = round（kruskal.test（value 〜）$ p.value，2），
y = min（value），
x = 1）％>％
ungroup（）％>％
mutate（y = min（y））
＃6.有两个方面的情节
P < - df％>％
ggplot（aes（x = key， y = value））+
geom_boxplot（aes（fill =物种））+
facet_grid（处理〜物种）+
ylim（-5,7）+
他们e（axis.text.x = element_text（angle = 45，hjust = 1））+
geom_signif（比较= combn（水平（df $ key），2，简化= F），
map_signif_level = T）+
stat_summary（fun.y = mean，geom =point，shape = 5，size = 4）+
xlab（）+
geom_text（data = KW， aes（x，y = y，label = paste0（KW p =，p）），hjust = 0）+
ggtitle（Plot）+ ylab（这是我自己的y实验室）
＃7.删除不重要的值并添加步骤增加
P_new < - ggplot_build（P）
P_new $ data [[2]] < - P_new $ data [ [％2]]％>％
filter（annotation！=NS。）％>％
group_by（PANEL）％>％
mutate（index =（as.numeric （group [drop = T]） - 1）* 0.5）％>％
mutate（y = y + index，
yend = yend + index）％>％
select指数）％>％
as.data.frame（）
＃最终情节
plot（ggplot_gtable（P_new））
编辑。
关于您的 p.adjust 需求，您可以自行设置一个函数并直接在函数内调用它 geom_signif（）。
wilcox.test.BH.adjusted< - 函数（x，y，n）{
tmp< -wilcox.test（x，y）
tmp $ p.value< - p.adjust（tmp $ p.value，n = n，method =BH）
tmp
}
geom_signif（comparisons = combn（levels（df $ both），2，simplified = F），
map_signif_level = T，test =wilcox.test.BH.adjusted，
test.args = list（n = 8））
面临的挑战是要知道最终会有多少独立测试。然后你可以自己设置 n 。在这里我使用了 8 。但这可能是错误的。
Following up on this question and for the sake of completeness, I modified the accepted answer and customized the resulting plot, but I am still facing some important problems.
To sum up, I am doing boxplots reflecting significance of Kruskal-Wallis and pairwise Wilcoxon test comparisons.
I want to replace the p-value numbers with asterisks, and show only the significant comparisons, reducing vertical spacing to the max.
Basically I want to do this, but with the added problem of facets, that messes everything up.
So far I have worked on a very decent MWE, but it still shows problems...
library(reshape2)
library(ggplot2)
library(gridExtra)
library(tidyverse)
library(data.table)
library(ggsignif)
library(RColorBrewer)
data(iris)
iris$treatment <- rep(c("A","B"), length(iris$Species)/2)
mydf <- melt(iris, measure.vars=names(iris)[1:4])
mydf$treatment <- as.factor(mydf$treatment)
mydf$variable <- factor(mydf$variable, levels=sort(levels(mydf$variable)))
mydf$both <- factor(paste(mydf$treatment, mydf$variable), levels=(unique(paste(mydf$treatment, mydf$variable))))
# Change data to reduce number of statistically significant differences
set.seed(2)
mydf <- mydf %>% mutate(value=rnorm(nrow(mydf)))
##
##FIRST TEST BOTH
#Kruskal-Wallis
addkw <- as.data.frame(mydf %>% group_by(Species) %>%
summarize(p.value = kruskal.test(value ~ both)$p.value))
#addkw$p.adjust <- p.adjust(addkw$p.value, "BH")
a <- combn(levels(mydf$both), 2, simplify = FALSE)
#new p.values
pv.final <- data.frame()
for (gr in unique(mydf$Species)){
for (i in 1:length(a)){
tis <- a[[i]] #variable pair to test
as <- subset(mydf, Species==gr & both %in% tis)
pv <- wilcox.test(value ~ both, data=as)$p.value
ddd <- data.table(as)
asm <- as.data.frame(ddd[, list(value=mean(value)), by=list(both=both)])
asm2 <- dcast(asm, .~both, value.var="value")[,-1]
pf <- data.frame(group1=paste(tis[1], gr), group2=paste(tis[2], gr), mean.group1=asm2[,1], mean.group2=asm2[,2], FC.1over2=asm2[,1]/asm2[,2], p.value=pv)
pv.final <- rbind(pv.final, pf)
}
}
#pv.final$p.adjust <- p.adjust(pv.final$p.value, method="BH")
pv.final$map.signif <- ifelse(pv.final$p.value > 0.05, "", ifelse(pv.final$p.value > 0.01,"*", "**"))
cols <- colorRampPalette(brewer.pal(length(unique(mydf$Species)), "Set1"))
myPal <- cols(length(unique(mydf$Species)))
#Function to get a list of plots to use as "facets" with grid.arrange
plot.list=function(mydf, pv.final, addkw, a, myPal){
mylist <- list()
i <- 0
for (sp in unique(mydf$Species)){
i <- i+1
mydf0 <- subset(mydf, Species==sp)
addkw0 <- subset(addkw, Species==sp)
pv.final0 <- pv.final[grep(sp, pv.final$group1), ]
num.signif <- sum(pv.final0$p.value <= 0.05)
P <- ggplot(mydf0,aes(x=both, y=value)) +
geom_boxplot(aes(fill=Species)) +
stat_summary(fun.y=mean, geom="point", shape=5, size=4) +
facet_grid(~Species, scales="free", space="free_x") +
scale_fill_manual(values=myPal[i]) + #WHY IS COLOR IGNORED?
geom_text(data=addkw0, hjust=0, size=4.5, aes(x=0, y=round(max(mydf0$value, na.rm=TRUE)+0.5), label=paste0("KW p=",p.value))) +
geom_signif(test="wilcox.test", comparisons = a[which(pv.final0$p.value<=0.05)],#I can use "a"here
map_signif_level = F,
vjust=0,
textsize=4,
size=0.5,
step_increase = 0.05)
if (i==1){
P <- P + theme(legend.position="none",
axis.text.x=element_text(size=20, angle=90, hjust=1),
axis.text.y=element_text(size=20),
axis.title=element_blank(),
strip.text.x=element_text(size=20,face="bold"),
strip.text.y=element_text(size=20,face="bold"))
} else{
P <- P + theme(legend.position="none",
axis.text.x=element_text(size=20, angle=90, hjust=1),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
axis.title=element_blank(),
strip.text.x=element_text(size=20,face="bold"),
strip.text.y=element_text(size=20,face="bold"))
}
#WHY USING THE CODE BELOW TO CHANGE NUMBERS TO ASTERISKS I GET ERRORS?
#P2 <- ggplot_build(P)
#P2$data[[3]]$annotation <- rep(subset(pv.final0, p.value<=0.05)$map.signif, each=3)
#P <- plot(ggplot_gtable(P2))
mylist[[sp]] <- list(num.signif, P)
}
return(mylist)
}
p.list <- plot.list(mydf, pv.final, addkw, a, myPal)
y.rng <- range(mydf$value)
# Get the highest number of significant p-values across all three "facets"
height.factor <- 0.3
max.signif <- max(sapply(p.list, function(x) x[[1]]))
# Lay out the three plots as facets (one for each Species), but adjust so that y-range is same for each facet. Top of y-range is adjusted using max_signif.
png(filename="test.png", height=800, width=1200)
grid.arrange(grobs=lapply(p.list, function(x) x[[2]] +
scale_y_continuous(limits=c(y.rng[1], y.rng[2] + height.factor*max.signif))),
ncol=length(unique(mydf$Species)), top="Random title", left="Value") #HOW TO CHANGE THE SIZE OF THE TITLE AND THE Y AXIS TEXT?
#HOW TO ADD A COMMON LEGEND?
dev.off()
It produces the following plot:
As you can see there are some problems, most obviously:
1- Coloring does not work for some reason
2- I do not seem to be able to change the annotation with the asterisks
I want something more like this (mockup):
So we need to:
1- Make coloring work
2- Show asterisks instead of numbers
...and for the win:
3- Make a common legend
4- Place Kruskal-Wallis line on top
5- Change the size (and alignment) of the title and y axis text
IMPORTANT NOTES
I would appreciate my code is left as intact as possible even if it isn't the prettiest, cause I still have to make use of intermediate objects like "CNb" or "pv.final".
The solution should be easily transferable to other cases; please consider testing "variable" alone, instead of "both"... In this case we have 6 "facets" (vertically and horizontally) and everything gets even more screwed up...
I made this other MWE:
##NOW TEST MEASURE, TO GET VERTICAL AND HORIZONTAL FACETS
addkw <- as.data.frame(mydf %>% group_by(treatment, Species) %>%
summarize(p.value = kruskal.test(value ~ variable)$p.value))
#addkw$p.adjust <- p.adjust(addkw$p.value, "BH")
a <- combn(levels(mydf$variable), 2, simplify = FALSE)
#new p.values
pv.final <- data.frame()
for (tr in levels(mydf$treatment)){
for (gr in levels(mydf$Species)){
for (i in 1:length(a)){
tis <- a[[i]] #variable pair to test
as <- subset(mydf, treatment==tr & Species==gr & variable %in% tis)
pv <- wilcox.test(value ~ variable, data=as)$p.value
ddd <- data.table(as)
asm <- as.data.frame(ddd[, list(value=mean(value, na.rm=T)), by=list(variable=variable)])
asm2 <- dcast(asm, .~variable, value.var="value")[,-1]
pf <- data.frame(group1=paste(tis[1], gr, tr), group2=paste(tis[2], gr, tr), mean.group1=asm2[,1], mean.group2=asm2[,2], FC.1over2=asm2[,1]/asm2[,2], p.value=pv)
pv.final <- rbind(pv.final, pf)
}
}
}
#pv.final$p.adjust <- p.adjust(pv.final$p.value, method="BH")
# set signif level
pv.final$map.signif <- ifelse(pv.final$p.value > 0.05, "", ifelse(pv.final$p.value > 0.01,"*", "**"))
plot.list2=function(mydf, pv.final, addkw, a, myPal){
mylist <- list()
i <- 0
for (sp in unique(mydf$Species)){
for (tr in unique(mydf$treatment)){
i <- i+1
mydf0 <- subset(mydf, Species==sp & treatment==tr)
addkw0 <- subset(addkw, Species==sp & treatment==tr)
pv.final0 <- pv.final[grep(paste(sp,tr), pv.final$group1), ]
num.signif <- sum(pv.final0$p.value <= 0.05)
P <- ggplot(mydf0,aes(x=variable, y=value)) +
geom_boxplot(aes(fill=Species)) +
stat_summary(fun.y=mean, geom="point", shape=5, size=4) +
facet_grid(treatment~Species, scales="free", space="free_x") +
scale_fill_manual(values=myPal[i]) + #WHY IS COLOR IGNORED?
geom_text(data=addkw0, hjust=0, size=4.5, aes(x=0, y=round(max(mydf0$value, na.rm=TRUE)+0.5), label=paste0("KW p=",p.value))) +
geom_signif(test="wilcox.test", comparisons = a[which(pv.final0$p.value<=0.05)],#I can use "a"here
map_signif_level = F,
vjust=0,
textsize=4,
size=0.5,
step_increase = 0.05)
if (i==1){
P <- P + theme(legend.position="none",
axis.text.x=element_blank(),
axis.text.y=element_text(size=20),
axis.title=element_blank(),
axis.ticks.x=element_blank(),
strip.text.x=element_text(size=20,face="bold"),
strip.text.y=element_text(size=20,face="bold"))
}
if (i==4){
P <- P + theme(legend.position="none",
axis.text.x=element_text(size=20, angle=90, hjust=1),
axis.text.y=element_text(size=20),
axis.title=element_blank(),
strip.text.x=element_text(size=20,face="bold"),
strip.text.y=element_text(size=20,face="bold"))
}
if ((i==2)|(i==3)){
P <- P + theme(legend.position="none",
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.title=element_blank(),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
strip.text.x=element_text(size=20,face="bold"),
strip.text.y=element_text(size=20,face="bold"))
}
if ((i==5)|(i==6)){
P <- P + theme(legend.position="none",
axis.text.x=element_text(size=20, angle=90, hjust=1),
axis.text.y=element_blank(),
#axis.ticks.y=element_blank(), #WHY SPECIFYING THIS GIVES ERROR?
axis.title=element_blank(),
axis.ticks.y=element_blank(),
strip.text.x=element_text(size=20,face="bold"),
strip.text.y=element_text(size=20,face="bold"))
}
#WHY USING THE CODE BELOW TO CHANGE NUMBERS TO ASTERISKS I GET ERRORS?
#P2 <- ggplot_build(P)
#P2$data[[3]]$annotation <- rep(subset(pv.final0, p.value<=0.05)$map.signif, each=3)
#P <- plot(ggplot_gtable(P2))
sptr <- paste(sp,tr)
mylist[[sptr]] <- list(num.signif, P)
}
}
return(mylist)
}
p.list2 <- plot.list2(mydf, pv.final, addkw, a, myPal)
y.rng <- range(mydf$value)
# Get the highest number of significant p-values across all three "facets"
height.factor <- 0.5
max.signif <- max(sapply(p.list2, function(x) x[[1]]))
# Lay out the three plots as facets (one for each Species), but adjust so that y-range is same for each facet. Top of y-range is adjusted using max_signif.
png(filename="test2.png", height=800, width=1200)
grid.arrange(grobs=lapply(p.list2, function(x) x[[2]] +
scale_y_continuous(limits=c(y.rng[1], y.rng[2] + height.factor*max.signif))),
ncol=length(unique(mydf$Species)), top="Random title", left="Value") #HOW TO CHANGE THE SIZE OF THE TITLE AND THE Y AXIS TEXT?
#HOW TO ADD A COMMON LEGEND?
dev.off()
That produces the following plot:
Now the color problem becomes more striking, the facet heights are uneven, and something should be done with the redundant facet strip texts too.
I am stuck at this point, so would appreciate any help. Sorry for the long question, but I think it is almost there! Thanks!!
解决方案
You can try following. As your code is really busy and for me too complicated to understand, I suggest a different approach. I tried to avoid loops and to use the tidyverse as much as possible. Thus, first I created your data. Then calculated kruskal wallis tests as this was not possible within ggsignif. Afterwards I will plot all p.values using geom_signif. Finally, insignificant ones will be removed and a step increase is added.
1- Make coloring work done
2- Show asterisks instead of numbers done
...and for the win:
3- Make a common legend done
4- Place Kruskal-Wallis line on top done, I placed the values at the bottom
5- Change the size (and alignment) of the title and y axis text done
library(tidyverse)
library(ggsignif)
# 1. your data
set.seed(2)
df <- as.tbl(iris) %>%
mutate(treatment=rep(c("A","B"), length(iris$Species)/2)) %>%
gather(key, value, -Species, -treatment) %>%
mutate(value=rnorm(n())) %>%
mutate(key=factor(key, levels=unique(key))) %>%
mutate(both=interaction(treatment, key, sep = " "))
# 2. Kruskal test
KW <- df %>%
group_by(Species) %>%
summarise(p=round(kruskal.test(value ~ both)$p.value,2),
y=min(value),
x=1) %>%
mutate(y=min(y))
# 3. Plot
P <- df %>%
ggplot(aes(x=both, y=value)) +
geom_boxplot(aes(fill=Species)) +
facet_grid(~Species) +
ylim(-3,7)+
theme(axis.text.x = element_text(angle=45, hjust=1)) +
geom_signif(comparisons = combn(levels(df$both),2,simplify = F),
map_signif_level = T) +
stat_summary(fun.y=mean, geom="point", shape=5, size=4) +
xlab("") +
geom_text(data=KW,aes(x, y=y, label=paste0("KW p=",p)),hjust=0) +
ggtitle("Plot") + ylab("This is my own y-lab")
# 4. remove not significant values and add step increase
P_new <- ggplot_build(P)
P_new$data[[2]] <- P_new$data[[2]] %>%
filter(annotation != "NS.") %>%
group_by(PANEL) %>%
mutate(index=(as.numeric(group[drop=T])-1)*0.5) %>%
mutate(y=y+index,
yend=yend+index) %>%
select(-index) %>%
as.data.frame()
# the final plot
plot(ggplot_gtable(P_new))
and similar approach using two facets
# --------------------
# 5. Kruskal
KW <- df %>%
group_by(Species, treatment) %>%
summarise(p=round(kruskal.test(value ~ both)$p.value,2),
y=min(value),
x=1) %>%
ungroup() %>%
mutate(y=min(y))
# 6. Plot with two facets
P <- df %>%
ggplot(aes(x=key, y=value)) +
geom_boxplot(aes(fill=Species)) +
facet_grid(treatment~Species) +
ylim(-5,7)+
theme(axis.text.x = element_text(angle=45, hjust=1)) +
geom_signif(comparisons = combn(levels(df$key),2,simplify = F),
map_signif_level = T) +
stat_summary(fun.y=mean, geom="point", shape=5, size=4) +
xlab("") +
geom_text(data=KW,aes(x, y=y, label=paste0("KW p=",p)),hjust=0) +
ggtitle("Plot") + ylab("This is my own y-lab")
# 7. remove not significant values and add step increase
P_new <- ggplot_build(P)
P_new$data[[2]] <- P_new$data[[2]] %>%
filter(annotation != "NS.") %>%
group_by(PANEL) %>%
mutate(index=(as.numeric(group[drop=T])-1)*0.5) %>%
mutate(y=y+index,
yend=yend+index) %>%
select(-index) %>%
as.data.frame()
# the final plot
plot(ggplot_gtable(P_new))
Edit.
Regarding to your p.adjust needs, you can set up a function on your own and calling it directly within geom_signif().
wilcox.test.BH.adjusted <- function(x,y,n){
tmp <- wilcox.test(x,y)
tmp$p.value <- p.adjust(tmp$p.value, n = n,method = "BH")
tmp
}
geom_signif(comparisons = combn(levels(df$both),2,simplify = F),
map_signif_level = T, test = "wilcox.test.BH.adjusted",
test.args = list(n=8))
The challenge is to know how many independet tests you will have in the end. Then you can set the n by your own. Here I used 8. But this is maybe wrong.
这篇关于R ggplot2：具有wilcoxon显着性水平和方面的箱子图。只显示与星号的重要比较的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
