星期日, 十二月 04, 2011

在R语言中用自助法求统计量置信区间

当样本不符合理论分布假设时,求样本统计量的置信区间就成为一个难题。而自助法(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,即对偏差进行了调整。结果如下:

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 )  

没有评论:

发表评论