|
對于正態(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)
 img#計(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)
 img## 計(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)哪天的推文就是專門解答你的問題哦！
|
