如何使用 R 中的“引导函数”计算置信区间

Question

我正在尝试计算 R 中的置信区间。由于一些特殊原因，我必须使用“bootstrap”包中的函数来计算。（这意味着我不能使用“boot”包中的函数。 )

这是我的代码。

而我正在做的是尝试计算 Pearson 相关系数，然后应用 Bootstrap 方法（B = 100）来获得相关系数的估计值。但我不知道如何构造 95% 的置信区间。

library(bootstrap) 
data('law')

set.seed(1)
theta <- function(ind) {
  cor(law[ind, 1], law[ind, 2], method = "pearson")
  }
law.boot <- bootstrap(1:15, 100, theta) 
# sd(law$thetastar)
percent.95 <- function(x) {
  quantile(x,  .95)
  }
law.percent.95 <- bootstrap(1:15, 100, theta, func=percent.95)

对不起，如果我没有让自己清楚或标记错误的标签。 很抱歉没有生成数据集（现在已提供），并感谢 Roland 教授指出。非常感谢！

Answer 1

通常，在引导之后，我们使用 2.5% 和 97.5% 的百分位数作为 95% 的置信区间（因为我们从每一边减去 α/2=.025）。另请参阅 @thothal 的 answer 和答案下的 cmets。

R <- 1e5 - 1  ## number of bootstrap replications
est <- with(law, cor(lsat, gpa))  ## naïve correlation

theta <- function(ind) cor(law[ind, 1], law[ind, 2], method="pearson")
set.seed(1)
B1 <- bootstrap::bootstrap(seq(nrow(law)), R, theta) 
(ci1 <- c(estimate=est, quantile(B1$thetastar, c(.025, .975))))
#  estimate      2.5%     97.5% 
# 0.7763745 0.4594845 0.9620884 

这里是一种从头开始的替代方法：

theta2 <- function(x) with(x, cor(lsat, gpa))
set.seed(1)
B2 <- replicate(R, theta2(law[sample(nrow(law), nrow(law), replace=TRUE), ]))
(ci2 <- c(estimate=est, quantile(B2, c(.025, .975))))
#  estimate      2.5%     97.5% 
# 0.7763745 0.4607644 0.9617970 

最后是使用具有boot.ci 函数的boot 包的方法：

theta3 <- function(data, k) cor(data[k, ])[1,2]
set.seed(1)
B3 <- boot::boot(law, theta3, R=R)
(ci3 <- c(est, boot::boot.ci(B3, type='perc')$percent[4:5]))
# [1] 0.7763745 0.4593727 0.9620923

Answer 2

有多种计算引导估计器 CI 的方法（参见to this Wikipedia article for instance。

最简单的方法是从自举系数（维基百科文章中的 Percentile Bootstrap）中找出 2.5% 和 97.5% 分位数：

quantile(law.boot$thetastar, c(0.025, 0.975))
#      2.5%     97.5% 
# 0.4528745 0.9454483 

Basic Bootstrap 将被计算为

2 * mean(law.boot$thetastar) - quantile(law.boot$thetastar, c(0.975, 0.025))
#     97.5%      2.5% 
# 0.5567887 1.0493625