對于正態(tài)分布或其他已知分布的數(shù)據(jù)喉脖，有相應(yīng)的假設(shè)檢驗(yàn)與置信區(qū)間的計(jì)算方法况鸣，但是當(dāng)數(shù)據(jù)抽樣自未知或混合分布何乎、樣本量過小昔搂、存在離群點(diǎn)玲销、基于理論分布設(shè)計(jì)合適的統(tǒng)計(jì)檢驗(yàn)過于復(fù)雜且數(shù)學(xué)上難以處理等情況，就需要使用基于隨機(jī)化和重抽樣的統(tǒng)計(jì)方法摘符。
上次推文主要我們介紹了置換檢驗(yàn)贤斜，本次推文主要介紹自助法。
自助法
自助法逛裤，即從初始樣本重復(fù)隨機(jī)替換抽樣瘩绒，生成一個(gè)或一系列待檢驗(yàn)統(tǒng)計(jì)量的經(jīng)驗(yàn)分布，無需假設(shè)一個(gè)特定的理論分布带族，便可生成統(tǒng)計(jì)量的置信區(qū)間锁荔，并能檢驗(yàn)統(tǒng)計(jì)假設(shè)。
簡單來講蝙砌，自助法就是不依賴變量具體分布形式計(jì)算置信區(qū)間的方法阳堕。
自助法步驟：（以下引自R語言實(shí)戰(zhàn)）
（1）從樣本中隨機(jī)選擇10個(gè)觀測，抽樣后再放回择克。有些觀測可能會被選擇多次恬总，有些可能一直都不會被選中；
（2）計(jì)算并記錄樣本均值肚邢；
（3）重復(fù)1和2一千次越驻；
（4）將1000個(gè)樣本均值從小到大排序；
（5）找出樣本均值2.5％和97.5％的分位點(diǎn)道偷，此時(shí)即初始位置和最末位置的第25個(gè)數(shù)缀旁，它們就限定了95%的置信區(qū)間。
針對不服從正態(tài)分布的樣本均值勺鸦，自助法的優(yōu)勢十分明顯并巍，因此常用于潛在分布未知、出現(xiàn)離群值换途、樣本量過小懊渡、沒有可選的參數(shù)方法等情況。
對單個(gè)統(tǒng)計(jì)量使用自助法
rsq<-function(formula,data,indices){
d<-data[indices,]
fit<-lm(formula,data=d)
return(summary(fit)$r.square)
}
library(boot)
set.seed(1234)
results<-boot(data=mtcars,statistic=rsq,R=1000,formula=mpg~wt+disp)
print(results)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mtcars, statistic = rsq, R = 1000, formula = mpg ~
## wt + disp)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.7809306 0.01379126 0.05113904
plot(results)
#計(jì)算95%的置信區(qū)間
boot.ci(results,type=c("perc","bca"))
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = results, type = c("perc", "bca"))
##
## Intervals :
## Level Percentile BCa
## 95% ( 0.6753, 0.8835 ) ( 0.6344, 0.8561 )
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable
我們來分析一下這段代碼：
首先军拟，我們定義了一個(gè)rsq函數(shù)剃执，這個(gè)函數(shù)表示對于formula計(jì)算R方值。此處主要需要修改的代碼是 return(summary(fit)$r.square)
其次懈息，我們進(jìn)行自主法肾档，重復(fù)1000次。此處主要需要修改的代碼是data和formula
最后，計(jì)算95%置信區(qū)間怒见，使用“perc”和“bca”的方法俗慈。
從圖中，我們可以看出自助的R方不呈正態(tài)分布遣耍。
依據(jù)分位數(shù)法計(jì)算置信區(qū)間：Percentile( 0.6753, 0.8835 )
依據(jù)偏差對置信區(qū)間進(jìn)行調(diào)整：BCa ( 0.6344, 0.8561 )
對多個(gè)統(tǒng)計(jì)量使用自助法
bs<-function(formula,data,indices){
d<-data[indices,]
fit<-lm(formula,data=d)
return(coef(fit))
}
library(boot)
set.seed(1234)
results<-boot(data=mtcars,statistic=bs,R=1000,formula=mpg~wt+disp)
print(results)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mtcars, statistic = bs, R = 1000, formula = mpg ~
## wt + disp)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 34.96055404 4.715497e-02 2.546106756
## t2* -3.35082533 -4.908125e-02 1.154800744
## t3* -0.01772474 6.230927e-05 0.008518022
boot(data=mtcars,statistic=bs,R=1000,formula=mpg~wt+disp)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mtcars, statistic = bs, R = 1000, formula = mpg ~
## wt + disp)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 34.96055404 0.2734867271 2.523874212
## t2* -3.35082533 -0.1145055892 1.171195882
## t3* -0.01772474 0.0002178158 0.008599433
plot(results,index=2)
## 計(jì)算車重回歸系數(shù)的置信區(qū)間
boot.ci(results,type="bca",index=2)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = results, type = "bca", index = 2)
##
## Intervals :
## Level BCa
## 95% (-5.477, -0.937 )
## Calculations and Intervals on Original Scale
## 計(jì)算發(fā)動機(jī)排量回歸系數(shù)的置信區(qū)間
boot.ci(results,type="bca",index=3)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = results, type = "bca", index = 3)
##
## Intervals :
## Level BCa
## 95% (-0.0334, -0.0011 )
## Calculations and Intervals on Original Scale
本例計(jì)算的是截距闺阱、車重和發(fā)動機(jī)排量的回歸系數(shù)的置信區(qū)間。
注意事項(xiàng)：
①初始樣本大小10~30即可得到足夠好的結(jié)果舵变。
②重復(fù)1000次一般可滿足需求酣溃。
注意，如果樣本代表性不佳纪隙，或者樣本量過小無法反映總體情況赊豌，使用基于隨機(jī)化和重抽樣的統(tǒng)計(jì)方法也無法將其轉(zhuǎn)化為好數(shù)據(jù)，因此分析數(shù)據(jù)的準(zhǔn)確是結(jié)論準(zhǔn)確的前提瘫拣。
好了亿絮，今天的R語言實(shí)現(xiàn)統(tǒng)計(jì)方法系列推文暫時(shí)告一段落，我們下次再見吧麸拄！ 小伙伴們?nèi)绻惺裁唇y(tǒng)計(jì)上的問題派昧，或者如果想要學(xué)習(xí)什么方面的生物信息內(nèi)容，可以在微信群或者知識星球提問拢切，沒準(zhǔn)哪天的推文就是專門解答你的問題哦蒂萎！
