《R語言實戰》第三部分第十二章-重抽樣與自助法學習筆記
以下我們將重溫一下用傳統方法分析過的問題，然後對比看看本章介紹的新方法如何解決它們。
1 置換檢驗
置換檢驗（permutation test），也稱隨機化檢驗或重隨機化檢驗。書本上舉了一個例子，如果看得不是特別明白可以再看看以下這個例子
Permutation Test 置換檢驗
首先，它也與之前的參數方法類似，構成並計算一個觀測數據的t統計量，並計算出t0，後面就有所不同，也是它的本質所在。它將所有的樣本放在一起，重新隨機分配的兩個不同的組，計算每次隨機到的新觀測的統計量，如果有N種可能的隨機分配組織，則有N個新觀測的統計量，我們把每次計算得到的新觀測的統計量的分布記為經驗分布。最後根據t0落在經驗分布中的位置來確定是否拒絕之前的零假設。
需要注意的是，雖然置換方法和參數方法都計算了相同的t統計量，但是參數方法是將統計量與理論分布進行比較，置換方法是將統計量與置換觀測數據後獲得的經驗分布進行比較，再根據統計量值得極端性來判斷是否有足夠的理由拒絕零假設。
根據經驗分布所依據的數據是否是所有可能的排列組合，置換檢驗可以分為「精確」檢驗和非「精確」檢驗，對於後者可以使用蒙特卡羅模擬（蒙特卡羅方法入門 - 阮一峰的網路日誌），從所有可能的排列當中進行抽樣，獲得一個近似的檢驗，特別是在樣本量很大，減少可能的計算開銷的情況下。
在介紹完了具體的原理以後，下面結合R中提供的兩個軟體包來具體實踐一下置換檢驗。
2 用coin包做置換檢驗
包的安裝：
>install.packages(c("coin"))n> install.packages("lmPerm)n
幸運的是，lmPerm包已經可以更新，而且版本也更新到了2.1，目前是Marco Torchiano在維護。在此向原作者Bob Wheeler致意！
coin包提供了一個針對獨立性問題進行置換檢驗的一般性框架，主要回答以下問題：
a: 響應值與組的分配獨立嗎？
b: 兩個數值變數獨立嗎？
c: 兩個類型變數獨立嗎？
以下為可選置換檢驗的coin函數：
function_name(formula, data, distribution = )n
其中：
a: formula描述待檢驗變數間的關係，而且類別型變數和序數變數必須分別轉化為因子和有序因子
b: data是一個數據框，而且必須得是數據框
c: distribution指定經驗分布在零假設條件下的形式，有exact、asymptotic和approximate三個選項。
若distribution="exact"，那麼在零假設條件下，分布的計算是精確的（即依據所有可能的排列組合）。當然，也可以根據它的漸進分布（distribution="asymptotic"）或蒙特卡洛重抽樣（distribution="approxiamate(B=#)"）來做近似計算，其中#指所需重複的次數。distribution="exact"當前僅可用於兩樣本問題。
2.1 獨立兩樣本和K樣本檢驗
> library(coin)n載入需要的程輯包：survivaln> score <- c(40, 57, 45, 55, 58, 57, 64, 55, 62, 65)n> treatment <- factor(c(rep("A", 5),rep("B",5)))n> treatmentn [1] A A A A A B B B B BnLevels: A Bn> mydata <- data.frame(treatment, score)n> t.test(score~treatment,data = mydata, var.equal = TRUE)nntTwo Sample t-testnndata: score by treatmentnt = -2.345, df = 8, p-value = 0.04705nalternative hypothesis: true difference in means is not equal to 0n95 percent confidence interval:n -19.0405455 -0.1594545nsample estimates:nmean in group A mean in group B n 51.0 60.6 nn> oneway_test(score~treatment, data = mydata, distribution = "exact")nntExact Two-Sample Fisher-Pitman Permutation Testnndata: score by treatment (A, B)nZ = -1.9147, p-value = 0.07143nalternative hypothesis: true mu is not equal to 0nn> n
傳統t檢驗表明存在顯著性差異（p<0.05），而精確檢驗卻表明差異並不顯著（p>0.072）。
> library(MASS)n> UScrime <- transform(UScrime, So = factor(So))n> wilcox_test(Prob ~ So, data = UScrime, distribution = "exact")nntExact Wilcoxon-Mann-Whitney Testnndata: Prob by So (0, 1)nZ = -3.7493, p-value = 8.488e-05nalternative hypothesis: true mu is not equal to 0nn> wilcox.test(Prob ~ So, data = UScrime)nntWilcoxon rank sum testnndata: Prob by SonW = 81, p-value = 8.488e-05nalternative hypothesis: true location shift is not equal to 0n
此處wilcox_test()計算結果與第七章的wilcox.test()函數是一致的，這是因為exact也是wilcox.test()函數的默認演算法。
> library(multcomp)n載入需要的程輯包：mvtnormn載入需要的程輯包：survivaln載入需要的程輯包：TH.datan載入需要的程輯包：MASSnn載入程輯包：『MASS』nnThe following object is masked _by_ 『.GlobalEnv』:nn UScrimennn載入程輯包：『TH.data』nnThe following object is masked from 『package:MASS』:nn geysernn> set.seed(1234)n> library(coin)n> oneway_test(response ~ trt, data = cholesterol, distribution = approximate(B=9999))nntApproximative K-Sample Fisher-Pitman Permutation Testnndata: response bynt trt (1time, 2times, 4times, drugD, drugE)nchi-squared = 36.381, p-value < 2.2e-16nn> n
此處需要注意兩點，一個是此時的參考分布得自於數據9999次的置換，另外設置了隨機數種子，確保每次結果的一致性。
2.2 列聯表的獨立性
通過chisq_test()或cmh_test()函數，我們可用置換檢驗判斷兩類別型變數的獨立性。
> library(coin)n載入需要的程輯包：survivaln> library(vcd)n載入需要的程輯包：gridn> Arthritis <- transform(Arthritis, Improved = as.factor(as.numeric(Improved)))n> set.seed(134)n> set.seed(1234)n> chisq_test(Treatment ~ Improved, data = Arthritis, distribution = approximate(B=9999))nntApproximative Pearson Chi-Squared Testnndata: Treatment by Improved (1, 2, 3)nchi-squared = 13.055, p-value = 0.0018nn> n
經過9999次的置換，獲得一個近似的卡方檢驗。在這個例子當中，將Improved從一個有序因子變成一個分類因子，這是因為如果你用有序因子，則coin()將生成一個線性也線性趨勢檢驗，而不是卡方檢驗。Rstudio 提示如下：
> class(Arthritis$Improved)n[1] "factor"n> Arthritis <- transform(Arthritis, Improved = as.numeric(Improved))n> class(Arthritis$Improved)n[1] "numeric"n> chisq_test(Treatment ~ Improved, data = Arthritis, distribution = approximate(B=9999))nError in check(object) : n 『object』 does not represent a contingency problemn> n
2.3 數值變數間的獨立性
spearman_test()函數提供了兩數值變數的獨立性置換檢驗。
> states <- as.data.frame(state.x77)n> set.seed(1234)n> spearman_test(Illiteracy ~ Murder,data = state.x77, distribution = approximate(B=9999))nError in model.frame.default(formula = ~Illiteracy, data = c(3615, 365, : n data must be a data.frame, not a matrix or an arrayn> spearman_test(Illiteracy ~ Murder,data = state, distribution = approximate(B=9999))nError in terms.formula(formula, data = data) : object state not foundn> spearman_test(Illiteracy ~ Murder,data = states, distribution = approximate(B=9999))nntApproximative Spearman Correlation Testnndata: Illiteracy by MurdernZ = 4.7065, p-value < 2.2e-16nalternative hypothesis: true rho is not equal to 0nn> n
基於9999次的重複的近似檢驗可知，獨立性並不被滿足。另外，細心的讀者發現，我在輸入的時候犯了兩次錯誤，不難發現，在coin包中，可接受的數據類型是數據框。
2.4 兩樣本和K樣本相關性檢驗
> library(coin)n> library(MASS)n> wilcoxsign_test(U1~U2, data = UScrime, distribution = "exact")nntExact Wilcoxon-Pratt Signed-Rank Testnndata: y bynt x (pos, neg) nt stratified by blocknZ = 5.9691, p-value = 1.421e-14nalternative hypothesis: true mu is not equal to 0nn> n
當處於不同組的觀測已經被分配得當，或者使用了重複測量時，樣本相關檢驗便可派上用場。對於兩配對組的置換檢驗， 可使用wilcoxsign_test() 函數； 多於兩組時， 使用friedman_test()函數。
2.5 深入研究
coin包提供了一個置換檢驗的一般性框架，可以分析一組變數相對於其他任意變數，是否與第二組變數（可根據一個區組變數分層）相互獨立。而且independence_test()函數可以讓我們從置換角度來思考大部分傳統檢驗，進而在面對無法用傳統方法解決的問題時，使用戶可以自己構建新的統計檢驗。當然，這種靈活性也是有門檻的：要正確使用該函數必須具備豐富的統計知識。更多函數細節請參閱包附帶的文檔（運行vignette("coin")即可）。vignette("包的名稱")可以下載並用pdf閱讀軟體打開各個包的文檔。
下一節學習的lmPerm包，提供了線性模型（包括回歸和方差分析）的置換方法。
3 lmPerm包的置換檢驗
在lmPerm包中有lm()函數和aov()函數的修改版：lmp()和aovp()，進行置換檢驗，而非正態檢驗。參數也僅在原函數的基礎上添加了perm=參數。perm參數一共有三個選項，意義分別如下：
- Exact: 根據所有可能的排列組合生成精確檢驗；
- Prob: 從所有可能的排列中不斷抽樣，直至估計的標準差在估計的p值0.1之下，判停準則由可選的Ca參數控制；
- SPR: 使用貫序概率比檢驗來判斷何時停止抽樣。
為了在樣本數過大時減小計算機運行時間，精確建議僅在樣本數小於10才使用。因此，如果觀測數大於10，perm="Exact"將自動默認轉為perm="Prob"。
3.1 簡單回歸和多項式回歸
> library(lmPerm)n> set.seed(1234)n> fitPerm <- lmp(weight ~ height, data = women, perm = "Prob")n[1] "Settings: unique SS : numeric variables centered"n> summary(fitPerm)nnCall:nlmp(formula = weight ~ height, data = women, perm = "Prob")nnResiduals:n Min 1Q Median 3Q Max n-1.7333 -1.1333 -0.3833 0.7417 3.1167 nnCoefficients:n Estimate Iter Pr(Prob) nheight 3.45 5000 <2e-16 ***n---nSignif. codes: n0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnResidual standard error: 1.525 on 13 degrees of freedomnMultiple R-Squared: 0.991,tAdjusted R-squared: 0.9903 nF-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14 nn> n
> fitPerm <- lmp(weight~height + I(height^2), data=women, perm="Prob")n[1] "Settings: unique SS : numeric variables centered"n> summary(fitPerm)nnCall:nlmp(formula = weight ~ height + I(height^2), data = women, perm = "Prob")nnResiduals:n Min 1Q Median 3Q n-0.509405 -0.296105 -0.009405 0.286151 n Max n 0.597059 nnCoefficients:n Estimate Iter Pr(Prob) nheight -7.34832 5000 <2e-16 ***nI(height^2) 0.08306 5000 <2e-16 ***n---nSignif. codes: n0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnResidual standard error: 0.3841 on 12 degrees of freedomnMultiple R-Squared: 0.9995,tAdjusted R-squared: 0.9994 nF-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16 nn> n
以上分別為簡單回歸和多項式回歸的代碼，用置換檢驗來檢驗這些回歸是非常容易的，換一個僅多一個Perm參數的函數即可。從summary()函數的輸出可以看到增添的Iter列給出了達到判停標準所需要的迭代次數。
3.2 多元回歸
> library(lmPerm)n> set.seed(1234)n> states <- as.data.frame(state.x77)n> fitM <- lmp(Murder ~ Population + Illiteracy + Income + Frost, data = states, perm = "Prob")n[1] "Settings: unique SS : numeric variables centered"n> summary(fitM)nnCall:nlmp(formula = Murder ~ Population + Illiteracy + Income + Frost, n data = states, perm = "Prob")nnResiduals:n Min 1Q Median 3Q Max n-4.79597 -1.64946 -0.08112 1.48150 7.62104 nnCoefficients:n Estimate Iter Pr(Prob) nPopulation 2.237e-04 51 1.0000 nIlliteracy 4.143e+00 5000 0.0004 ***nIncome 6.442e-05 51 1.0000 nFrost 5.813e-04 51 0.8627 n---nSignif. codes: n0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnResidual standard error: 2.535 on 45 degrees of freedomnMultiple R-Squared: 0.567,tAdjusted R-squared: 0.5285 nF-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08 nn> n
與第八章的結果相比，Population不再顯著。當兩種方法的結果不一致時，需要更加謹慎的審視數據，可能是違反了相關前提假設或者存在離群點。
3.3 單因素方差分析和協方差分析
> library(lmPerm)n> library(multcomp)n載入需要的程輯包：mvtnormn載入需要的程輯包：survivaln載入需要的程輯包：TH.datan載入需要的程輯包：MASSnn載入程輯包：『MASS』nnThe following object is masked _by_ 『.GlobalEnv』:nn UScrimennn載入程輯包：『TH.data』nnThe following object is masked from 『package:MASS』:nn geysernn> set.seed(1234)n> fit <- aovp(response ~ trt, data = cholesterol, perm = "Prob")n[1] "Settings: unique SS "n> anova(fit)nAnalysis of Variance TablennResponse: responsen Df R Sum Sq R Mean Sq Iter Pr(Prob) ntrt 4 1351.37 337.84 5000 < 2.2e-16 ***nResiduals 45 468.75 10.42 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n> n
以上為9.1節中單因素方差分析問題，即各種療法對降低膽固醇的影響的置換檢驗，結果表明療法的效果不全相同。
> fit <- aovp(weight ~ gesttime + dose, data = litter, perm = "Prob")n[1] "Settings: unique SS : numeric variables centered"n> anova(fit)nAnalysis of Variance TablennResponse: weightn Df R Sum Sq R Mean Sq Iter Pr(Prob) ngesttime 1 161.49 161.493 5000 0.00060 ***ndose 3 137.12 45.708 4436 0.04531 * nResiduals 69 1151.27 16.685 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n
以上為第9章中的控制妊娠期時間相同時，觀測四種藥物劑量對鼠崽體重的影響問題，結果表明影響不太相同。
3.4 雙因素方差分析
> fit <- aovp(len~supp*dose, data = ToothGrowth, perm = "Prob")n[1] "Settings: unique SS : numeric variables centered"n> anova(fit)nAnalysis of Variance TablennResponse: lenn Df R Sum Sq R Mean Sq Iter Pr(Prob) nsupp 1 205.35 205.35 5000 < 2e-16 ***ndose 1 2224.30 2224.30 5000 < 2e-16 ***nsupp:dose 1 88.92 88.92 2276 0.04218 * nResiduals 56 933.63 16.67 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n> n
以上為第9章中的維生素C對豚鼠牙齒生長的影響問題，該問題有兩個可操作的因子：餵食方式（兩水平）和劑量（三水平）。
在不同的顯著性水平下，效應的影響水平也不盡相同。在0.05的顯著性水平下，三種效應都不等於0；在0.01的水平下，只有主效應顯著。
另外，當將aovp()應用到方差分析設計中時，它默認使用唯一平方和法（也稱為SAS類型Ⅲ平方和）。每種效應都會依據其他效應做相應調整。R中默認的參數化方差分析設計使於平衡設計，兩種方法結果相同，但是對於每個單元格觀測數不同的不平衡設計，兩種方法結果則不同。不平衡性越大，結果分歧越大。若在aovp()函數中設定seqs=TRUE，可以生成你想要的序貫平方和。
4 置換檢驗點評
R語言中還有其它可做置換檢驗的包。比如：
- perm包能夠實現coin包的部分功能，可以作為coin包所得結果的驗證；
- corrperm包提供了有重複測量的相關性的置換檢驗；
- logregperm包提供了Logistic回歸的置換檢驗；
- glmperm它涵蓋了廣義線性模型的置換檢驗。
從以上實例可以發現，基於正態理論的檢驗與上面置換檢驗的結果非常接近，這也說明正態理論適用於這些實例。另外一方面，置換檢驗也提供了另外一個十分強大可選檢驗思路。但是，當你的初始樣本對感興趣的總體情況代表性很差時，置換檢驗也無法提高推斷效果。
置換檢驗的優勢在於回答「效應是否存在」這也的問題，對於於獲取置信區間和估計測量精度是比較困難的。幸運的是，這正是接下來介紹的自助法大顯神通的地方。
5 自助法
所謂自助法，即從初始樣本重複隨機替換抽樣，生成一個或一系列待檢驗統計量的經驗分布。無需假設一個特定的理論分布，便可生成統計量的置信區間，並能檢驗統計假設。主要步驟如下：
- 從樣本中隨機選擇n個觀測，計算並記錄樣本均值，記為
，抽樣後將觀測放回，有些觀測可能被選擇多次，有些可能一直都不會被選中；
- 將步驟1重複N次，比如1000次，此時會有1000個樣本均值，記為
；
- 將所有的樣本均值
從下到大排序；
- 找出樣本均值2.5%和97.5%的分位點。此時即初始位置和最末位置的第25個數，它們就限定了95%的置信區間。
在潛在分布未知、樣本出現離群點、或者樣本量過小以及參數檢驗行不通的時候，自助法將是你用於生成置信區間和做假設檢驗的一個利器。
6 boot包中的自助法
自助法過程看起來非常複雜，但是對著例子走一遍你會明白，當然，首先還是需要安裝boot包。
install.packages("boot")n
通常來說，自助法有三個主要步驟：
- 寫一個能夠返回待研究統計量值的函數；
- 使用boot()函數對步驟1中產生的函數進行處理；
- 使用boot.ci()函數獲得步驟2生成的統計量的置信區間。
其中boot()函數的格式如下：
bootobject <- boot(data=, statistic=, R=, ...)n
生成了自助的樣本以後，可以通過print()函數和plot()函數來檢查結果。假如結果看起來合理，使用boot.ci()函數獲取統計量的置信區間。格式如下：
boot.ci(bootobject, conf = ,type = )n
參數如下：
以mtcars數據集為例，假設你想獲得95%的R平方值的置信區間。首先，寫一個獲取R平方值的函數：
> rsq <- function(formula, data, indices){n+ d <- data[indices,]n+ fit <- lm(formula, data = d)n+ return(summary(fit)$r.square)n+ }n
然後做大量的自助抽樣，比如1000次，代碼如下：
> library(boot)nn載入程輯包：『boot』nnThe following object is masked from 『package:survival』:nn amlnn> set.seed(1234)n> results <- boot(data = mtcars, statistic = rsq, R = 1000, formula= mpg~wt+disp)n> print(results)nnORDINARY NONPARAMETRIC BOOTSTRAPnnnCall:nboot(data = mtcars, statistic = rsq, R = 1000, formula = mpg ~ n wt + disp)nnnBootstrap Statistics :n original bias std. errornt1* 0.7809306 0.0133367 0.05068926n
也可將結果繪製出來：
> plot(results)n
> boot.ci(results, type = c("perc", "bca"))nBOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSnBased on 1000 bootstrap replicatesnnCALL : nboot.ci(boot.out = results, type = c("perc", "bca"))nnIntervals : nLevel Percentile BCa n95% ( 0.6838, 0.8833 ) ( 0.6344, 0.8549 ) nCalculations and Intervals on Original ScalenSome BCa intervals may be unstablen> n
生成置信區間的不同方法（type不同）將會獲得不同的區間。不管是那種方法，在本例中，0都不在置信區間內，因此拒絕零假設。
6.2 多個統計量的自助法
繼續前例，讓我們獲得一個統計量向量的（三個回歸係數（截距項、車重和發動機排量））的95%的置信區間。
首先，創建一個返回回歸係數向量的函數：
> bs <- function(formula, data, indices){n+ d <- data[indices,]n+ fit <- lm(formula, data = d)n+ return(coef(fit))n+ }n
然後，用該函數自助抽樣1000次：
> library(boot)n> set.seed(1234)n> results <- boot(data = mtcars, statistic = bs,n+ R = 1000, formula = mpg~wt+disp)n> print(results)nnORDINARY NONPARAMETRIC BOOTSTRAPnnnCall:nboot(data = mtcars, statistic = bs, R = 1000, formula = mpg ~ n wt + disp)nnnBootstrap Statistics :n original bias std. errornt1* 34.96055404 0.1378732942 2.485756673nt2* -3.35082533 -0.0539040965 1.170427539nt3* -0.01772474 -0.0001208075 0.008786803n
當對多個統計量自助抽樣時，在繪圖的時候需要添加一個索引參數，指明plot()和boot.ci()函數所分析的bootobject$t的列。比如，在本例中，索引1指截距項，索引2指車重，索引3指發動機排量。以下代碼繪製車重的結果：
plot(results, index = 2)n
> boot.ci(results,type = "bca", index = 1)nBOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSnBased on 1000 bootstrap replicatesnnCALL : nboot.ci(boot.out = results, type = "bca", index = 1)nnIntervals : nLevel BCa n95% (30.39, 40.25 ) nCalculations and Intervals on Original Scalen> boot.ci(results,type = "bca", index = 2)nBOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSnBased on 1000 bootstrap replicatesnnCALL : nboot.ci(boot.out = results, type = "bca", index = 2)nnIntervals : nLevel BCa n95% (-5.655, -1.194 ) nCalculations and Intervals on Original Scalen> boot.ci(results,type = "bca", index = 3)nBOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSnBased on 1000 bootstrap replicatesnnCALL : nboot.ci(boot.out = results, type = "bca", index = 3)nnIntervals : nLevel BCa n95% (-0.0331, 0.0010 ) nCalculations and Intervals on Original Scalen> n
有個問題，我本來想用以下代碼的，但是沒有輸出，也沒有提示，不知道是為何？
for(i in 1:3){n boot.ci(results,type = "bca", index = i)n }n
最後我們討論下，在使用自助法時經常被提起的兩個問題：
- 初始樣本需要多大？
- 應該重複多少次？
對於第一個問題，我們無法給出簡單的回答。有些人認為只要樣本能夠較好地代表總體，初始樣本大小為20～30即可得到足夠好的結果。從感興趣的總體中隨機抽樣的方法可信度最高，它能夠保證初始樣本的代表性。對於第二個問題，我發現1000次重複在大部分情況下都可滿足要求。由於計算機資源變得廉價，如果你願意，也可以增加重複的次數。
7 小結
本章我們介紹了一系列基於隨機化和重抽樣的計算機密集型方法，它們使你無需理論分布的知識便能夠進行假設檢驗，獲得置信區間。當數據來自未知分布，或者存在嚴重的離群點，或者樣本量過小，又或者沒有參數方法可以回答你感興趣的假設問題時，這些方法是非常實用的。
置換檢驗和自助法並不是萬能的，它們無法將爛數據轉化為好數據。如果初始樣本對於總體情況的代表性不佳，或者樣本量過小而無法準確地反映總體情況，這些方法也是愛莫能助。
推薦閱讀：
※《R語言實戰》第四部分第十五章-時間序列學習筆記（II）
※Python數據分析及可視化實例之Bokeh與Jupyter生成可視化圖表（8）
※市場規模是怎麼估算出來的？都有哪些依據？
※老鹿玩數據——不光是求婚神器（二）
※學會處理你的視頻數據，做粉絲更喜歡的內容