多重检验标字母 ｜ R
2024-05-09 本文已影响0人
kkkkkkang
两两比较显著性，虽然可以用ggpubr（可参考https://www.jianshu.com/p/f53a05da7745），可是分组一多起来就比较麻烦了，看着也很乱。用字母标注比较好。
在多重比较检验时需要根据数据特征选择参数或者非参数检验方法，简单说来就是服从正态分布的可以用参数检验方法如Tukey法或非参数检验方法如Wilcox法，不服从正态分布的则用非参数检验法如Wilcox法
如果比较的两组中有任何一组不服从正态分布，非官方建议：使用非参数检验方法。
首先说如何检验是否正态分布，并对服从正态分布的数据进行检验和画图
library(devtools)
remotes::install_github("taowenmicro/EasyStat")
library(EasyStat)
library(ggplot2)
library(dplyr)
library(ggsci) #配色
NorNorCVTest(data = data_wt, i= 6, method_cv = "leveneTest")[[1]]
# No Name W p.value norm.test
# 1 1 CF 0.8087115 0.07029506 Norm
# 2 2 CK 0.8087115 0.07029506 Norm
# 3 3 Rhi 0.9399859 0.65908550 Norm
# 4 4 WT 0.8087115 0.07029506 Norm
# 5 Test Method: Shapiro-Wilk NA NA <NA>
# 可以看出各组均服从正态分布
# 作者这里检验正态性是调用了shapiro.test()函数，感兴趣的可以查看源码
result <- aovMcomper(data = data_wt, i= 6,method_Mc = "Tukey")
sig <- result[[1]]
sig
# groups group
# CF a CF
# CK a CK
# Rhi b Rhi
# WT a WT
# 加上每个组别的平均值信息，方便后面字母标位置
fun_mean <- data_wt %>% group_by(group) %>% summarise(f_mean = mean(fun)) %>% as.data.frame()
fun_mean
# group f_mean
# 1 CF 46500
# 2 CK 58500
# 3 Rhi 243000
# 4 WT 43000
label_data <- merge(sig, fun_mean)
label_data
# group groups f_mean
# 1 CF a 46500
# 2 CK a 58500
# 3 Rhi b 243000
# 4 WT a 43000
ggplot(data_wt) + geom_boxplot(aes(x = group, y = fun, color = group)) + geom_text(data = label_data, aes(x = group, y = f_mean+ 40000, label = groups)) + scale_color_d3() + theme_classic() + xlab("") + theme(axis.title.y = element_text(face = "bold"), legend.title = element_blank())
image.png
另外再看一下这个用非参数检验结果如何
res <- KwWlx2(data = data_wt, i= 6)[[1]]
res
# groups group
# Rhi a Rhi
# CK b CK
# CF bc CF
# WT c WT
# 这里KwWlx函数有点问题，统计检验和字母标记对不上，看了一下源码应该是ord_sig函数那里有问题，已提交issue。作者修复前不要用。
label2 <- merge(res, fun_mean)
label2
# group groups f_mean
# 1 CF bc 46500
# 2 CK b 58500
# 3 Rhi a 243000
# 4 WT c 43000
ggplot(data_wt) + geom_boxplot(aes(x = group, y = fun, color = group)) + geom_text(data = label2, aes(x = group, y = f_mean+ 40000, label = groups)) + scale_color_d3() + theme_classic() + xlab("") + theme(axis.title.y = element_text(face = "bold"), legend.title = element_blank())
image.png
总体上和参数检验相似
然后我又自己编了一组数据，有不服从正态分布的
data_wt$new <- c(1,1,1,12,13,14,5,6,7,8,4,4,5,9,6,8,10,9,10,20,59,61,62,63)
NorNorCVTest(data = data_wt, i= 13, method_cv = "leveneTest")[[1]]
# No Name W p.value norm.test
# 1 1 CF 0.7396239 0.01573710 Other_situation
# 2 2 CK 0.9199447 0.50496262 Norm
# 3 3 Rhi 0.9123751 0.45220542 Norm
# 4 4 WT 0.7283279 0.01212355 Other_situation
# 5 Test Method: Shapiro-Wilk NA NA <NA>
res2 <- KwWlx2(data = data_wt, i = 13)[[1]]
res2
# groups group
# WT a WT
# Rhi b Rhi
# CK b CK
# CF b CF
new_mean <- data_wt %>% group_by(group) %>% summarise(n_mean = mean(new)) %>% as.data.frame()
new_mean
# group n_mean
# 1 CF 7.000000
# 2 CK 5.666667
# 3 Rhi 7.833333
# 4 WT 45.833333
label3 <- merge(res2, new_mean)
label3
# group groups n_mean
# 1 CF b 7.000000
# 2 CK b 5.666667
# 3 Rhi b 7.833333
# 4 WT a 45.833333
ggplot(data_wt) + geom_boxplot(aes(x = group, y = new, color = group)) + geom_text(data = label3, aes(x = group, y = n_mean+ 20, label = groups)) + scale_color_d3() + theme_classic() + xlab("") + theme(axis.title.y = element_text(face = "bold"), legend.title = element_blank())
image.png
参考：https://www.guhejk.com/wordpress/?p=7051
https://mp.weixin.qq.com/s/Rqlr0WqutKCfdz7Ij2R2Lw
https://mp.weixin.qq.com/s/LV2W5aen9QSRcs4wB7MgDA
https://github.com/taowenmicro/EasyStat/issues/3
