在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 )