R统计-单因素ANOVA/Kruskal-Wallis置换检验

1 基本信息

许多时候数据并不能满足许多统计假设，比如数据抽样于未知或混合分布，样本量过小、存在离群点、基于理论分布设计的合适的统计检验过于复杂且数学上难以处理的情况。这时可以使用基于随机化和重抽样的统计方法进行检验。

本文介绍一种应用广泛的依据随机化思想的统计方法：置换检验。置换方法与参数方法都需要计算检验统计量，但是置换方法并不是将统计量与理论分布进行比较，而是将其与置换观测数据后获得的经验分布进行比较，根据统计量值的极端性判断是否有足够的理由拒绝零假设。

置换检验(随机化检验或重随机化检验)的研究思路：

如果处理之间是等价的，那么样本分配的处理标签就应该是任意的，检验处理间是否存在差异的步骤：

1）与参数方法类似，计算观测数据的检验统计量，比如t0;

2) 将所有样本数据归为1组；

3）随机分配与原始各分组样本相同的样本给各个处理；

4）计算并记录新观测数据的检验统计量；

5）对每一种可能的随机分配重复3）-4）步，存在多种可能的分配组合；

6）将每个分配组合的检验统计量按照升序排列，这便是基于样本数据的检验分布；

7）如果实际的检验统计量，如t0,落在经验分布中间95%部分的外面，则在0.05的显著性水平下，拒绝组间总体均值相等的零假设。

2 分析流程

2.1 设置工作路径并调用R包

# 设置工作路径
#knitr::opts_knit$set(root.dir="D:\\EnvStat\\公众号文件\\diff_stat")# 使用Rmarkdown进行程序运行
setwd("D:\\EnvStat\\公众号文件\\diff_stat") # 设置工作目录
Sys.setlocale('LC_ALL','C') # Rmarkdown全局设置
​
# 调用R包
library(coin)

2.2 数据准备

使用虚构数据环境因子数据进行分析。

# 导入数据
data = read.csv("data.csv",check.names = FALSE,header = TRUE,row.names = 1)
data
​
# 将分类变量设置为factor
data$depth = factor(data$depth)
data$T = factor(data$T)
data$N = factor(data$N)

图1|data数据。三个分类变量，depth包含3个水平，T和N包含2个水平。

2.3 单变量两样本和K样本置换检验

下面介绍使用coin包进行两样本和K样本置换检验。

置换检验的经验分布如果依据的是数据所有可能的排列组合，则称为"精确"检验(两样本置换检验才能使用)，如果样本量很大，为减少计算量，可以使用蒙特卡洛模拟，从所有可能的排列中进行抽样，获得一个近似的检验，但其是使用伪随机数从所有排列组合中进行抽样，因此，每次检验的结果都有所不同，大样本近似随机检验，可使用固定随机值，保证结果的可重现性。

2.3.1 单变量两样本和K样本置换检验

这里计算统计量的方法可以选择参数或非参数检验方法，置换检验计算p值。coin包中有三种方法在零假设条件下形成经验分布：1）exact：依据所有可能的排列组合形成经验分布，仅在两样本检验时可用；2）asymptotic：使用渐进分布形成经验分布；3）approxiamate:使用蒙特卡洛重抽样形成近似经验分布。

# 两样本参数置换检验
##Exact Two-Sample Fisher-Pitman Permutation Test
a1 =coin::oneway_test(v1 ~ T,data=data,
                      distribution = exact())
a1
statistic(a1) # 输出统计值
pvalue(a1) # 输出p值及其置信区间

图2|两样本参数置换检验结果。

# 两样本非参数置换检验
##Exact Wilcoxon-Mann-Whitney Permutation Test
a2 =coin::wilcox_test(v1 ~ T,data=data,
                      #distribution = approximate(nresample = 9999), # 蒙特卡洛重抽样
                      distribution =  "exact")
a2
​
statistic(a2) # 输出统计值
pvalue(a2) # 输出p值及其置信区间

图3|两样本Exact Wilcoxon-Mann-Whitney置换检验结果。"Here, the conditional Wilcoxon-Mann-Whitney test was performed via a rank transformation of the response, employing the exact distribution for obtaining the p-value。"

2.3.2 单变量单因素置换检验及多重比较

# 单因素方差置换检验
##Permutation test for One-Way ANOVA
##Monte Carlo resampling using 10000 replicates  
set.seed(12345) # 设计随机重抽样，需要设置随机种子，保证结果可重现。
b1 =coin::oneway_test(v1 ~ depth,data=data,
                      distribution = approximate(nresample = 10000))
b1
b1@statistic@df # 自由度
statistic(b1) # 输出统计值
pvalue(b1) # 输出p值及其置信区间

图4|One-Way ANOVA置换检验结果。

# 单因素非参数置换检验
##Kruskal-Wallis Tests Permutation Test
##Monte Carlo resampling using 10000 replicates  
set.seed(12345)
b2 =coin::kruskal_test(v1 ~ depth,data=data,
                      distribution = approximate(nresample = 10000))
b2
b2@statistic@df # 自由度
statistic(b2) # 输出统计值
pvalue(b2,
         method = "step-down",
         distribution = "marginal",
         type="Bonferroni") # 输出p值及其置信区间

图5|Kruskal-Wallis置换检验结果。"Here, the exact null distribution of the Kruskal-Wallis test is approximated by Monte Carlo resampling using 10000 replicates via"

# 参数单因素置换检验多重比较
set.seed(12345)
it1 = independence_test(v1 ~ depth, data = data,
xtrafo = mcp_trafo(depth = "Tukey"),
distribution = approximate(nresample = 10000))
p.value = pvalue(it1,method = "step-down",
         distribution = "marginal",
         type="Bonferroni") # Bonferroni step-down (Holm) adjust p-values
res1 = data.frame(comparison = rownames(statistic(it1,type="standardized")),
           statistic= statistic(it1,type="standardized"),
           p.value) # 提取统计结果
colnames(res1) <- c("comparison","statistic","p.value")
res1

图6|单因素ANOVA置换检验多重比较结果。

# 非参数单因素置换检验多重比较
set.seed(12345)
it2 = independence_test(v1 ~ depth, data = data,
xtrafo = mcp_trafo(depth = "Tukey"),
ytrafo = function(data)
  trafo(data, numeric_trafo = rank_trafo),
distribution = approximate(nresample = 10000))
​
res2 = data.frame(comparison = rownames(statistic(it2,type="standardized")),
           statistic= statistic(it2,type="standardized"),
           p.value=pvalue(it2,method = "step-down",
         distribution = "marginal",
         type="Bonferroni") 
         ); # 提取统计结果
res2

图7|单因素Kruskal-Wallis置换检验多重比较结果。

图8|使用coin包进行多重比较会报警告。介意的话，可以进行自行对样本进行两两间比较，然后校正p值。

微信公众号后台回复"单因素差异置换检验"或QQ群文件获取数据和代码。

R统计-单因素ANOVA/Kruskal-Wallis置换检验

参考资料

[1] Hothorn T, Hornik K, van de Wiel MA, Zeileis A (2006). "A Lego system for conditional inference." _The American Statistician_,*60*(3), 257-263. doi: 10.1198/000313006X118430 (URL: https://doi.org/10.1198/000313006X118430).

[2] Hothorn T, Hornik K, van de Wiel MA, Zeileis A (2008). "Implementing a class of permutation tests: The coin package." _Journal of Statistical Software_, *28*(8), 1-23. doi: 10.18637/jss.v028.i08 (URL: https://doi.org/10.18637/jss.v028.i08).

[3]《R语言实战》

---

推荐阅读

R绘图-物种、环境因子相关性网络图(简单图、提取子图、修改图布局参数、物种-环境因子分别成环径向网络图)

R统计绘图-分子生态相关性网络分析(拓扑属性计算，ggraph绘图)

R中进行单因素方差分析并绘图

R统计绘图-多变量单因素非参数差异检验及添加显著性标记图

R统计绘图-单因素Kruskal-Wallis检验

R统计绘图-单、双、三因素重复测量方差分析[Translation]

R统计-多变量单因素参数、非参数检验及多重比较

R统计-多变量双因素参数、非参数检验及多重比较

R绘图-KEGG功能注释组间差异分面条形图

R绘图-相关性分析及绘图

R绘图-相关性系数图

R统计绘图-环境因子相关性热图

R绘图-RDA排序分析

R统计-VPA分析(RDA/CCA)

R统计绘图-RDA分析、Mantel检验及绘图

R统计-PCA/PCoA/db-RDA/NMDS/CA/CCA/DCA等排序分析教程

R统计-正态性分布检验[Translation]

R统计-数据正态分布转换[Translation]

R统计-方差齐性检验[Translation]

R统计-Mauchly球形检验[Translation]

R统计绘图-混合方差分析[Translation]

R统计绘图-协方差分析[Translation]

R统计绘图-RDA分析、Mantel检验及绘图

R统计绘图-factoextra包绘制及美化PCA结果图

R统计绘图-环境因子相关性+mantel检验组合图(linkET包介绍1)

R统计绘图-随机森林分类分析及物种丰度差异检验组合图

机器学习-分类随机森林分析(randomForest模型构建、参数调优、特征变量筛选、模型评估和基础理论等)

---