R:从引导结果向量计算 BCa

Question

我正在寻找一种方法，使用引导结果向量(这是人口增长率的引导估计 - lambda)来计算 R 中的偏差校正加速置信区间。然而，我发现的包要么是使用特定的对象类型(如“boot”包中)，要么不计算 BCa 类型置信区间。我使用 for 循环引导结果然后将结果存储在向量中的原因是，对于每个引导重新采样，我首先得到一个 80 x 33 的结果矩阵，它定义了每年采样中每个总体的参数，而这些参数又为每个群体定义 lambda。据我所知，这在引导包中会很麻烦，并且很容易作为 for 循环进行编程。实际的功能集相当复杂，无法包含在这里。

我确实尝试使用这个问题作为伪造“启动”对象的指南，但它不起作用:How can I use pre bootstrapped data to obtain a BCa confidence interval? 。

假设我有观察到的 lambda 估计

lambda = 1.18

并且我们模拟引导估计的向量

library(fGarch)
lambdaBS = rsnorm(999,mean=lambda-0.04,sd=0.11,xi=2.5)
plot(density(lambdaBS))

这是右偏和有偏见的。

我希望，使用这些信息，当前存在一个函数可以计算 BCa 置信区间，否则很容易编写一个函数来执行此操作。到目前为止，我还没有发现这种情况。

Answer 1

与 R 的情况一样，某些实用程序在包之间广泛分布，有一个简单的解决方案，但我花了几个小时的搜索时间才找到，所以我将为任何可能的人回答我自己的问题寻找类似的东西。

使用问题中的示例数据，“coxed”R 包中的 bca 函数为引导结果向量提供偏差校正和加速置信区间。我们可以将它们与其他置信区间进行比较。

library(fGarch)
library(coxed)

set.seed(15438)

#simulate bootstrap statistics
lambdaBS = rsnorm(9999,mean=lambda-0.04,sd=0.11,xi=2.5)

#bias-corrected and accelerated
bca(lambdaBS)

1.002437 1.452525

#confidence intervals using standard error (inappropriate)
c(lambda-(sd(lambdaBS)*2),lambda+(sd(lambdaBS)*2))

0.9599789 1.4000211

#percentile confidence intervals
quantile(lambdaBS, c(0.025,0.975))

    2.5%     97.5% 
0.9895892 1.4016528 

这似乎运作良好。我不确定它如何在不需要对相关统计数据进行初始估计的情况下纠正偏差，但尚未阅读该方法所基于的论文。

另一个模拟展示了它与使用 boot 和 boot.ci 的结果的比较。

library(boot)

#generate data
set.seed(12345)
dat = rsnorm(500,mean=1.6,sd=0.5,xi=3.0)

#bootstrap the median
meanfun = function(x,id){ mean(x[id])}
test = boot(data=dat,R=999,statistic=meanfun)

#BCa using boot.ci
boot.ci(test,type="bca")

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = test, type = "bca")

Intervals : 
Level       BCa          
95%   ( 1.537,  1.626 )  
Calculations and Intervals on Original Scale


#BCa using bca function from coxed package
bca(test$t)

 1.536888 1.625524

在这种情况下，两个函数给出相同的结果。