R之重抽樣與自助法（一）

04-18

先來回顧一下幾種統計檢驗的假設條件：

（1）t檢驗的假設條件是結果變數為連續型的組間比較，且呈正態分布。

（2）回歸OLS(普通最小二乘, ordinary least square)的假設條件是結果變數呈正態性，獨立性，同方差性。

（3）方差分析假設：

每組樣本取自一個正太分布的總體（正態性）
所有總體具有相同的方差（同方差性）
每組樣本中的觀測值是隨機選取的，且相互獨立
因子效應是可加的

離群點會導致問題。另外，我們也需要保證F檢驗表現良好。特別的，即使樣本不滿足正太分布，F統計檢驗也是魯棒的：

樣本是對稱的且單峰的
每組樣本容量相同且大於10個

好了，我們看到上面不論 t 檢驗、OSL回歸還是方差分析ANOVA，都要求觀測數據抽樣自正態分布或者其他性質較好的理論分布，但在實際情況中統計假設不一定滿足，比如數據抽樣自未知或混合分布、樣本量過小、存在離群點、基於理論分布設計合適的統計檢驗過於複雜且數學上難以處理等情況，這是就可以用基於隨機化和重抽樣的統計方法，這裡介紹重抽樣與自助法。

置換檢驗的步驟（如果兩組數據的均值相等，那麼分配給觀測得分的標籤（A或者B）就應當是任意的）：

（1）與參數方法類似，計算觀測數據的 t 統計量，稱為t0；

（2）將所有觀測值放在一起，隨機分配一般給A, 一般給B；
（3）計算並記錄新觀測的 t 統計量；
（4）重複（3）~（4），假設樣本有n個，那一共有 $C_{n}^{[n/2]}$ 種可能性，有這麼多 t 值(這裡也看到了為什麼叫Permutation Test的原因，利用了組合數）；
（5）將這 $C_{n}^{[n/2]}$ 個 t 統計量按升序排列，這便是基於樣本數據的經驗分布;
（6）如果t0 落在經驗分布的95%部分的外面，則在0.05水平

下面這個例子雖然不是用R語言來寫的，但是可以幫助理解置換檢驗:

Permutation tests

——by Saunak Sen

高血壓數據示例
我們測量了兩組小鼠的收縮壓。為了簡單起見，我們就看D4Mit214標記基因分型的50個（隨機選擇的）小鼠。我們想檢測D4Mit214標記基因型與血壓之間的關係。下表以增序顯示了標記基因型BA（雜合）或BB（純合子）小鼠的收縮壓（mm Hg）。
BA: 86 88 89 89 92 93 94 94 94 95 95 96 96 97 97 98 98 99 99 101 106 107 110 113 116 118BB: 89 90 92 93 93 96 99 99 99 102 103 104 105 106 106 107 108 108 110 110 112 114 116 116
粗略地一看好像BB組的小鼠血壓更高。我們可以使用t檢驗來檢驗兩個基因型組中的平均血壓相等的零假設。
首先在Stata中讀入數據：
/* read data */. insheet using hyper.csv(22 vars, 50 obs)
進行 t 檢驗：
/* perform t-test allowing for unequal variances */. ttest bp, by(d4mit214) unequalTwo-sample t test with unequal variances------------------------------------------------------------------------------ Group | Obs Mean Std. Err. Std. Dev. [95% Conf. Interval]---------+-------------------------------------------------------------------- BA | 26 98.48077 1.657301 8.45061 95.06749 101.894 BB | 24 103.1833 1.637388 8.021529 99.79614 106.5705---------+--------------------------------------------------------------------combined | 50 100.738 1.202249 8.501185 98.32199 103.154---------+-------------------------------------------------------------------- diff | -4.702564 2.329739 -9.386925 -.0182032------------------------------------------------------------------------------ diff = mean(BA) - mean(BB) t = -2.0185Ho: diff = 0 Satterthwaites degrees of freedom = 47.958 Ha: diff < 0 Ha: diff != 0 Ha: diff > 0 Pr(T < t) = 0.0246 Pr(|T| > |t|) = 0.0492 Pr(T > t) = 0.9754
t 檢驗表明，BA基因型組的平均血壓比BB基因型組低4.7個單位。雙尾檢驗的p值為0.049 - 剛剛達到「顯著」水平。
但是注意，t 檢驗假設 (a) 觀測值之間獨立 (b) 數據正態分布。如果分布是對稱單峰且連續的，那檢驗是魯棒的。我們沒有先驗經驗來確定血壓滿足正態分布，而且樣本量太小，無法可靠地檢查假設（低功率）。而且在這裡數據可能不呈現對稱單峰的分布，所以我們用置換檢驗：

「洗」數據

如果我們像洗牌那樣隨機將血壓分配給基因型組，也就是將一般的血壓值分配給a，另外一半給b，這樣就可以計算新的 t 統計量。這樣的分配一共有 $C_{50}^{25}$ 種（也就是為什麼叫置換檢驗的原因了），將這 $C_{50}^{25}$ 個 t 統計量按升序排列，我們就得到了基於樣本數據的經驗分布。因此，從重複隨機洗數據獲得的 t 檢驗統計量的分布可以用於近似檢驗統計量的零分布。

我們列舉其中一種情形 (將原始BA組中的值標記為a，將BB組中的值標記為b):
BA: 101a 104b 98a 94a 92a 99b 116b 103b 86a 90b 88a 93b 112b 106b 107b 97a 93a 113a 110a 118a 99b 95a 96a 99a 110b 89aBB: 105b 108b 116a 94a 106b 95a 94a 106a 108b 114b 99b 89b 99a 96b 102b 110b 97a 98a 93b 96a 89a 92b 116b 107a
BA和BB基因型組的平均血壓分別為100.3和101.2。現在的差值相對於實際觀察到的差值更小，為（0.9）-4.7。現在我們可以提出這樣一個問題：如果血壓和基因型之間沒有關聯，血壓的差值為多少？與樣品中的差值相比如何？為了回答這個問題，我們重複洗牌，記錄兩個基因型組之間的差值。在400次洗牌之後，BA和BB組中的手段之間的差值分布如下：
/* make a new variable that is 1 if d4mit214 is BA and 0 otherwise */generate geno=0 if d4mit214=="BB"replace geno=1 if d4mit214=="BA"/* we can get the difference between the means in the two genotype groups using the regress command */. regress bp geno Source | SS df MS Number of obs = 50-------------+------------------------------ F( 1, 48) = 4.06 Model | 275.984052 1 275.984052 Prob > F = 0.0496 Residual | 3265.25335 48 68.0261114 R-squared = 0.0779-------------+------------------------------ Adj R-squared = 0.0587 Total | 3541.2374 49 72.270151 Root MSE = 8.2478------------------------------------------------------------------------------ bp | Coef. Std. Err. t P>|t| [95% Conf. Interval]-------------+---------------------------------------------------------------- geno | -4.702564 2.334697 -2.01 0.050 -9.396787 -.0083409 _cons | 103.1833 1.683574 61.29 0.000 99.79828 106.5684------------------------------------------------------------------------------. permute bp "regress bp geno" _b, saving(perm400) reps(400). clear. use perm400. stem b_geno, line(1)Stem-and-leaf plot for b_geno (_b[geno])b_geno rounded to nearest multiple of .1plot in units of .1 -7* | 7 -6* | 0 -5* | 533 -4* | 988776643332200 -3* | 777776655555544444332221111000 -2* | 99887775544443332222211111111100000 -1* | 999987777766655554444433332222221111100000000 -0* | 99999998888888888877777777666555555444333332222222111 0* | 0000000000011112222222233333334444445555555555667777788888899999999999 1* | 0000011111112222222333334444445555555556666667777888999999 2* | 0000000122222233344445555556666667777889999 3* | 0000001123333344566679 4* | 000112223333566789 5* | 1456 6* | 35
經過更多次之後（十萬次洗牌），我們從Stata得到這個。
. permute bp "regress bp geno" _b, reps(100000)command: regress bp genostatistics: b_geno = _b[geno] b_cons = _b[_cons]permute var: bpMonte Carlo permutation statistics Number of obs = 50 Replications = 100000------------------------------------------------------------------------------T | T(obs) c n p=c/n SE(p) [95% Conf. Interval]-------------+----------------------------------------------------------------b_geno | -4.702564 4986 1.0e+05 0.0499 0.0007 .0485191 .0512271b_cons | 103.1833 2485 1.0e+05 0.0249 0.0005 .023894 .0258337------------------------------------------------------------------------------Note: confidence intervals are with respect to p=c/nNote: c = #{|T| >= |T(obs)|}
由上可知，洗完後大約4.99％的差值比觀察所得的差值大-4.7。這意味著p值為0.0499，但注意到這個p值存在置信區間。為什麼？因為這個p值是根據隨機洗牌來估計的。一般來說，p值越小，準確估算所需要的洗牌次數越多（事件越少，查找時間越長）。
本次案例中，置換測試的p值與從t檢驗獲得的p值幾乎相同。然而，這不一定是這種情況（這也表明了t檢驗的魯棒性）。

R語言中常用的有coin包和Imperm包，coin包用來解決響應值與組的分配的獨立性檢驗、兩個數值變數的獨立性判斷和兩個類別型變數的獨立性判斷，而Imperm包提供了線性模型的置換方法，可用於回歸和方差分析。

（一）獨立性及相關性檢驗（coin包）

> library(coin)載入需要的程輯包：survival> score<- c(40,57,45,55,58,57,64,55,62,65)> treatment <- c(rep("a",5),rep("b",5))> mydata <- data.frame(treatment, score)> mydata treatment score1 a 402 a 573 a 454 a 555 a 586 b 577 b 648 b 559 b 6210 b 65> t.test(scsore ~ treatment, data = mydata, var.equal = TRUE)Error in eval(expr, envir, enclos) : 找不到對象scsore> t.test(score ~ treatment, data = mydata, var.equal = TRUE) Two Sample t-testdata: score by treatmentt = -2.345, df = 8, p-value = 0.04705alternative hypothesis: true difference in means is not equal to 095 percent confidence interval: -19.0405455 -0.1594545sample estimates:mean in group a mean in group b 51.0 60.6 > oneway_test(score ~ treatment, data = mydata, distribution = "exact") Exact Two-Sample Fisher-Pitman Permutation Testdata: score by treatment (a, b)Z = -1.9147, p-value = 0.07143alternative hypothesis: true mu is not equal to 0

exact表示分布是精確的，即依據所有可能的排列組合。傳統 t 檢驗表明存在顯著性差異（p<0.05），而精確檢驗卻顯示差異不顯著。由於數據只有10個，置換檢驗的結果相對更為可信。我們再用下秩和檢驗：

> wilcox_test(score~ treatment, data= mydata) Asymptotic Wilcoxon-Mann-Whitney Testdata: score by treatment (a, b)Z = -1.7865, p-value = 0.07403alternative hypothesis: true mu is not equal to 0

列表聯中的獨立性：

> chisq_test(Treatment~Improved, data = Arthritis, daistribution = approximate(B=9999)) Asymptotic Pearson Chi-Squared Testdata: Treatment by Improved (1, 2, 3)chi-squared = 13.055, df = 2, p-value = 0.001463

數值變數間的獨立性：

> spearman_test(Illiteracy ~ Murder, data = states, distribution = approximate(B=9999)) Approximative Spearman Correlation Testdata: Illiteracy by MurderZ = 4.7065, p-value < 2.2e-16alternative hypothesis: true rho is not equal to 0

兩樣本和K樣本相關性檢驗：

> wilcoxsign_test(U1~U2, data= UScrime, distribution = "exact") Exact Wilcoxon-Pratt Signed-Rank Testdata: y by x (pos, neg) stratified by blockZ = 5.9691, p-value = 1.421e-14alternative hypothesis: true mu is not equal to 0

（二）回歸（lmPerm包）

lmp()和avop()函數的參數與lm()和aov()函數類似，只額外添加了perm = 參數。perm可以選擇的值有 Exact , Prob , SPR。Exact根據所有可能的排列組合生成精確檢驗，但是只適用於小樣本，注意，若觀測值大於10， perm = "Exact"將自動默認轉為「Prob"。Prob不斷從所有可能的排列中不斷抽樣，直到估計的標準差在估計的p值0.1之下，判停準則由可選的Ca參數控制。SPR使用貫序概率比檢驗來判斷何時停止抽樣。

> library(lmPerm)> set.seed(1234)> fit <- lmp(Murder ~ Population + Illiteracy + Income + Frost, data = states, perm = "Prob")[1] "Settings: unique SS : numeric variables centered"> summary(fit)Call:lmp(formula = Murder ~ Population + Illiteracy + Income + Frost, data = states, perm = "Prob")Residuals: Min 1Q Median 3Q Max -4.79597 -1.64946 -0.08112 1.48150 7.62104 Coefficients: Estimate Iter Pr(Prob) Population 2.237e-04 51 1.0000 Illiteracy 4.143e+00 5000 0.0004 ***Income 6.442e-05 51 1.0000 Frost 5.813e-04 51 0.8627 ---Signif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『』 1Residual standard error: 2.535 on 45 degrees of freedomMultiple R-Squared: 0.567, Adjusted R-squared: 0.5285 F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08

如果用lm回歸模型Population和Illiteracy均顯著，而在置換檢驗中Population不再顯著。當兩種方法所得結果不一致時，你需要更加謹慎地審視數據，這很可能因為違反了正態性假設或者存在離群點。

（三）方差分析（lmPerm包）

> nrow(ToothGrowth)[1] 60> set.seed(1234)> fit <- aovp( len ~ supp*dose, data = ToothGrowth, perm = "Prob")[1] "Settings: unique SS : numeric variables centered"> anova(fit)Analysis of Variance TableResponse: len Df R Sum Sq R Mean Sq Iter Pr(Prob) supp 1 205.35 205.35 5000 < 2e-16 ***dose 1 2224.30 2224.30 5000 < 2e-16 ***supp:dose 1 88.92 88.92 2032 0.04724 * Residuals 56 933.63 16.67 ---Signif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『』 1

在0.05水平下，三種效應都不等於0；在0.001水平下，只有主效應顯著。

置換檢驗主要用於生成檢驗的零假設的p值，它有助於回答「效應是否存在」這樣的問題。不過，置換方法對於獲取置信區間和估計測量精度是比較困難的。這是就可以用自助法來獲取置信區間和估計測量精度。

參考自：R語言實戰