R 語言之數據分析「Resampling」

‍‍‍‍‍‍‍‍‍

本節主要總結「數據分析」的「Resampling」重抽樣思想，並通過 R 語言實現。

有一種東西叫作「傳統」，它在很多時候很有用，但會讓你思維固化，在新的環境下讓你出錯。

在總結回歸分析和方差分析的時候 ④Ｒ語言之數據分析「初章」，我總是會在模型的建立之前提到「統計假設」，在模型建立之後進行「假設檢驗」，原因想必大家都能理解，就是因為這些「統計假設」是我們模型建立思想的基礎，是支撐我們模型正確性的「必要條件」。但是，不可否認的是，這些「必要條件」最終會成為我們「數據分析」的局限，讓我們對「不滿足條件的數據集」束手無策。

本節，我們就來解決這個「必要條件」中的其中一條假設，從特例到普遍。

統計假設中有一條，叫做「假定觀測數據抽樣自正態分布或者是是其他性質較好的分布」，那麼當數據集抽樣自「未知分布、混合分布」、樣本容量過小或存在離群點時，傳統的統計方法所得到的模型可能就會不那麼準確，原因之前已經講過，這個時候「Resampling」 的思想就出現了。它拋棄了分布的理論，而是完全基於同一個樣本，在這個樣本中多次重複抽樣，然後將所有抽樣的結果作為總體，將原樣本放到其中去檢驗，判定其顯著性。因為需要多次重複抽樣，所有它被稱為重抽樣「Resampling」。

1. Resampling 的分類

1.1. 置換檢驗（ permutation test ）

這個方法是傳統統計方法的創建者 R. A. Fisher 建立的，但是由於這個方法的計算量過大、且計算機技術也未成熟，他最後放棄了這個方法。但是，數十年後的今天，計算機技術的高速發展，這個方法終於能夠實現並發揮其價值。

為了更清楚的說明置換檢驗的思想，我舉一個「兩總體均值之差」的推斷問題：

現在有兩套學習方案 A 和 B，在 10 個受試者中隨機抽取 5 個按照方案 A 學習，另外 5 個按照方案 B 學習，在學習完畢後對 10 個受試者進行測試，得到分數如下：              方案 A                 方案 B          91                    89          88                    78          76                    93          79                    81          82                    77原假設H_{0} ：兩種方案的總體均值相等 ；備擇假設H_{a}：兩種方案的總體均值不等

這種問題，用參數方法來解決你應該是很熟悉的，現在我來解釋置換檢驗的思想是怎樣的。置換檢驗的思想下，由於我們原假設兩個方案的總體均值是相等的，那麼如果兩種方案真的是等價的話，則 10 個受試者的分數的分配應該是任意的。這麼講你可能還是不太理解，其實在原假設成立的前提下，我們把兩種方案的總體視為等價的，那麼我們同樣也可將這兩個方案的總體視為同一個總體，這樣的話，10 個受試者者的分數在可以在任意位置，而不用在乎這個位置是屬於方案 A 還是方案 B。在這樣的前提下，我們利用重抽樣得到一個「源於已有樣本的抽樣分布」，再利用 t 統計量在這個「源於已有樣本的抽樣分布」上判別初始觀測數據的顯著性，具體步驟如下：

1. 計算初始觀測數據的 t 統計量，記為 t0
2. 將初始觀測數據的 10 個得分，作為一個「重抽樣的總體」，從中抽取 5 個放到 A 組，另五個放到 B 組，並計算此種情況下的 t 統計量
3. 重複步驟 2，直到窮盡所有可能。其實就是一個組合問題，一共有252種可能
4. 將 252 種可能的 t 統計量，按照從小到大排序，得到了「源於已有樣本的抽樣分布」
5. 根據t0在此分布中的位置，和我們確定的顯著性水平，判別原假設我們可否拒絕

1.2. 交叉驗證（ Cross validation ）

在交叉驗證的思想中，一個樣本被隨機的分成兩個或 K 個子樣本，將其中一個作為驗證集，其他 K-1 個子樣本組合成訓練集，在訓練集上獲得回歸方程，而後在驗證集上做出預測，這為一重驗證，而後所有子樣本輪流充當驗證集，如此往復，直至遍歷所有可能。這樣我們就可以得到 K 個「預測誤差」，求其平均值即為所謂的「 CV 值」，所以常說的 CV 值實際上是預測誤差期望Err的一個估計值。

交叉驗證的思想比較容易理解，也是一種十分符合我們直覺的驗證方法，這裡就不舉例說明了。
但值得注意的一點是，K 的取值會影響我們預測的準確性，具體是怎麼影響的涉及到「偏誤」、「方差」以及「計算效率」的權衡，主要是「偏誤」與「計算效率」的權衡，這裡就不深究，感興趣的同學可以看看 「模型選擇的一些基本思想和方法」這篇文章。

不過，我們可以粗略的設置 K 的取值：

大樣本，5 重交叉驗證，對預測誤差的估計已經足夠準確，這裡偏向於增大計算效率

小樣本，10重交叉驗證，為了得到相對準確的估計，這裡考慮損失計算效率

1.3. 刀切法（ Jackknife ）

刀切法的思想，十分適用於分布的分散過寬或者出現異常值的數據，通俗來講就是「既然我們所擁有的數據都只是樣本，都是抽樣得到的，那麼我們在做推斷和估計的時候，扔掉幾個觀測，看看效果如何。」在刀切法中我們反覆的行為就是從樣本中選取一個觀測，把它拋棄，而後利用剩餘的數據做出估計，如此往復，直到每一個觀測都被拋棄過一次。根據每次估計的統計量，取均值，而後將這個均值與原樣本（沒有觀測沒拋棄）上估計到到的統計量想比較，就可以得到「偏差校正刀切估計值」。我的表述可能你能聽懂，但是還理解不深，下面我舉例說明。

因為一個理科生的各科成績之間可能存在某種相關性，現在我們想用用學生的物理成績和化學成績來預測他的數學成績，現在有 50 個觀測：         序號      物理      化學      數學      1        78       87        90      2        88       79        83      3        67       71        78      …        …        …         …      50       89       88        92

1.4. Bootstrap

與刀切法相比，Bootstrap 法的重抽樣思想更加透徹。刀切法的重抽樣次數與觀測數量有關，而 Bootstrap 法的抽樣次數理論上是沒有上限的，將原樣本在計算機性能允許的基礎上儘可能重複抽樣更多次，得到了一個「virtual population 」，從其中我們可以得到統計量的「虛擬潛在分布」，並根據此分布估計置信區間。

舉例說明：

現在有一個容量為 10 的樣本，均值 20，標準差 2，希望根據樣本估計總體均值的置信區間。

按照傳統思想，我們可以假設樣本分布為正態分布，而後根據正態分布計算總體均值的置信區間。但當樣本分布不服從正態分布時，可以利用 Bootstrap 法：
1. 有放回地從樣本中抽取 10 個觀測，並計算其均值（注意由於抽樣是有放回的，重抽樣得到的樣本與原樣本雖然容量相同，但是觀測不一定相同，因為有的觀測可能被多次選擇，有的可能一直都不會被選中）
2. 重複 1，次數為 2000
3. 將 2000 個均值升序排列
4. 確定顯著性水平，計算得到置信區間（例如：要求 95% 的置信區間，找出 2.5% 和 97.5% 的分位點，其中間部分就為 95% 的置信區間）

2. R 語言實現重抽樣

2.1. 置換檢驗

格式：

> function_name( formula , data , distribution = )# formula 是被檢驗變量之間的關係， data 是數據框，distribution 指定重抽樣的方式

這裡的重抽樣方式，除了之前將的按照完全排列組合抽樣之外（ distribution = "exact" ），為了節省計算機資源，可以選擇根據漸進分布抽樣（ distribution = "asymptotic" ），和蒙特卡洛重抽樣（ distribution = "approxiamate(B=#)" ）# 為重複次數，一般要設定隨機種子，以便將來重現。

2.1.1. 兩樣本和 K 樣本的獨立性檢驗（ coin 包 ）

兩樣本獨立性的參數檢驗法（ 置換檢驗 ）

兩樣本獨立性的非參數檢驗法（ 置換檢驗 ）

> wilcox_test(Prob ~ So ,data=UScrime,  distribution="exact")# 這裡表達式 ~ 的兩邊分別為應變量和分組變量# 傳統為 wilcox.test() 函數

K 樣本獨立性的參數檢驗法（ 置換檢驗 ）

> oneway_test(response ~ trt ,data=cholesterol,distribution="approxiamate(B=9999)")# 這裡表達式 ~ 的兩邊分別為應變量和分組變量# 傳統為 單因素方差分析

2.1.2. 列聯表中兩分組變量之間的獨立性檢驗（ coin 包 ）

> Arthritis <- transform(Arthritis,                       Improved=as.factor(as.numeric(Improved)))> set.seed(1234)> chisq_test(Treatment~Improved,data=Arthritis,           distribution=approximate(B=9999))# 這裡表達式 ~ 的兩邊分別為兩個分組變量# 如果使用有序因子，則次函數為線性趨勢檢驗，因此第一行代碼將 Improved 轉換成無序

2.1.3. 數值型變量之間的獨立性檢驗（ coin 包 ）

> states <- as.data.frame(state.x77)> set.seed(1234)> spearman_test(Illiteracy ~ Murder , data=states,              distribution=approximate(B=9999))# 這裡表達式 ~ 的兩邊分別為兩個數值型變量# coin 包中必須為數據框

2.1.4. 兩樣本和 K 樣本的相關性檢驗（ coin包 ）

兩樣本

> wilcoxsign_test(U1~U2,data=UScrime,distribution="exact")

K 樣本

> friedman_test(times ~ methods | block, data = rounding)

2.1.5. 簡單回歸和多項式回歸（ lmPerm 包 ）

簡單線性回歸（置換檢驗）

> set.seed(1234)> fit1 <- lmp(weight~height,data = women, perm="Prob")> summary(fit)# lmp() 函數與 lm() 函數類似，只是多了 perm 參數，其可選值有 Exact、Prob 或SPR# Exact 根據所有可能的排列組合生成精確檢驗。Prob 從所有可能的排列中不斷抽樣，直至估計的標準差在估計的 p 值 0.1 之下，判停準則由可選的 Ca 參數控制。SPR 使用貫序概率比檢驗來判斷何時停止抽樣。#  觀測數大於10，perm="Exact"將 自動默認轉為 perm="Prob"，因為精確檢驗只適用於小樣本問題。

多項式回歸（置換檢驗）

> set.seed(1234)> fit2 <- lmp(weight~height+I(height^2),data=women,perm="Prob")> summary(fit)

2.1.6. 多元回歸（置換檢驗）

> set.seed(1234)> states <- as.data.frame(state.x77)> fit <- lmp(Murder ~ Population + Illiteracy + Income + Frost,           data=states, perm = "Prob")> summary(fit)

2.1.7. 單因素方差分析和協方差分析

單因素方差分析（置換檢驗）

> set.seed(1234)> fit <- aovp(response ~ trt, data = cholesterol, perm="Prob")> summary(fit)# aovp() 函數與 aov() 函數類似，只是多了 perm 參數，其可選值有 Exact、Prob 或 SPR

單因素協方差分析（置換檢驗）

> set.seed(1234)> fit <- aovp(weight ~ gesttime + dose, data=litter, perm="Prob")> summary(fit)

2.1.8. 雙因素方差分析

> set.seed(1234)> fit <- aovp(len~supp*dose, data = ToothGrowth, perm="Prob")> summary(fit)

2.2. 交叉驗證

《 R 語言實戰 》中的自編函數 shrinkage ( ) 可以實現 K 重交叉驗證，默認為 10 重，可以通過參數 k 修改。

> states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])     > fit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data=states)     > shrinkage(fit)

2.3. 刀切法

《 R 語言實踐 》中並未介紹刀切法的實現，可以參考以下代碼

# 抽樣得到樣本> set.seed(1234)> x <- rnorm(100,0,5)> # 運用刀切法重抽樣> jack <- function(x){+   jackknife <- 0+   for(i in 1:length(x)){+     jackknife[i]=length(x)*var(x)-(length(x)-1)/length(x)*sum(var(x[-i]))+   }+   return(jackknife)+ }> # 計算刀切法得到的抽樣方差> mean(jack(x))/length(x)[1] 24.97106> # 計算抽樣方差> var(x)[1] 25.22075

2.4. Bootstrap

格式：

# Bootstrap 重抽樣，並返回統計量> bootobject <- boot(data=, statistic=, R=, ...)# data 為原樣本，也是我們重抽樣的數據來源，可以是向量、矩陣或者數據框# statistic 為生成 k 個統計量的函數，這個函數的參數中必須包括 indices，indices 參數用於 boot() 對原樣本重抽樣# R 為重抽樣的次數# 獲取置信區間> boot.ci(bootobject, conf=, type= )# bootobject 為存儲重抽樣結果的對象# conf 為置信區間的置信度# type 為置信區間的類型

2.4.1. 對單個統計量使用自助法

> #獲取 R 平方值的自編函數> rsq <- function(formula, data, indices){+           d <- data[indices,]+           fit <- lm(formula, data=d)+           return(summary(fit)$r.square)+ }> #自助抽樣> set.seed(1234)> results <- boot(data=mtcars, statistic = rsq,+                 R=1000, formula=mpg~wt+disp)> #輸出 boot 對象> print(results)> plot(results)> #計算置信區間> boot.ci(results, type=c("perc","bca"))# type 參數為 perc 時展示的是樣本均值，為 bca 時將根據偏差對區間做簡單調整

2.4.2. 多個統計量的自助法

> #返回回歸係數向量的自編函數>  bs <- function(formula, data, indices){+    d <- data[indices,]+    fit <- lm(formula,data=d)+    return(coef(fit))+  }> #自助抽樣>  set.seed(1234)>  results <- boot(data=mtcars, statistic = bs,+                  R=1000, formula=mpg~wt+disp)>  #輸出 boot 對象>  print(results)>  plot(results,index=2)>  #計算置信區間>  boot.ci(results,type = "bca", index=2)>  boot.ci(results,type = "bca", index=3)# index 參數用於指明需要分析的重抽樣所得樣本的列

公眾號後臺回復關鍵字即可學習

回復 R                  R語言快速入門及數據挖掘 
回復 Kaggle案例  Kaggle十大案例精講（連載中）
回復 文本挖掘      手把手教你做文本挖掘
回復 可視化          R語言可視化在商務場景中的應用 
回復 大數據         大數據系列免費視頻教程 
回復 量化投資      張丹教你如何用R語言量化投資 
回復 用戶畫像      京東大數據，揭秘用戶畫像
回復 數據挖掘     常用數據挖掘算法原理解釋與應用
回復 機器學習     人工智慧系列之機器學習與實踐
回復 爬蟲            R語言爬蟲實戰案例分享