置換檢驗概述及其在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

##使用 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)