单因素方差分析的非参数方法在R中实现
对于多组数据间比较差异，我们常会想到使用方差分析（ANOVA）来实现。不过由于ANOVA的前提假设条件比较严格，要求数据必须满足正态性、方差齐性等，而很多情况下我们的数据并不符合方差分析的条件。通常情况下，我们可以考虑转换数据，例如，使用log转换等或许可以使非正态分布的原始数据转变为正态分布类型（当然，我们需要确保转换后的数据能够被合理解释，否则将无意义）；另一种合适的选择是使用非参数的检验方法，替代ANOVA。
由于非参数的检验方法普遍保守，所以当数据满足ANOVA的条件时，尽可能使用ANOVA分析，能够有效地鉴别出非参数检验鉴别不到的差异；当无法适用ANOVA时，再考虑非参数的方法。
在R语言中，常见ANOVA的非参数替代方法例如：
单因素方差分析（单因素ANOVA） => Kruskal-Wallis检验，kruskal.test()；Friedman检验，friedman.test()；
单因素协方差分析（单因素ANCOVA） => sm包sm.ancova()；
双因素方差分析（双因素ANOVA） => Scheirer-Ray-Hare检验，rcompanion包scheirerRayHare()；
多元方差分析（MANOVA） => 置换多元方差分析（PERMANOVA），vegan包adonis()。
本文首先以单因素ANOVA的非参数检验替代方法Kruskal-Wallis检验、Friedman检验为例，简述其在R中的执行过程。其它几种类型，会在后续的几篇文章中再一一阐述。
本文使用的作图数据的网盘链接（提取码fjuq）：
https://pan.baidu.com/s/1I2dX58-q5XJgexxfYTvpmg
Kruskal-Wallis检验
示例数据说明
我们首先将示例数据读到R中，并从中挑选部分数据作为演示（测试数据同前文中的单因素ANOVA，以方便二者比较）。
soil <- read.table('soil.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
soil <- merge(soil, group, by = 'sample')
#以 chao1 指数为例，同时将分组列转换为因子变量
chao1 <- soil[ ,c('sample', 'site', 'chao1')]
chao1$site <- factor(chao1$site)
str(chao1)
head(chao1)
假设存在这么一个关注的问题：
我们在3个地域（A、B、C）分别采集了土壤样本，即获得了3种类型的土壤，并通过16S测序，获得了每种类型土壤中细菌群落的Alpha多样性指数。我们想要得知，三种土壤环境下的细菌群落的Alpha多样性指数是否存在显著不同。
对应于上述挑选出的测试数据“chao1”：sample，采集的土壤样本名称；site，土壤样本来源的环境（A、B、C），这列作为分组列，需要转换为因子变量类型；chao1，Alpha多样性指数中的Chao1指数，数值变量。
在前文方差分析中，我们已经使用了单因素ANOVA对3种土壤环境下的细菌群落的Chao1指数进行比较。在这里，我们更换为Kruskal-Wallis非参数检验比较差异。
Kruskal-Wallis检验
单因素ANOVA要求数据满足正态性和方差齐性，如前文单因素ANOVA中的方法所述。若二者之一不符合，如果各组独立，则Kruskal-Wallis检验将会是一种实用的方法。（尽管本示例数据符合方差分析的前提假设，但这里作为演示Kruskal-Wallis检验，故请忽略这点）
Kruskal-Wallis检验的调用格式为：kruskal.test(y~A, data)。其中y是数值型结果变量（因变量），A是因子类型的分组变量（因子变量）。
#Kruskal-Wallis Test，可使用 ?kruskal.test 查看帮助kruskal.test(chao1~site, data = chao1)
Kruskal-Wallis检验结果表明，3种土壤环境下的细菌群落的Chao1指数具有显著差异，p值远低于0.05水平。
非正态性的数据，不太适合以均值及标准差的形式作为展示，可以使用中位数等代替。
#若想查看各组中位数、极大（小）值、（四）分位数，可使用 aggregate()
aggregate(chao1$chao1, by = list(chao1$site), FUN = max)
aggregate(chao1$chao1, by = list(chao1$site), FUN = function(x) quantile(x, 0.75))
aggregate(chao1$chao1, by = list(chao1$site), FUN = median)
aggregate(chao1$chao1, by = list(chao1$site), FUN = function(x) quantile(x, 0.25))
aggregate(chao1$chao1, by = list(chao1$site), FUN = min)
事后两两比较
上述Kruskal-Wallis检验告诉我们3种土壤环境下的细菌群落的Chao1指数具有显著差异，这种差异是在整体水平而言的，并没有告诉我们究竟谁和谁存在差异。在ANOVA中，我们常使用Tukey HSD检验基于ANOVA的结果继续执行组间两两比较；而在非参数检验中，如果各组独立，则可以考虑使用Wilcoxon秩和检验。
以下使用Wilcoxon秩和检验继续探究两两差异。当分组水平较多时，多次的Wilcoxon秩和检验容易提升I类错误的风险，因此推荐引入p值校正过程。
##此时若想继续探寻两两分组间的差异，可使用 Wilcoxon 秩和检验，如果分组较多，可以使用循环来完成#继续在上述 Kruskal-Wallis 检验的基础上，探究两两分组间差异，一个示例如下
group <- levels(chao1$site)
group1 <- NULL
group2 <- NULL
median1 <- NULL
median2 <- NULL
p <- NULL
for (i in 1:(length(group) - 1)) {
for (j in (i + 1):length(group)) {
group1 <- c(group1, group[i])
group2 <- c(group2, group[j])
group_ij <- subset(chao1, site %in% c(group[i], group[j]))
group_ij$site <- factor(group_ij$site, levels = c(group[i], group[j]))
wilcox_test <- wilcox.test(chao1~site, data = group_ij, alternative = 'two.sided', conf.level = 0.95)
p <- c(p, wilcox_test$p.value)
median1 <- c(median1, median(subset(group_ij, site == group[i])$chao1))
median2 <- c(median2, median(subset(group_ij, site == group[j])$chao1))
}
}
result <- data.frame(group1, group2, median1, median2, p)
result$padj <- p.adjust(result$p, method = 'BH') #推荐加上 p 值校正，这里使用 Benjamini 方法校正 p 值
result
#write.table(result, 'Wilcoxon.txt', sep = '\t', row.names = FALSE, quote = FALSE)
示例结果如下，展示两组比较的组名、中位数、Wilcoxon检验p值以及Benjamini校正后的p值。这样，我们就获得了两组间的差异了。你可以在后续根据p值（或校正后p值）以及中位数（非参数方法比较的普遍是中位数的差异，而非均值差异）数值，手动标注显著性“abc”或者“*”等，不再多说。
可视化示例
stat_compare_means()支持T检验（t.test）、ANOVA（anova）、Wilcoxon秩和检验（wilcox.test）以及Kruskal-Wallis检验（kruskal.test）。若使用T检验或ANOVA，一定要在数据满足前提假设时使用。#使用 ggpubr（ggplot2 的扩展包），作图展示 wilcox 比较的结果
library(ggpubr)
ggboxplot(data = chao1, x = 'site', y = 'chao1', color = 'site') +
stat_compare_means(method = 'wilcox.test', comparisons = list(c('A', 'B'), c('A', 'C'), c('B', 'C')))
Friedman检验
与Kruskal-Wallis检验相比，如果各组不独立（如重复测量设计或随机区组设计），那么Friedman检验可能会更合适。尽管此时你仍然使用Kruskal-Wallis检验也可以。
Friedman检验的调用格式为：friedman.test(y~A|B, data)。其中y是数值型结果变量（因变量），A是因子类型的分组变量（因子变量），而B是一个用以认定匹配观测的区组变量（blocking variable）。
#使用 ?friedman.test 查看帮助文档，这里以帮助文档中最下方提供的示例为例wb <- aggregate(warpbreaks$breaks, by = list(w = warpbreaks$wool, t = warpbreaks$tension), FUN = mean)
friedman.test(x ~ w | t, data = wb)
这个示例的p值不显著。
同样地，关于事后两两比较，我们也可以通过Wilcoxon秩和检验来完成。这里考虑到各组是不独立的，可尝试Wilcoxon符号秩和检验的方法。方法参考上述循环。
对于可视化，我们可以首先计算均值、中位数、显著性p值等，然后通过这些数据再单独使用作图包（如ggplot2）作图展示等。不再多说。