R:在 boxplot/GGPLOT 上添加 p 值(emmeans 的 p 值与 "emmeans" 函数形成对比)

R: Add p-values on boxplot/GGPLOT (p-values of emmeans contrasts from "emmeans" function)

我拟合了一个线性混合模型(带有 lmer 函数)。 “比例”作为我的响应变量,“阶段”作为固定因素(有 5 个级别),主题作为随机因素。

LMM.fit<-lmer(Proportion~Stage+(1|Subject),data=tab)

       Chisq Df Pr(>Chisq)    
Stage 290.07  4  < 2.2e-16 ***

Post 使用 emmeans 函数进行临时成对比较以确定不同阶段发生差异的位置:

emmeans(LMM.fit, pairwise ~ Stage)$contrasts

 contrast        estimate     SE     df t.ratio p.value
 Stage1 - Stage2    0.178 0.0505 105087   3.528  0.0038
 Stage1 - Stage3    0.601 0.0363 565787  16.529  <.0001
 Stage1 - Stage4    0.438 0.0496 224258   8.833  <.0001
 Stage1 - Stage5    0.272 0.0497 295762   5.472  <.0001
 Stage2 - Stage3    0.422 0.0514 103631   8.216  <.0001
 Stage2 - Stage4    0.260 0.0615  95267   4.224  0.0002
 Stage2 - Stage5    0.094 0.0620 107696   1.517  0.5516
 Stage3 - Stage4   -0.163 0.0496 229001  -3.279  0.0092
 Stage3 - Stage5   -0.328 0.0491 314661  -6.683  <.0001
 Stage4 - Stage5   -0.166 0.0599 190576  -2.769  0.0446

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 5 estimates

我基本上想在上面显示的箱线图上添加 emmeans 结果中显示的 p 值(在同一图中两两分组之间)。我知道有函数 stat_pvalue_manual() 但我很难知道如何将它与 emmeans contrasts output

一起使用

这是我设法做到的:

source("http://news.mrdwab.com/install_github.R")
install_github("mrdwab/overflow-mrdwab") 
install_github("mrdwab/SOfun")
library(SOfun)

LMM.fit<-lmer(Proportion~Stage+(1|Subject),data=tab)

p.val.test<-pwpm(emmeans(LMM.fit,  ~ Stage),means = FALSE, flip = TRUE,reverse = TRUE)
p.val.test<-sub("[<>]", "", p.val.test)
p.matx<-matrix(as.numeric((p.val.test)),nrow = length(p.val.test[,1]),ncol = length(p.val.test[,1])) #if your factor has 5 levels ncol and nrow=5
rownames(p.matx) <- colnames(p.matx) <-colnames(p.val.test)
p.matx[upper.tri(p.matx, diag=FALSE)] <- NA
stat.test<-subset(melt(p.matx),!is.na(value))
names(stat.test)<-c("group1","group2","p.adj")
stat.test[stat.test$p.adj<=0.001,"p.adj.signif"]<-"***"
stat.test[stat.test$p.adj>0.001 & stat.test$p.adj<=0.01,"p.adj.signif"]<-"**"
stat.test[stat.test$p.adj>0.01 & stat.test$p.adj<=0.05,"p.adj.signif"]<-"*"
stat.test[ stat.test$p.adj>0.05,"p.adj.signif"]<-"ns"
stat.test<-mc_tribble(stat.test) # copy & paste the result of this line before the line of ggboxplot!!

stat.test <- tribble(
        ~group1, ~group2, ~p.adj, ~p.adj.signif,
        "Stage2","Stage1",0.0038,"**",
        "Stage3","Stage1",1e-04,"***",
        "Stage4","Stage1",1e-04,"***",
        "Stage5","Stage1",1e-04,"***",
        "Stage3","Stage2",1e-04,"***",
        "Stage4","Stage2",2e-04,"***",
        "Stage5","Stage2",0.5516,"ns",
        "Stage4","Stage3",0.0092,"**",
        "Stage5","Stage3",1e-04,"***",
        "Stage5","Stage4",0.0446,"*")

bxp<-ggboxplot(tab, x = "Stage", y = "Proportion",color = "Stage", palette = "jco",add = "jitter",bxp.errorbar=T) 
bxp + stat_pvalue_manual(stat.test,
                         y.position = max(tab$Proportion)+sd(tab$Proportion),
                         color = "midnightblue")

PS:this website 帮助 (“GGPUBR:如何将在其他地方生成的 P 值添加到 GGPLOT”)