置換檢驗概述及其在R中的實現

某些差異分析方法建立在觀測數據服從較好的特定理論分布的基礎上，如常見的參數檢驗T檢驗、方差分析等，假定觀測數據抽樣自正態分布；再如分析差異表達基因時最常用的負二項分布模型，假定整體基因表達具有負二項分布特徵。在這些情況下，根據數據的理論分布，即可比較容易建立假設檢驗和總體參數的置信區間估計。但在許多情況中，統計檢驗並不一定滿足，比如數據抽樣於未知或混合分布、樣本量過小、存在離群點、基於理論分布設計合適的統計檢驗過於複雜且數學上難以處理等情況，這時基於隨機化的統計方法就可以派上用場。置換檢驗（Permutation Test），也稱隨機化檢驗或重隨機化檢驗，提供了穩健的檢驗框架。

為理解置換檢驗的邏輯，可首先考慮如下問題。

在兩種處理條件的試驗中，10個個體被隨機分配到其中一種條件（A或B）中，相應的結果變量也被記錄，如下所示。

如果使用T檢驗分析，可假設數據抽樣自正態分布，然後使用雙尾T檢驗來驗證結果。此時，零假設為A處理的總體均值與B處理的總體均值相等，根據數據計算t統計量，將其與理論分布進行比較，如果觀測到的t統計值落在理論分布值的95%置信區間外，將會拒絕零假設，斷言在0.05顯著性水平下兩組的總體均值不相等。

置換檢驗的思路與之不同，為了檢驗兩種處理方式的差異，其工作方式如下：

（1）首先計算觀測數據的t統計量，稱為t0；

（2）隨機置換數據，即將10個觀測變量取出，再隨機分配5個變量給組A，另外5個給組B，並計算置換後數據的t統計量；如此隨機置換N次，獲得N個置換後數據的t統計量，分別稱為tn（n為第N次置換次數）；

（3）如此可獲得tn的經驗分布，如果t0落在tn經驗分布中間95%區域的外側，或者說t0大於（或小於）95%的tn，則在0.05的顯著性水平下，拒絕兩個處理組的總體均值相等的零假設。

以該圖為例，Observed t-ratio就是t0，橫軸的t-ratio就代表了tn的頻數分布（經驗分布），考慮到t0出現在經驗分布的右側，並由於tn代表一種隨機化的統計量，因此通常認為若t0>95%的tn，則表明t0是顯著的。

並且此時t0<tn的概率就是p值，即

其中nx是置換數據後所獲得統計量的零分布期望值高於觀測值（上圖紅色區域）的置換次數，N是置換總次數。該公式還解釋了為什麼通常在置換檢驗中，使用的置換次數總設置以「9」結尾（例如99、199、499，而不使用100、200、500）：由觀測數據所得的初始統計量也被認為是零分布的一部分，因此「+1」。

 

二者相比，儘管傳統的T檢驗和置換檢驗過程都計算了相同的t統計量，但置換方法並不是將統計量與理論分布進行比較，而是將其與置換觀測數據後獲得的經驗分布進行比較，根據統計量的極端性判斷是否有足夠的理由拒絕零假設。如果數據偏離正態分布，存在離群值，或者感覺對於標準的參數方法來說數據集太小，那麼置換檢驗則是一個不錯的選擇。

由於置換檢驗通過隨機置換數據的方式實現評估觀測數據結構，因此可知每次置換後數據集的變量分布都是不一致的，即每次所得的隨機統計量都有所不同，所以對於同一數據，每次通過置換檢驗所得的p值會略有差異。再次以該圖為例，也就是每次所得的紅色區域範圍會有差別，但是當置換次數足夠大時，p值的差異將很小。

此外，對於執行置換檢驗的軟體，如R，也可以通過設置隨機數種子，實現結果重現的目的。

 

當置換次數達到一定數值時，數據中所有可能的隨機排列情況均已被考慮完全，即置換次數達到飽和。此時經驗分布依據的是數據所有可能的排列組合，這種情況下的置換檢驗也常被稱為「精確檢驗」。

隨著樣本量的增加，獲取所有可能排列的時間開銷會非常大。這種情況下，可以使用蒙特卡洛模擬，從所有可能的排列中進行抽樣，獲得一個近似值。

不難看出，除了上述演示的代替T檢驗，置換檢驗的邏輯還可以延伸至很多統計檢驗或線性模型等方法上來，為它們提供非參數的檢驗框架，用以評估結果是否是合理的。

依靠基礎的抽樣分布理論知識，置換檢驗提供了另外一種強大的可選檢驗思路，由於其不受數據分布特徵的限制，因此它的應用更為靈活。

儘管如此，置換檢驗主要發揮作用的地方仍是處理數據偏倚程度較大（如偏離非正態性）、存在離群點、樣本很小或無法做參數檢驗等情況。並且如果初始樣本對總體情況代表性很差，即使置換檢驗也無法提高推斷效果。此外，置換檢驗主要用於生成檢驗零假設的p值，有助於回答「效應是否存在」，但對於獲取置信區間和估計測量精度則比較困難。

置換檢驗的思想在數十年前就已被提出，但由於置換數據的計人工計算比較費力，限制的它的使用。直到計算機普及後，該方法才得到了真正的應用價值。

在前述已寫過的推文中，很多分析中都借鑑了置換檢驗的思想。

例如在基於距離的差異檢驗方法中，提到的幾種方法均首先計算數據的一種初始統計量，然後再隨機置換數據計算置換後數據的統計量獲得經驗分布，最後比較初始統計量在經驗分布中的區間位置，評估檢驗的顯著性。

再例如，非約束排序中被動添加解釋變量一文中，將外部解釋變量以多元回歸的方式被動擬合到非約束排序軸中，可獲得擬合的R2，其代表了外部解釋變量與非約束軸的關聯程度。之後再通過置換檢驗的方式，獲得隨機置換後數據的擬合R2，比較觀測R2與隨機R2的關係，若觀測R2明顯大於95%以上的隨機R2，則代表這種關聯是可靠的。與此類似，約束排序中確立約束軸的有效性時，也廣泛使用了置換檢驗。

再例如協慣量分析，統計量RV代表了兩個數據矩陣之間的pearson相關性，對於兩矩陣結構相似程度的有效性的評估，置換檢驗就提供了一種便捷的手段。通過比較觀測數據集的RV和置換數據後所得的RV，若觀測RV明顯大於95%以上的隨機RV，則代表兩組數據集的潛在相關是合理的。

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

接下來展示置換檢驗在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.