多重校正p值与post-hoc多重比较p值-美图全攻略
连续变量的组间差异比较是科研中经常遇到的统计问题，事后post-hoc 两两比较的校正p值及p值自动绘图展示是门大学问。箱线图、小提琴图、云雨图、蜂群图等，都是此类分析常见的可视化形式。小编之前在p值FDR多重校正与可视化,实战医学统计—R语言数据可视化之美读书笔记中简要介绍过相应的绘图过程，但是并未全面介绍多重比较实操，涉及多组比较的图形回溯如下。
一、*与abc标记法
对于多组的组间差异分析，小提琴图不仅清楚地展示了各组的数据分布，还可展示组间两两事后多重比较的显著性，然而在采用的于组间比较的“ *,**,***”显著性标记方式，在分组过多时会就会显得过于紧凑、冗余，例如下图。
abc字母标记法是另一种常见的显著性标记方法。但这种方法也存在一定的问题，因为字母无法对统计差异的强弱进行半定量的表现。标记字母法标记法则：凡是存在相同字母的组别之间都不存在显著差异，字母完全不同的才存在显著差异。
小狐狸文献中包括了上面两种，感兴趣可以围观。
上述两种经典标记方法，都比较适合那么，有没有其他更好的显著性标记方法，能够明确标记事后n组多重比较的显著性，但又让画面看上去更加干净清爽呢？
二、p与adj p
R语言可以进行多种多重比较方法，两两比较的P分散到各包。R内置了一些方法来调整一系列p值，以控制多重比较谬误FWE或控制错误发现率FDR。Holm、Hochberg、Hommel和Bonferroni方法控制了多重比较谬误。这些方法试图限制错误发现的概率（I型错误，在没有实际效果时错误地拒绝无效假设），因此都是相对较保守的。
SAS中提供了方差分析系统的两两比较说明，SAS推荐如果你对几个单独的比较感兴趣，并且不关心多重推论的影响，你可以使用重复t检验或Fisher无保护LSD；如果您对所有成对比较或与对照组的所有比较感兴趣，则应分别使用Tukey或Dunnett检验，以做出最强的可能推断；如果您有较弱的推理要求，特别是如果您不需要均值差的置信区间，则应使用REGWQ方法。最后，如果您同意贝叶斯方法以及Waller和邓肯的假设，则应使用Waller-邓肯检验。如果比较次数相对于平均数较大，Scheffé方法可能比Bonferroni或Sidak方法更有效。对于成对比较，Sidak t检验通常更强大。小编汇总方法学后，总结如下。
三、R实战
从案例中可以看到字母标记法1很好的规避了5组内多重比较太冗余的问题，但是缺点是看不出来究竟是P<0.05 还是<0.01还是<0.001, 而星星**标记法虽然能很好的展示两两P值，但是一个图展示了10个p值，过于累赘。
使用for循环算法，完成显著性符号的添加。通过显著性符号的颜色，即可判断出是与哪一组数据进行比较。例如，TIME_IA列上方的3个显著性符号分别是粉红色、群青色、绿色，说明经校正假阳性的Tukey检验反映出，TIME_IA组与TIME_IR，TIME_IE，TIME_ISM这些组存在显著差异。至于是显著更高还是显著更低，可观察箱线图中间横线所代表的中位数。最后，我们将Dunn事后多重比较的结果保存为Excel文件。文件种包括差异分析的组，BH法矫正的p-value（p.adj），及显著性符号（p.adj.signif）。
R实战代码
#save(data,file="data.rdata")
load("data.rdata")
pacman::p_load(tidyverse,ggpubr,ggsignif,rstatix,car)
#pivot_longer函数将宽格式的数据转换为长格式。长格式更方便进行分组分析。
data_long <- pivot_longer(data, cols = everything(),
names_to = "group", values_to = "Score") %>% as.data.frame()
#保存数据
#write.csv(data_long,"data_long.csv")
#读取数据
#data_long=read.csv("data_long.csv")
#数据进行因子化，赋值为原name
data_long$group <- factor(data_long$group, levels = colnames(data))
#----1 SW QQ图正态----
#正态分布且方差齐-方差分析;不满足可对数据进行预处理或使用非参数检验
pacman::p_load(broom,dplyr)
#各组正态SW检验与p>0.05正态
data_long %>% group_by(group) %>% do(tidy(shapiro.test(.$Score))) %>% select(group, statistic, p.value)
## qq图 点沿着线分布 则满足正态分布
qqPlot <- data_long %>% car::qqPlot(lm(Score ~ group, data=.), simulate=TRUE, main="Q-Q Plot", lables=FALSE)
# 拟合线性模型
qqPlot <- data_long %>%
lm(Score ~ group, data = .) %>%
residuals() %>%
qqPlot(., main = "Q-Q Plot")
#----2 levene_Test 方差齐----
## Bartlett检验，p > 0.05 则方差齐
bartlett_test <- data_long %>% stats::bartlett.test(Score ~ group, data = .)
bartlett_test
# Levene检验,对原始数据的正态性不敏感
levene_Test <- data_long %>% leveneTest(Score ~ group, data = .)
levene_Test
#总体F与p
anova_test <- data_long %>% anova_test(Score ~ group)
anova_test
# 推荐方差分析方法2
##----2.1 ??F齐 F检验----
variance <- data_long %>% aov(Score ~ group, data=.)
summary(variance)
fit <- aov(Score~group, data = data_long)
summary(fit)
#对于 ANOVA 外的其他常见的参数检验方法（这些方法均可继承后续的多重比较）
#fit <- lm(Score ~ group, data = data_long) #线性模型（结果同 ANOVA）
#fit <- glm(Score ~ group, data = data_long, family = 'gaussian')
#广义线性模型（除正态数据外，如果是其他分布的情形，在 family 中更改）
#fit <- lme(Score ~ group, random = ~1|random, data = data_long)
#混合线性模型（如想纳入随机效应的话，在 random 中指定随机效应，nlme 包）
#校正p方法有3套，c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr", "none"
stat.test <- data_long %>%
rstatix::tukey_hsd(
Score ~ group,
p.adjust.method = "BH") %>%
rstatix::add_significance("p.adj", #根据校正后的p value换算显著性符号
#对于显著性符号，小编一般只用到以下几个层级（做多3个星号，默认最多是4个星号）
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))
stat.test
#保存Dunn test结果
write.csv(stat.test,"dunn.csv")
##----2.2 F不齐 非参数检验----
#Kruskal-Wallis 检验，整体差异
fit <- kruskal.test(Score~group, data_long)
fit
###----2.21 F不齐 wilcox.test 两两连线----
library(ggpubr)
# data_long$group 分5组，TIME_IA TIME_ISM TIME_ISS TIME_IE TIME_IR
ggboxplot(data = data_long, x = 'group', y = 'Score') +
stat_compare_means(method = 'wilcox.test',
comparisons = list(
c('TIME_IA', 'TIME_ISM'), c('TIME_IA', 'TIME_ISS'), c('TIME_IA', 'TIME_IE'))
)
#comparisons 中可指定任意组，需要写出分组两两配对
###----2.22 ??F不齐 BH分色----
library(agricolae)
kw <- agricolae::kruskal(data_long$Score,data_long$group, p.adj = 'BH',group=TRUE)
kw
plot(kw)
#----3 标记字母法----
##----3.1 LSD多重比较----
pacman::p_load(agricolae,ggplot2)
#F分析
variance <- data_long %>% aov(Score ~ group, data=.)
#summary(variance)
MC <- agricolae::LSD.test(variance,"group", p.adj="bonferroni")
print(MC$group)
mark <- data.frame(MC$group)
mark$group = row.names(mark)
head(mark)
# 绘制
p.zimu <- ggplot(data_long,aes(group, Score, color = group))+
geom_boxplot()+
geom_jitter()+
scale_color_manual(values = c('#eb4b3a', "#48bad0", "#1a9781", "#355783", "#d4a017"))+
# 添加了第五种颜色
geom_text(data=mark,
aes(x=group,y=Score+2,label=groups),
color="black",
size = 5,
fontface="bold")+#添加字母标记
xlab("")+
theme_classic()+
theme(axis.text.y = element_text(size = 12),
axis.text.x = element_text(size = 12),
axis.title.y = element_text(size = 14,face = 'bold'))
p.zimu
ggsave('zimu.pdf', plot = p.zimu, width = 6, height = 6, family="Times")
##----3.2 ??TukeyHSD多重比较----
#总体fit方差
fit <- aov(Score~group, data = data_long)
res <- TukeyHSD(variance)
plot(res)
#tukey方法2
library(multcomp)
tuk <- cld(glht(fit, alternative = 'two.sided', linfct = mcp(group = 'Tukey')), decreasing = TRUE)
tuk
plot(tuk)
##----3.3 ??Dunnett-t多重比较----
library(multcomp)
#group分组
res <- multcomp::glht(variance, linfct = mcp(group = "Dunnett"))
summary(res)
plot(res)
##----3.4 SNK-q多重比较----
#Student-Newman-Keuls test
res <- agricolae::SNK.test(variance, "group")
res
plot(res)
#snk方法2
library(PMCMRplus)
snk <- snkTest(fit)
summaryGroup(snk)
plot(snk)
##----3.5 Scheffe多重比较----
library(DescTools)
sf <- ScheffeTest(variance)
sf
par(las = 1)
plot(sf)小狐狸款代码请见文末
#----4 **标记法----
##----4.1 bonferroni多重比较----
#两两比较tukey_hsd t_test等等
#校正p method中c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr", "none"
stat.test <- data_long %>%
rstatix::tukey_hsd(
Score ~ group,
p.adjust.method = "BH") %>%
rstatix::add_significance("p.adj", #根据校正后的p value换算显著性符号
#对于显著性符号，小编一般只用到以下几个层级（做多3个星号，默认最多是4个星号）
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))
stat.test
stat.test <- stat.test %>%
rstatix::add_xy_position(x='group',dodge = 1)
# 绘制**标记
p.xx <- ggplot(data_long,aes(group, Score, color = group))+
geom_boxplot()+
geom_jitter()+
scale_color_manual(values = c('#eb4b3a', "#48bad0", "#1a9781","#355783", "#d4a017"))+
stat_pvalue_manual(
stat.test, label = "p = {p.adj}", tip.length = 0
) +
xlab("")+
theme_classic()+
theme(axis.text.y = element_text(size = 12),
axis.text.x = element_text(size = 12),
axis.title.y = element_text(size = 14,face = 'bold'))
p.xx
ggsave('xingxing.pdf', plot = p.xx, width = 6, height = 6, family="Times")
小狐狸款代码请见文末
四、发散思维
FWER控制方法（如Bonferroni、Holm）适用于多种统计检验场景，非常通用但较为保守。FDR控制方法（如Benjamini-Hochberg）适用于处理大量假设检验，如生物信息学中的高通量数据分析，平衡发现真实发现与控制假阳性之间的关系。Tukey的HSD方法及其他事后分析方法专门用于ANOVA之后的两两比较，它们在特定条件下提供了更高的统计功效和适当的错误率控制。
小结
“1.小样本成对比较或与对照组推荐Tukey或Dunnett检验，大样本多重比较推荐BH多重校正，容易出阳性adj P值。2.多组组间差异p值绘图，分组≤3组推荐用**标记款，≥4组可用小狐狸款。
本公众号建立了学术交流群(群4)，仅供SCI学术交流，人数有限需要实名制。二维码7天有效期，过期请加小编微信 popnie，请备注说明：姓名-学校(单位)-专业。格式不对不会通过好友哦！
PS1：需要小狐狸款代码请公众号后台回复 20240219
PS2：小狐狸同学是数据分析师，全职代做图表，需要的可以微信@ jijitoutou