R ggplot2 - 在一个方面执行每对配对测试，并用ggsignif显示p值 [英] R ggplot2 - perform pairwise tests per pair in a facet and show the p-values with ggsignif
本文介绍了R ggplot2 - 在一个方面执行每对配对测试，并用ggsignif显示p值的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
继续阅读
我想利用ggsignif的geom_signif（）来比较DATABASE1到DATABASE2每组（每对在一个方面）并显示显着性水平（ *对于p-值<0.05，**对于p-值<0.01）进行两两比较。
任何帮助将不胜感激。谢谢!!
编辑！
我设法将显着性水平正确，但是当我添加方面的时候，一切都崩溃了......
查看我的新MWE
## MWE
库（ggplot2）
库（ggsignif）
库（tidyverse）
库（扫帚）
set.seed（1）
alpha.subA < - data.frame（Sample.ID = paste（sample（LETTERS，163，replace = TRUE）），sample（1：1000,163，replace = FALSE）， （'C'，10），代表（'FH'，10），代表（'I'，19），代表（'IF'，42）， rep（'NA'，14），rep（'NF'，42），rep（'NI'，15），rep（'NS'，10），rep（'PGMC4'，1）），
值= rnorm（n = 163，平均值= 2.3，sd = 0.45））
alpha.subA $ DB < - DATABASE1
alpha.subB < - data.frame（Sample.ID = （样本（LETTERS，163，replace = TRUE），样本（1：1000,163，replace = FALSE），sep =''），
Group = c（rep（'C'，10），rep （ 'FH'，10），代表（ 'I'，19），代表（'I F，42），代表（ 'NA'，14），代表（ 'NF'，42），代表（ 'NI'，15），代表（ 'NS'，10），代表（ 'PGMC4'，1） ），
Value = rnorm（n = 163，mean = 2，sd = 0.5））
alpha.subB $ DB < - DATABASE2
alpha.sub < - rbind alpha.subA，alpha.subB）
alpha.sub $ DB< - as.factor（alpha.sub $ DB）
alpha.sub $ both< - factor（paste（paste alpha.sub $ Group，alpha.sub $ DB），levels = paste（rep（levels（alpha.sub $ Group），each = length（levels（alpha.sub $ DB））），rep（levels（alpha.sub $ DB） $ db），长度（水平（alpha.sub $ Group）））））
v1 < - grep（DATABASE1，levels（alpha.sub $ both），val = TRUE）[ -9]
v2 < - grep（DATABASE2，levels（alpha.sub $ both），val = TRUE）[ - 9]
CNb < - mapply（c，v1，v2， SIMPLIFY = FALSE）
CNb < - unname（CNb）
pv < - tidy（with（alpha.sub [alpha.sub $ Group！=PGMC4，]，pairwise .wilcox.test（Value，both，p.adjust.method =BH）））＃ADJUSTED PVALUES
＃data preparation
CNb2 < - do.call（rbind.data.frame，CNb ）
colnames（CNb2）< - colnames（pv）[ - 3]
＃通过合并CN列表
pv.final< - merge（CNb2，pv，by.x = c（group2，group1），by.y = c（group1 ，group2））
＃修复订单
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，*，**））
＃子集
gr <-pv.final $ p.value< ; = 0.05
CNb [gr]
＃plot
png（filename =test.png，height = 1000，width = 2000）
print（＃or ggsave（ ）
ggplot（alpha.sub，aes（x = both，y = Value，fill = Group））+ geom_boxplot（）+
#facet_grid（〜Group，scales =free，space = （比较= CNb [gr]，
vjust = 0.7，
annotation = pv.final $ map.signif [gr]，
textsize = 10，
size = 1）
）
dev.off（）
产生：
我如何重现上面的这个图，但是在第一个图的方面？谢谢！
解决方案
我使用了
＃构建图，例如获得数据帧列表
myplot2< - ggplot_build（myplot）
＃在列表2中，您可以访问每个注释
head（myplot2 $ data [[2]]）
x xend y yend注释组PANEL形状颜色文字大小角度hjust vjust alpha系列
1 1 1 3.968526 4.063608 0.11 DATABASE1-DATABASE2-1 1 19黑色3.88 0 0.5 0 NA
2 1 2 4.063608 4.063608 0.11 DATABASE1- DATABASE2-1 1 19黑色3.88 0 0.5 0 NA
3 2 2 4.063608 3.968526 0.11 DATABASE1-DATABASE2-1 1 19黑色3.88 0 0.5 0 NA
4 1 1 3.968526 4.063608 0.035 DATABASE1-DATABASE2-1 2 19黑色3.88 0 0.5 0 NA
5 1 2 4.063608 4.063608 0.035数据库1-数据库2-1 2 19黑色3.88 0 0.5 0 NA
6 2 2 4.063608 3.968526 0.035数据库1-数据库2-1 2 19黑色3.88 0 0.5 0 NA
fontface lineheight线型尺寸
1 1 1.2 1 0.5
2 1 1.2 1 0.5
3 1 1.2 1 0.5
4 1 1.2 1 0.5
5 1 1.2 1 0.5
6 1 1.2 1 0.5
＃现在您可以删除或更新注释
＃获取全部ADJUSTED PVALUES
pv < - tidy（使用（alpha.sub，pairwise.wilcox.test（Value，interaction（Group，DB），p.adjust.method =BH）））
＃使用dplyr创建最终数据集
pv_final< - pv％>％
separate（group1，c（g1_1，g1_2））％>％
separate（group2，c（g2_1，g2_2））％> ;％
filter（g1_1 == g2_1）％>％
mutate（p = ifelse（p.value> 0.05，，ifelse（p.value> 0.01，*，**）））
＃每个pvalue在该数据集中重复三次。
myplot2 $ data [[2]] $ annotation< - rep（pv_final $ p，each = 3）
＃删除非有效值
myplot2 $ data [[2]]< - myplot2 $ data [[2]] [myplot2 $ data [[2]] $ annotation！=，]
＃和最终的绘图
绘图（ggplot_gtable（myplot2））
Following up on this question I posted some days ago, I want to perform something similar.
Given the following MWE:
##############################
##MWE
library(ggplot2)
library(ggsignif)
set.seed(1)
alpha.subA <- data.frame(Sample.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))
alpha.subA$DB <- "DATABASE1"
set.seed(2)
alpha.subB <- data.frame(Sample.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))
alpha.subB$DB <- "DATABASE2"
alpha.sub <- rbind(alpha.subA, alpha.subB)
alpha.sub$DB <- as.factor(alpha.sub$DB)
alpha.sub$both <- factor(paste(alpha.sub$Group, alpha.sub$DB), levels=paste(rep(levels(alpha.sub$Group), each=length(levels(alpha.sub$DB))), rep(levels(alpha.sub$DB), length(levels(alpha.sub$Group)))))
png(filename="test.png", height=1000, width=2000)
print(#or ggsave()
ggplot(alpha.sub, aes(x=both, y=Value, fill=Group)) + geom_boxplot() +
facet_grid(~Group, scales="free", space="free_x") +
stat_summary(fun.y=mean, geom="point", shape=5, size=4)
# + geom_signif() ##HOW TO TEST EACH PAIR IN EACH FACET (DATABASE1 vs DATABASE2 PER GROUP)?
)
dev.off()
##############################
which produces:
I would like to make use of ggsignif's geom_signif() to compare DATABASE1 to DATABASE2 per Group (each pair in a facet) and show the significance level (* for p-value<0.05, ** for p-value<0.01) for each pairwise comparison.
Any help would be greatly appreciated. Thanks!!
EDIT!
I have managed to place the significance levels properly, but the moment I add facets, everything crumbles...
See my new MWE
##MWE
library(ggplot2)
library(ggsignif)
library(tidyverse)
library(broom)
set.seed(1)
alpha.subA <- data.frame(Sample.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, mean=2.3, sd=0.45))
alpha.subA$DB <- "DATABASE1"
alpha.subB <- data.frame(Sample.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, mean=2, sd=0.5))
alpha.subB$DB <- "DATABASE2"
alpha.sub <- rbind(alpha.subA, alpha.subB)
alpha.sub$DB <- as.factor(alpha.sub$DB)
alpha.sub$both <- factor(paste(alpha.sub$Group, alpha.sub$DB), levels=paste(rep(levels(alpha.sub$Group), each=length(levels(alpha.sub$DB))), rep(levels(alpha.sub$DB), length(levels(alpha.sub$Group)))))
v1 <- grep("DATABASE1", levels(alpha.sub$both), val=TRUE)[-9]
v2 <- grep("DATABASE2", levels(alpha.sub$both), val=TRUE)[-9]
CNb <- mapply(c, v1, v2, SIMPLIFY=FALSE)
CNb <- unname(CNb)
pv <- tidy(with(alpha.sub[ alpha.sub$Group != "PGMC4", ], pairwise.wilcox.test(Value, both, p.adjust.method = "BH")))#ADJUSTED PVALUES
# data preparation
CNb2 <- do.call(rbind.data.frame, CNb)
colnames(CNb2) <- colnames(pv)[-3]
# subset the pvalues, by merging the CN list
pv.final <- merge(CNb2, 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,"*", "**"))
# subset
gr <- pv.final$p.value <= 0.05
CNb[gr]
# the plot
png(filename="test.png", height=1000, width=2000)
print(#or ggsave()
ggplot(alpha.sub, aes(x=both, y=Value, fill=Group)) + geom_boxplot() +
#facet_grid(~Group, scales="free", space="free_x") +
stat_summary(fun.y=mean, geom="point", shape=5, size=4) +
geom_signif(comparisons=CNb[gr],
vjust=0.7,
annotation=pv.final$map.signif[gr],
textsize=10,
size=1)
)
dev.off()
which produces:
How can I reproduce this plot above, but with the facets in the first plot? Thanks!
解决方案
I used the idea from here. In brief, plot the plot, get the underlying data, update the annotations and plot the final picture with facets.
# first save the plot in variable
myplot <- ggplot(alpha.sub, aes(x=DB, y=Value)) +
geom_boxplot(aes(fill=Group)) +
facet_grid(~Group) +
geom_signif(test="wilcox.test", comparisons = list(c("DATABASE1", "DATABASE2")), map_signif_level = F)
# a plot with all significance layer per facet group.
myplot
# build the plot, e.g. get a list of data frames
myplot2 <- ggplot_build(myplot)
# in list 2 you have access to each annotation
head(myplot2$data[[2]])
x xend y yend annotation group PANEL shape colour textsize angle hjust vjust alpha family
1 1 1 3.968526 4.063608 0.11 DATABASE1-DATABASE2-1 1 19 black 3.88 0 0.5 0 NA
2 1 2 4.063608 4.063608 0.11 DATABASE1-DATABASE2-1 1 19 black 3.88 0 0.5 0 NA
3 2 2 4.063608 3.968526 0.11 DATABASE1-DATABASE2-1 1 19 black 3.88 0 0.5 0 NA
4 1 1 3.968526 4.063608 0.035 DATABASE1-DATABASE2-1 2 19 black 3.88 0 0.5 0 NA
5 1 2 4.063608 4.063608 0.035 DATABASE1-DATABASE2-1 2 19 black 3.88 0 0.5 0 NA
6 2 2 4.063608 3.968526 0.035 DATABASE1-DATABASE2-1 2 19 black 3.88 0 0.5 0 NA
fontface lineheight linetype size
1 1 1.2 1 0.5
2 1 1.2 1 0.5
3 1 1.2 1 0.5
4 1 1.2 1 0.5
5 1 1.2 1 0.5
6 1 1.2 1 0.5
# now you can remove or update the annotation
# get all ADJUSTED PVALUES
pv <- tidy(with(alpha.sub, pairwise.wilcox.test(Value, interaction(Group,DB), p.adjust.method = "BH")))
# create final dataset using dplyr
pv_final <- pv %>%
separate(group1, c("g1_1", "g1_2")) %>%
separate(group2, c("g2_1", "g2_2")) %>%
filter(g1_1 == g2_1) %>%
mutate(p=ifelse(p.value > 0.05, "", ifelse(p.value > 0.01,"*", "**")))
# each pvalue is repeated three times in this dataset.
myplot2$data[[2]]$annotation <- rep(pv_final$p, each=3)
# remove non significants
myplot2$data[[2]] <- myplot2$data[[2]][myplot2$data[[2]]$annotation != "",]
# and the final plot
plot(ggplot_gtable(myplot2))
这篇关于R ggplot2 - 在一个方面执行每对配对测试，并用ggsignif显示p值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
