|
本節(jié)主要總結(jié)「數(shù)據(jù)分析」的「Resampling」重抽樣思想,，并通過 R 語言實(shí)現(xiàn),。
有一種東西叫作「?jìng)鹘y(tǒng)」,，它在很多時(shí)候很有用,，但會(huì)讓你思維固化，在新的環(huán)境下讓你出錯(cuò),。
在總結(jié)回歸分析和方差分析的時(shí)候 ④Ｒ語言之?dāng)?shù)據(jù)分析「初章」,，我總是會(huì)在模型的建立之前提到「統(tǒng)計(jì)假設(shè)」，在模型建立之后進(jìn)行「假設(shè)檢驗(yàn)」,，原因想必大家都能理解,，就是因?yàn)檫@些「統(tǒng)計(jì)假設(shè)」是我們模型建立思想的基礎(chǔ)，是支撐我們模型正確性的「必要條件」。但是,，不可否認(rèn)的是,，這些「必要條件」最終會(huì)成為我們「數(shù)據(jù)分析」的局限，讓我們對(duì)「不滿足條件的數(shù)據(jù)集」束手無策,。 本節(jié),，我們就來解決這個(gè)「必要條件」中的其中一條假設(shè)，從特例到普遍,。 統(tǒng)計(jì)假設(shè)中有一條,，叫做「假定觀測(cè)數(shù)據(jù)抽樣自正態(tài)分布或者是是其他性質(zhì)較好的分布」，那么當(dāng)數(shù)據(jù)集抽樣自「未知分布,、混合分布」,、樣本容量過小或存在離群點(diǎn)時(shí)，傳統(tǒng)的統(tǒng)計(jì)方法所得到的模型可能就會(huì)不那么準(zhǔn)確,，原因之前已經(jīng)講過,，這個(gè)時(shí)候「Resampling」 的思想就出現(xiàn)了。它拋棄了分布的理論,，而是完全基于同一個(gè)樣本,，在這個(gè)樣本中多次重復(fù)抽樣，然后將所有抽樣的結(jié)果作為總體,，將原樣本放到其中去檢驗(yàn),，判定其顯著性。因?yàn)樾枰啻沃貜?fù)抽樣,，所有它被稱為重抽樣「Resampling」,。
1. Resampling 的分類1.1. 置換檢驗(yàn)（ permutation test ） 這個(gè)方法是傳統(tǒng)統(tǒng)計(jì)方法的創(chuàng)建者 R. A. Fisher 建立的,，但是由于這個(gè)方法的計(jì)算量過大,、且計(jì)算機(jī)技術(shù)也未成熟，他最后放棄了這個(gè)方法,。但是,，數(shù)十年后的今天，計(jì)算機(jī)技術(shù)的高速發(fā)展,，這個(gè)方法終于能夠?qū)崿F(xiàn)并發(fā)揮其價(jià)值,。 為了更清楚的說明置換檢驗(yàn)的思想，我舉一個(gè)「兩總體均值之差」的推斷問題： 現(xiàn)在有兩套學(xué)習(xí)方案 A 和 B,，在 10 個(gè)受試者中隨機(jī)抽取 5 個(gè)按照方案 A 學(xué)習(xí),，另外 5 個(gè)按照方案 B 學(xué)習(xí)，在學(xué)習(xí)完畢后對(duì) 10 個(gè)受試者進(jìn)行測(cè)試,，得到分?jǐn)?shù)如下： 方案 A 方案 B 91 89 88 78 76 93 79 81 82 77原假設(shè)H_{0} ：兩種方案的總體均值相等 ,；備擇假設(shè)H_{a}：兩種方案的總體均值不等
這種問題，用參數(shù)方法來解決你應(yīng)該是很熟悉的,，現(xiàn)在我來解釋置換檢驗(yàn)的思想是怎樣的,。置換檢驗(yàn)的思想下,，由于我們?cè)僭O(shè)兩個(gè)方案的總體均值是相等的，那么如果兩種方案真的是等價(jià)的話,，則 10 個(gè)受試者的分?jǐn)?shù)的分配應(yīng)該是任意的,。這么講你可能還是不太理解，其實(shí)在原假設(shè)成立的前提下,，我們把兩種方案的總體視為等價(jià)的,，那么我們同樣也可將這兩個(gè)方案的總體視為同一個(gè)總體，這樣的話,，10 個(gè)受試者者的分?jǐn)?shù)在可以在任意位置,，而不用在乎這個(gè)位置是屬于方案 A 還是方案 B。在這樣的前提下,，我們利用重抽樣得到一個(gè)「源于已有樣本的抽樣分布」,，再利用 t 統(tǒng)計(jì)量在這個(gè)「源于已有樣本的抽樣分布」上判別初始觀測(cè)數(shù)據(jù)的顯著性，具體步驟如下：
1. 計(jì)算初始觀測(cè)數(shù)據(jù)的 t 統(tǒng)計(jì)量,，記為 t0
2. 將初始觀測(cè)數(shù)據(jù)的 10 個(gè)得分,，作為一個(gè)「重抽樣的總體」，從中抽取 5 個(gè)放到 A 組,，另五個(gè)放到 B 組,，并計(jì)算此種情況下的 t 統(tǒng)計(jì)量
3. 重復(fù)步驟 2，直到窮盡所有可能,。其實(shí)就是一個(gè)組合問題,，一共有252種可能
4. 將 252 種可能的 t 統(tǒng)計(jì)量，按照從小到大排序,，得到了「源于已有樣本的抽樣分布」
5. 根據(jù)t0在此分布中的位置,，和我們確定的顯著性水平，判別原假設(shè)我們可否拒絕 1.2. 交叉驗(yàn)證（ Cross validation ） 在交叉驗(yàn)證的思想中,，一個(gè)樣本被隨機(jī)的分成兩個(gè)或 K 個(gè)子樣本,，將其中一個(gè)作為驗(yàn)證集，其他 K-1 個(gè)子樣本組合成訓(xùn)練集,，在訓(xùn)練集上獲得回歸方程,，而后在驗(yàn)證集上做出預(yù)測(cè)，這為一重驗(yàn)證,，而后所有子樣本輪流充當(dāng)驗(yàn)證集,，如此往復(fù)，直至遍歷所有可能,。這樣我們就可以得到 K 個(gè)“預(yù)測(cè)誤差”,，求其平均值即為所謂的“ CV 值”，所以常說的 CV 值實(shí)際上是預(yù)測(cè)誤差期望Err的一個(gè)估計(jì)值。 交叉驗(yàn)證的思想比較容易理解,，也是一種十分符合我們直覺的驗(yàn)證方法,，這里就不舉例說明了。
但值得注意的一點(diǎn)是,，K 的取值會(huì)影響我們預(yù)測(cè)的準(zhǔn)確性,，具體是怎么影響的涉及到“偏誤”、“方差”以及“計(jì)算效率”的權(quán)衡,，主要是“偏誤”與“計(jì)算效率”的權(quán)衡,，這里就不深究，感興趣的同學(xué)可以看看 「模型選擇的一些基本思想和方法」這篇文章,。
不過,，我們可以粗略的設(shè)置 K 的取值：
大樣本，5 重交叉驗(yàn)證,，對(duì)預(yù)測(cè)誤差的估計(jì)已經(jīng)足夠準(zhǔn)確,，這里偏向于增大計(jì)算效率 小樣本，10重交叉驗(yàn)證,，為了得到相對(duì)準(zhǔn)確的估計(jì),，這里考慮損失計(jì)算效率
1.3. 刀切法（ Jackknife ） 刀切法的思想，十分適用于分布的分散過寬或者出現(xiàn)異常值的數(shù)據(jù),，通俗來講就是「既然我們所擁有的數(shù)據(jù)都只是樣本,，都是抽樣得到的，那么我們?cè)谧鐾茢嗪凸烙?jì)的時(shí)候,，扔掉幾個(gè)觀測(cè),，看看效果如何?！乖诘肚蟹ㄖ形覀兎磸?fù)的行為就是從樣本中選取一個(gè)觀測(cè),，把它拋棄，而后利用剩余的數(shù)據(jù)做出估計(jì),，如此往復(fù),，直到每一個(gè)觀測(cè)都被拋棄過一次,。根據(jù)每次估計(jì)的統(tǒng)計(jì)量,，取均值，而后將這個(gè)均值與原樣本（沒有觀測(cè)沒拋棄）上估計(jì)到到的統(tǒng)計(jì)量想比較,，就可以得到「偏差校正刀切估計(jì)值」,。我的表述可能你能聽懂，但是還理解不深,，下面我舉例說明,。 因?yàn)橐粋€(gè)理科生的各科成績之間可能存在某種相關(guān)性，現(xiàn)在我們想用用學(xué)生的物理成績和化學(xué)成績來預(yù)測(cè)他的數(shù)學(xué)成績，現(xiàn)在有 50 個(gè)觀測(cè)： 序號(hào) 物理 化學(xué) 數(shù)學(xué) 1 78 87 90 2 88 79 83 3 67 71 78 … … … … 50 89 88 92
1.4. Bootstrap 與刀切法相比,，Bootstrap 法的重抽樣思想更加透徹,。刀切法的重抽樣次數(shù)與觀測(cè)數(shù)量有關(guān)，而 Bootstrap 法的抽樣次數(shù)理論上是沒有上限的,，將原樣本在計(jì)算機(jī)性能允許的基礎(chǔ)上盡可能重復(fù)抽樣更多次,，得到了一個(gè)「virtual population 」，從其中我們可以得到統(tǒng)計(jì)量的「虛擬潛在分布」,，并根據(jù)此分布估計(jì)置信區(qū)間,。 舉例說明： 現(xiàn)在有一個(gè)容量為 10 的樣本，均值 20,，標(biāo)準(zhǔn)差 2,，希望根據(jù)樣本估計(jì)總體均值的置信區(qū)間。 按照傳統(tǒng)思想,，我們可以假設(shè)樣本分布為正態(tài)分布,，而后根據(jù)正態(tài)分布計(jì)算總體均值的置信區(qū)間。但當(dāng)樣本分布不服從正態(tài)分布時(shí),，可以利用 Bootstrap 法：
1. 有放回地從樣本中抽取 10 個(gè)觀測(cè),，并計(jì)算其均值（注意由于抽樣是有放回的，重抽樣得到的樣本與原樣本雖然容量相同,，但是觀測(cè)不一定相同,，因?yàn)橛械挠^測(cè)可能被多次選擇，有的可能一直都不會(huì)被選中）
2. 重復(fù) 1,，次數(shù)為 2000
3. 將 2000 個(gè)均值升序排列
4. 確定顯著性水平,，計(jì)算得到置信區(qū)間（例如：要求 95% 的置信區(qū)間，找出 2.5% 和 97.5% 的分位點(diǎn),，其中間部分就為 95% 的置信區(qū)間）
2. R 語言實(shí)現(xiàn)重抽樣2.1. 置換檢驗(yàn) 格式： > function_name( formula , data , distribution = )# formula 是被檢驗(yàn)變量之間的關(guān)系,， data 是數(shù)據(jù)框，distribution 指定重抽樣的方式
這里的重抽樣方式,，除了之前將的按照完全排列組合抽樣之外（ distribution = 'exact' ）,，為了節(jié)省計(jì)算機(jī)資源，可以選擇根據(jù)漸進(jìn)分布抽樣（ distribution = 'asymptotic' ）,，和蒙特卡洛重抽樣（ distribution = 'approxiamate(B=#)' ）# 為重復(fù)次數(shù),，一般要設(shè)定隨機(jī)種子，以便將來重現(xiàn),。 2.1.1. 兩樣本和 K 樣本的獨(dú)立性檢驗(yàn)（ coin 包 ） 兩樣本獨(dú)立性的參數(shù)檢驗(yàn)法（ 置換檢驗(yàn) ） > oneway_test(score~treatment,data=mydata,distribution='exact')# 這里表達(dá)式 ~ 的兩邊分別為應(yīng)變量和分組變量# 傳統(tǒng)為 t.test() 函數(shù) 兩樣本獨(dú)立性的非參數(shù)檢驗(yàn)法（ 置換檢驗(yàn) ） > wilcox_test(Prob ~ So ,data=UScrime, distribution='exact')# 這里表達(dá)式 ~ 的兩邊分別為應(yīng)變量和分組變量# 傳統(tǒng)為 wilcox.test() 函數(shù)
K 樣本獨(dú)立性的參數(shù)檢驗(yàn)法（ 置換檢驗(yàn) ） > oneway_test(response ~ trt ,data=cholesterol,distribution='approxiamate(B=9999)')# 這里表達(dá)式 ~ 的兩邊分別為應(yīng)變量和分組變量# 傳統(tǒng)為 單因素方差分析
2.1.2. 列聯(lián)表中兩分組變量之間的獨(dú)立性檢驗(yàn)（ coin 包 ） > Arthritis <- transform(arthritis,="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" improved="as.factor(as.numeric(Improved)))"> set.seed(1234)> chisq_test(Treatment~Improved,data=Arthritis, distribution=approximate(B=9999))# 這里表達(dá)式 ~ 的兩邊分別為兩個(gè)分組變量# 如果使用有序因子,，則次函數(shù)為線性趨勢(shì)檢驗(yàn)，因此第一行代碼將 Improved 轉(zhuǎn)換成無序
2.1.3. 數(shù)值型變量之間的獨(dú)立性檢驗(yàn)（ coin 包 ） > states <- as.data.frame(state.x77)=""> set.seed(1234)> spearman_test(Illiteracy ~ Murder , data=states, distribution=approximate(B=9999))# 這里表達(dá)式 ~ 的兩邊分別為兩個(gè)數(shù)值型變量# coin 包中必須為數(shù)據(jù)框
2.1.4. 兩樣本和 K 樣本的相關(guān)性檢驗(yàn)（ coin包 ） 兩樣本 > wilcoxsign_test(U1~U2,data=UScrime,distribution='exact')
K 樣本 > friedman_test(times ~ methods | block, data = rounding)
2.1.5. 簡(jiǎn)單回歸和多項(xiàng)式回歸（ lmPerm 包 ） 簡(jiǎn)單線性回歸（置換檢驗(yàn)） > set.seed(1234)> fit1 <- lmp(weight~height,data="women," perm='Prob' )=""> summary(fit)# lmp() 函數(shù)與 lm() 函數(shù)類似,，只是多了 perm 參數(shù),，其可選值有 Exact,、Prob 或SPR# Exact 根據(jù)所有可能的排列組合生成精確檢驗(yàn)。Prob 從所有可能的排列中不斷抽樣,，直至估計(jì)的標(biāo)準(zhǔn)差在估計(jì)的 p 值 0.1 之下,，判停準(zhǔn)則由可選的 Ca 參數(shù)控制。SPR 使用貫序概率比檢驗(yàn)來判斷何時(shí)停止抽樣,。# 觀測(cè)數(shù)大于10,，perm='Exact'將 自動(dòng)默認(rèn)轉(zhuǎn)為 perm='Prob'，因?yàn)榫_檢驗(yàn)只適用于小樣本問題,。
多項(xiàng)式回歸（置換檢驗(yàn)） > set.seed(1234)> fit2 <- lmp(weight~height+i(height^2),data="women,perm='Prob')"> summary(fit)
2.1.6. 多元回歸（置換檢驗(yàn)） > 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. 單因素方差分析和協(xié)方差分析 單因素方差分析（置換檢驗(yàn)） > set.seed(1234)> fit <- aovp(response="" ~="" trt,="" data="cholesterol," perm='Prob' )=""> summary(fit)# aovp() 函數(shù)與 aov() 函數(shù)類似,，只是多了 perm 參數(shù)，其可選值有 Exact,、Prob 或 SPR
單因素協(xié)方差分析（置換檢驗(yàn)） > 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. 交叉驗(yàn)證 《 R 語言實(shí)戰(zhàn) 》中的自編函數(shù) shrinkage ( ) 可以實(shí)現(xiàn) K 重交叉驗(yàn)證,，默認(rèn)為 10 重，可以通過參數(shù) 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 語言實(shí)踐 》中并未介紹刀切法的實(shí)現(xiàn),，可以參考以下代碼 # 抽樣得到樣本> set.seed(1234)> x <- rnorm(100,0,5)=""> # 運(yùn)用刀切法重抽樣> 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)+="" }=""> # 計(jì)算刀切法得到的抽樣方差> mean(jack(x))/length(x)[1] 24.97106> # 計(jì)算抽樣方差> var(x)[1] 25.22075
2.4. Bootstrap 格式： # Bootstrap 重抽樣，并返回統(tǒng)計(jì)量> bootobject <- boot(data="," statistic="," r="," ...)#="" data="" 為原樣本,，也是我們重抽樣的數(shù)據(jù)來源,，可以是向量、矩陣或者數(shù)據(jù)框#="" statistic="" 為生成="" k="" 個(gè)統(tǒng)計(jì)量的函數(shù),，這個(gè)函數(shù)的參數(shù)中必須包括="" indices,，indices="" 參數(shù)用于="" boot()="" 對(duì)原樣本重抽樣#="" r="" 為重抽樣的次數(shù)#="" 獲取置信區(qū)間=""> boot.ci(bootobject, conf=, type= )# bootobject 為存儲(chǔ)重抽樣結(jié)果的對(duì)象# conf 為置信區(qū)間的置信度# type 為置信區(qū)間的類型
2.4.1. 對(duì)單個(gè)統(tǒng)計(jì)量使用自助法 > #獲取 R 平方值的自編函數(shù)> 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 對(duì)象> print(results)> plot(results)> #計(jì)算置信區(qū)間> boot.ci(results, type=c('perc','bca'))# type 參數(shù)為 perc 時(shí)展示的是樣本均值，為 bca 時(shí)將根據(jù)偏差對(duì)區(qū)間做簡(jiǎn)單調(diào)整
2.4.2. 多個(gè)統(tǒng)計(jì)量的自助法 > #返回回歸系數(shù)向量的自編函數(shù)> 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 對(duì)象> print(results)> plot(results,index=2)> #計(jì)算置信區(qū)間> boot.ci(results,type = 'bca', index=2)> boot.ci(results,type = 'bca', index=3)# index 參數(shù)用于指明需要分析的重抽樣所得樣本的列
|
