R ggplot2 - ファセットのペアごとにペアワイズテストを実行し、ggsignifでp値を表示します。

Question

フォローアップthis question私は数日前に投稿しました。

############################## 
##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() 
############################## 

生成します：

が

私はggsignifのgeom_signifの使用は（）で（グループあたりのデータベース2に各ペアをDATABASE1を比較するためにしたいと思い、次のMWE考える

各ペアワイズ比較のための有意水準（p値については< 0.05、** p値については< 0.01）を示す。

ご協力いただければ幸いです。ありがとう!!

EDIT！

が生成

##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() 

私の新しいMWEを参照してください...すべてが崩れ、私は適切に有意水準を配置するために管理しているが、今は私がファセットを追加します。

を

上記のプロットはどのように再現できますか？最初のプロットではファセットを使用していますか？ありがとう！

Answer 1

私はhereからアイディアを使用しました。要約すると、プロットをプロットし、基礎となるデータを取得し、注釈を更新し、ファセットで最終画像をプロットします。

# 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))