在R中进行自助法是利用boot扩展包,其流程如下:
- 编写一个求取统计量的自定义函数
- 将上面的函数放入boot()函数中进行运算,得到自助法的结果
- 用boot.ci()函数求取置信区间
首先定义求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,即对偏差进行了调整。结果如下:
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 1000 bootstrap replicates CALL : boot.ci(boot.out = results, conf = 0.95, type = c("perc", "bca")) Intervals : Level Percentile BCa 95% ( 0.6838, 0.8833 ) ( 0.6344, 0.8549 )
没有评论:
发表评论