如何用R的boot包实现野抽样并获取BCa置信区间？

核心问题澄清 #

1. 野抽样的本质：野抽样属于参数化Bootstrap的变体，它依赖原模型拟合得到的残差，通过对残差施加随机权重（如Mammen权重）生成新的响应变量，而非直接对观测数据重抽样。因此必须使用sim="parametric"实现，无法直接改为非参数抽样（非参数抽样是对原始(X,Y)观测对有放回重抽样，逻辑完全不同）。

2. 参数化vs非参数抽样的区别

   • 非参数Bootstrap：完全基于观测数据的经验分布，对原始观测对重抽样，无需依赖模型假设，适用于对数据分布无先验认知的场景。
   • 参数化Bootstrap：先拟合模型得到参数/残差，再基于模型设定生成新数据（如野抽样的残差权重变换、正态分布残差模拟），依赖模型的正确性，适用于模型假设合理的场景。

解决BCa置信区间报错的方案 #

报错参数化抽样无法获取经验影响值的原因是：BCa区间需要每个观测对估计量的经验影响值（即去掉单个观测后的估计量与全样本估计量的差值），而参数化Bootstrap生成的boot对象默认不会计算并存储这些值。解决方法是手动计算影响值，传递给boot.ci()。

步骤1：基础野抽样代码 #

先构造数据并实现野抽样：

library(boot)
set.seed(123)

# 模拟OLS数据
n <- 50
x <- rnorm(n)
y <- 1 + 2*x + rnorm(n, 0, 0.5)
fit <- lm(y ~ x)

# 野抽样统计量函数（parametric模式专用）
wild_stat <- function(data, indices, model) {
  # 生成Mammen权重（野抽样常用权重）
  mammen_weights <- sample(c((1+sqrt(5))/2, (1-sqrt(5))/2), 
                           nrow(data), replace = TRUE, 
                           prob = c(0.5, 0.5))
  # 生成新响应变量
  y_star <- fitted(model) + mammen_weights * residuals(model)
  # 返回新模型的系数
  coef(lm(y_star ~ data$x))
}

# 执行野抽样
boot_wild <- boot(data = data.frame(x, y),
                  statistic = wild_stat,
                  R = 1000,
                  sim = "parametric",
                  model = fit)

步骤2：手动计算经验影响值 #

通过留一法（Leave-One-Out, LOO）计算每个观测对系数的影响值：

# 计算留一法的系数估计
loo_coefs <- t(sapply(1:n, function(i) {
  coef(lm(y ~ x, data = data.frame(x, y)[-i, ]))
}))

# 计算影响值L_i = 留一法估计值 - 全样本估计值
theta_full <- coef(fit)
influence_vals <- t(t(loo_coefs) - theta_full)

步骤3：生成BCa置信区间 #

将手动计算的影响值通过L参数传入boot.ci()，或直接赋值给boot对象的influence元素：

# 方法1：直接传入L参数
bca_ci <- boot.ci(boot_wild, type = "bca", L = influence_vals)
print(bca_ci)

# 方法2：修改boot对象的influence元素，后续调用无需重复传参
boot_wild$influence <- influence_vals
bca_ci <- boot.ci(boot_wild, type = "bca")

关于修改empinf()的说明 #

Stack Overflow提到的修改empinf()参数本质上是绕开boot包默认的影响值计算逻辑，直接传入手动计算的结果。上述两种方法已经等价于这个思路——因为boot.ci()在检测到L参数或boot$influence存在时，会跳过调用empinf()，直接使用传入的影响值，因此是完全合法的解决方案。

内容的提问来源于stack exchange，提问作者hanna-without-h