置换检验概述及其在R中的实现

置换检验及其在R中的实现

某些差异分析方法建立在观测数据服从较好的特定理论分布的基础上，如常见的参数检验T检验、方差分析等，假定观测数据抽样自正态分布；再如分析差异表达基因时最常用的负二项分布模型，假定整体基因表达具有负二项分布特征。在这些情况下，根据数据的理论分布，即可比较容易建立假设检验和总体参数的置信区间估计。

但在许多情况中，统计检验并不一定满足，比如数据抽样于未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况，这时基于随机化的统计方法就可以派上用场。

置换检验概述

置换检验（Permutation Test），也称随机化检验或重随机化检验，提供了稳健的检验框架。

置换检验的基本思想

为理解置换检验的逻辑，可首先考虑如下问题。

在两种处理条件的试验中，10个个体被随机分配到其中一种条件（A或B）中，相应的结果变量也被记录，如下所示。

如果使用T检验分析，可假设数据抽样自正态分布，然后使用双尾T检验来验证结果。此时，零假设为A处理的总体均值与B处理的总体均值相等，根据数据计算t统计量，将其与理论分布进行比较，如果观测到的t统计值落在理论分布值的95%置信区间外，将会拒绝零假设，断言在0.05显著性水平下两组的总体均值不相等。

置换检验的思路与之不同，为了检验两种处理方式的差异，其工作方式如下：

（1）首先计算观测数据的t统计量，称为t₀；

（2）随机置换数据，即将10个观测变量取出，再随机分配5个变量给组A，另外5个给组B，并计算置换后数据的t统计量；如此随机置换N次，获得N个置换后数据的t统计量，分别称为t_n（n为第N次置换次数）；

（3）如此可获得t_n的经验分布，如果t₀落在t_n经验分布中间95%区域的外侧，或者说t₀大于（或小于）95%的t_n，则在0.05的显著性水平下，拒绝两个处理组的总体均值相等的零假设。

以该图为例，Observed t-ratio就是t₀，横轴的t-ratio就代表了t_n的频数分布（经验分布），考虑到t₀出现在经验分布的右侧，并由于t_n代表一种随机化的统计量，因此通常认为若t₀>95%的t_n，则表明t₀是显著的。

并且此时t₀<t_n的概率就是p值，即

其中n_x是置换数据后所获得统计量的零分布期望值高于观测值（上图红色区域）的置换次数，N是置换总次数。该公式还解释了为什么通常在置换检验中，使用的置换次数总设置以“9”结尾（例如99、199、499，而不使用100、200、500）：由观测数据所得的初始统计量也被认为是零分布的一部分，因此“+1”。

二者相比，尽管传统的T检验和置换检验过程都计算了相同的t统计量，但置换方法并不是将统计量与理论分布进行比较，而是将其与置换观测数据后获得的经验分布进行比较，根据统计量的极端性判断是否有足够的理由拒绝零假设。如果数据偏离正态分布，存在离群值，或者感觉对于标准的参数方法来说数据集太小，那么置换检验则是一个不错的选择。

关于随机化排列

由于置换检验通过随机置换数据的方式实现评估观测数据结构，因此可知每次置换后数据集的变量分布都是不一致的，即每次所得的随机统计量都有所不同，所以对于同一数据，每次通过置换检验所得的p值会略有差异。再次以该图为例，也就是每次所得的红色区域范围会有差别，但是当置换次数足够大时，p值的差异将很小。

此外，对于执行置换检验的软件，如R，也可以通过设置随机数种子，实现结果重现的目的。

当置换次数达到一定数值时，数据中所有可能的随机排列情况均已被考虑完全，即置换次数达到饱和。此时经验分布依据的是数据所有可能的排列组合，这种情况下的置换检验也常被称为“精确检验”。

随着样本量的增加，获取所有可能排列的时间开销会非常大。这种情况下，可以使用蒙特卡洛模拟，从所有可能的排列中进行抽样，获得一个近似值。

置换检验的广泛应用

不难看出，除了上述演示的代替T检验，置换检验的逻辑还可以延伸至很多统计检验或线性模型等方法上来，为它们提供非参数的检验框架，用以评估结果是否是合理的。

依靠基础的抽样分布理论知识，置换检验提供了另外一种强大的可选检验思路，由于其不受数据分布特征的限制，因此它的应用更为灵活。

尽管如此，置换检验主要发挥作用的地方仍是处理数据偏倚程度较大（如偏离非正态性）、存在离群点、样本很小或无法做参数检验等情况。并且如果初始样本对总体情况代表性很差，即使置换检验也无法提高推断效果。此外，置换检验主要用于生成检验零假设的p值，有助于回答“效应是否存在”，但对于获取置信区间和估计测量精度则比较困难。

置换检验的思想在数十年前就已被提出，但由于置换数据的计人工计算比较费力，限制的它的使用。直到计算机普及后，该方法才得到了真正的应用价值。

在前述已写过的推文中，很多分析中都借鉴了置换检验的思想。

例如在基于距离的差异检验方法中，提到的几种方法均首先计算数据的一种初始统计量，然后再随机置换数据计算置换后数据的统计量获得经验分布，最后比较初始统计量在经验分布中的区间位置，评估检验的显著性。

再例如，非约束排序中被动添加解释变量一文中，将外部解释变量以多元回归的方式被动拟合到非约束排序轴中，可获得拟合的R²，其代表了外部解释变量与非约束轴的关联程度。之后再通过置换检验的方式，获得随机置换后数据的拟合R²，比较观测R²与随机R²的关系，若观测R²明显大于95%以上的随机R²，则代表这种关联是可靠的。与此类似，约束排序中确立约束轴的有效性时，也广泛使用了置换检验。

再例如协惯量分析，统计量RV代表了两个数据矩阵之间的pearson相关性，对于两矩阵结构相似程度的有效性的评估，置换检验就提供了一种便捷的手段。通过比较观测数据集的RV和置换数据后所得的RV，若观测RV明显大于95%以上的随机RV，则代表两组数据集的潜在相关是合理的。

以及更多的方法，不再多说。

R语言中的置换检验

接下来展示置换检验在R语言中的实现方法示例。

使用置换检验替代传统分析方法的R包

R目前有一些非常全面的包可用来做置换检验。

coin包

例如，coin包对于独立性问题提供了一个非常全面的置换检验框架，相对于多种传统的统计检验方法，提供的可选置换检验替代函数。

lmPerm包

lmPerm包则专门用来做方差分析和回归分析的置换检验。它提供了线性模型的置换方法，包括回归（lmp()函数，使用方法类似线性回归函数lm()）和方差分析（aovp()函数，使用方法类似方差分析函数aov()），而非正态理论的检验。

其它常见R包

例如perm包，与coin包中的部分功能类似；corrperm包，提供了有重复测量的相关性的置换检验；logregperm包提供了Logistic回归的置换检验；glmperm，涵盖了广义线性模型的置换检验；以及其它一系列的更多R包等。

自写置换检验过程实现的一个示例

其实在理解了置换检验的原理后（不难理解，对不对），也可以根据这种思想，自写代码去解决一些问题，对它进行广泛的应用。

例如这里以计算某数据中，变量间的Pearson相关系数为例。R函数cor()只能计算相关系数，但无法给出相关系数是否显著，我们借助置换检验的思想，首先计算观测数据的Pearson相关系数（cor0），然后对数据集随机置换并计算置换后数据的Pearson相关系数(corN)，最后根据|corN|>|cor0|的概率获得p值。

同时也通过psych包中带显著性检验的相关系数计算方法，比较手写置换检验结果和psych包中函数的计算结果是否一致。

#数据获取自 http://blog.sciencenet.cn/blog-3406804-1173120.html
dat <- read.table('data.txt', sep = '\t', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)

#只包含细菌类群丰度数据，计算细菌类群之间的相关性
dat_phylum <- dat[8:17]
var_name <- names(dat_phylum)

##自写过程实现
#计算观测数据的相关系数
dat_phylum <- scale(dat_phylum)
cor0 <- cor(dat_phylum, method = 'pearson')

p_num <- cor0
p_num[abs(p_num)>0] <- 1

#随机置换数据 999 次，计算每次获得的相关系数，并统计 |corN|>|cor0| 的频数
set.seed(123)

for (i in 1:999) {
random <- matrix(sample(dat_phylum), ncol = 10)
colnames(random) <- var_name
corN <- cor(random, method = 'pearson')

corN[abs(corN) >= abs(cor0)] <- 1
corN[abs(corN) < abs(cor0)] <- 0
p_num <- p_num + corN
}

#p 值矩阵，即 |corN|>|cor0| 的概率
p <- p_num/1000
p

##使用 R 内置函数的相关性系数计算并获得检验
library(psych)

cor_test <- corr.test(dat_phylum, method = 'pearson')

cor_test$r
cor_test$p

##相关图比较
#左图为手写的置换检验结果，右图为 psych 包获得的结果
#仅显著的相关系数标以背景色
library(corrplot)

layout(matrix(c(1,2), 1, 2, byrow = TRUE))
corrplot(cor0, method = 'square', type = 'lower', p.mat = p, sig.level = 0.05, insig = 'blank', addCoef.col = 'black', diag = FALSE, number.cex = 0.8, tl.cex = 0.8)
corrplot(cor_test$r, method = 'square', type = 'lower', p.mat = cor_test$p, sig.level = 0.05, insig = 'blank', addCoef.col = 'black', diag = FALSE, number.cex = 0.8, tl.cex = 0.8)

手写置换检验结果（左图）和psych包中函数的计算结果（右图）是一致的。

参考资料

Robert I. Kabacoff. R语言实战（第二版）（王小宁 刘撷芯 黄俊文 等 译）. 人民邮电出版社, 2016.

友情链接