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”)。
我拟合了一个线性混合模型(带有 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”)。