当样本不符合理论分布假设时,求样本统计量的置信区间就成为一个难题。而自助法(Bootstrap)的思路是对原始样本重复抽样产生多个新样本,针对每个样本求取统计量,然后得到它的经验分布,再通过求经验分布的分位数来得到统计量的置信区间,这种方法不需要对统计量有任何理论分布的假设。一般认为,只要样本具有代表性,采用自助法需要的原始样本只要20-30个,重复抽样1000次就能达到满意的结果。
在R中进行自助法是利用boot扩展包,其流程如下:
-
编写一个求取统计量的自定义函数
-
将上面的函数放入boot()函数中进行运算,得到自助法的结果
-
用boot.ci()函数求取置信区间
让我们用mtcars数据集来作为例子,我们可以将wt和disp作为自变量,mpg 作为因变量,进行回归后能得到一系列回归统计量。其中我们感兴趣的是判定系数R-square,希望用自助法求它的95%置信区间。
首先定义求R-square的函数,注意其中的indices是必不可少的参数,另外一个参数代表样本数据
————————
rsq=function(data,indices){
d=data[indices,]
fit=lm(formula=mpg~wt+disp,data=d)
return(summary(fit)$r.square)
}
————————
载入boot扩展包,将随机种子设为1234,以方便得到相同的结果,再利用boot函数得到结果results,其中R表示重复抽样得到1000个样本
————————
library(boot)
set.seed(1234)
results=boot(data=mtcars,statistic=rsq,R=1000)
print(results)
plot(results)
————————
results这个数据结构中包括了原始样本的统计量(results$t0)和再抽样样本的统计量(results$t0),上图左侧的直方图表示了再抽样样本的统计量的经验分布,其中的虚线表示了原始样本的统计量,从中可以观察到偏差。右侧QQ图有助于判断经验分布是否正态。下面我们用boot.ci函数从结果中提取置信区间。
————————
boot.ci(results,conf=0.95,type=c('perc','bca'))
————————
其中conf表示置信水平,type表示了用何种算法来求区间,perc即使用百分位方法,bca表示adjusted bootstrap percentile,即对偏差进行了调整。结果如下:
BOOTSTRAPCONFIDENCEINTERVALCALCULATIONS Basedon1000bootstrapreplicates CALL: boot.ci(boot.out=results,conf=0.95,type=c("perc","bca")) Intervals: LevelPercentileBCa 95%(0.6838,0.8833)(0.6344,0.8549)
来源:http://www.r-bloggers.com/lang/chinese/540
1、回复“数据分析师”查看数据分析师系列文章
2、回复“案例”查看大数据案例系列文章
3、回复“征信”查看相关征信的系列文章
4、回复“可视化”查看可视化专题系列文章
5、回复“SPPS”查看SPSS系列文章
6、回复“答案”查看hadoop面试题题目及答案
7、回复“爱情”查看大数据与爱情的故事
8、回复“笑话”查看大数据系列笑话
9、回复“大数据1、大数据2、大数据3、大数据4”查看大数据历史机遇连载
原文始发于微信公众号(PPV课数据科学社区):【数据分析】R语言中用自助法求统计量置信区间
原创文章,作者:ppvke,如若转载,请注明出处:http://www.ppvke.com/archives/18031