R语言中分组Bootstrap抽样获取回归系数置信区间问题

我有一个包含多地点响应比（RR）的数据框，每个地点被划分到不同区域组。需要为每个组构建回归模型：以RR为因变量，地点作为重复单元，最初用土壤类型做预测因子，后续扩展为多元回归（新增温度作为预测变量），同时要通过Bootstrap重抽样生成各预测变量系数的置信区间，不清楚具体实现方法。

示例数据与初始代码 #

初始示例数据 #

df <- data.frame(
  group=rep(c('region1','region2'), 100),
  subgroup=rep(c('location1','location2',
                 'location2', 'location1'), 25),
  predictor = rep(c('soil1','soil2','soil3','soil4'), 25),
  RR=rnorm(200)
)

含Bootstrap的多元回归代码 #

library(boot)

bootfun <- function(data, i) {
  d <- data[i,]
  fit <- lm(RR ~ soil_type + temperature, data = d)
  coef(fit)
}

set.seed(2022)
set.seed(123)
df <- data.frame(
  group=rep(c('region1','region2'), 100),
  subgroup=rep(c('location1','location2',
                 'location2', 'location1'), 25),
  soil_type = rep(c('soil1','soil2','soil3','soil4'), 25),
  temperature = abs(rnorm(100, 2,1.75)),
  RR=rnorm(200),
stringsAsFactors = TRUE
)

R <- 1000
b_list <- by(df, df$group, \(X) {
  boot(X, bootfun, R, strata = X$subgroup)
})

b_list$region1

关键修正与说明 #

1. 数据一致性修正：原代码中temperature仅生成100个值，但数据框有200行，需调整为rnorm(200, 2, 1.75)避免维度不匹配。
2. 分层Bootstrap合理性：通过strata = X$subgroup按地点分层抽样，确保每个地点的样本被整体抽取，符合“地点为重复单元”的设计，避免破坏同一地点内的数据相关性。

完整实现代码（含置信区间提取） #

library(boot)

# 定义Bootstrap统计量函数：拟合多元回归并返回系数
bootfun <- function(data, i) {
  d <- data[i, ]
  fit <- lm(RR ~ soil_type + temperature, data = d)
  coef(fit)
}

# 生成规范数据
set.seed(123) # 统一随机种子保证可复现
df <- data.frame(
  group = rep(c('region1','region2'), 100),
  subgroup = rep(c('location1','location2', 'location2', 'location1'), 25),
  soil_type = rep(c('soil1','soil2','soil3','soil4'), 25),
  temperature = abs(rnorm(200, 2, 1.75)), # 匹配200行数据
  RR = rnorm(200),
  stringsAsFactors = TRUE
)

# 按区域分组执行Bootstrap，分层单位为地点
R <- 1000
b_list <- by(df, df$group, function(X) {
  boot(data = X, statistic = bootfun, R = R, strata = X$subgroup)
})

# 提取系数与95%置信区间（百分位数法）
get_coef_ci <- function(boot_result) {
  # 遍历每个系数的Bootstrap样本
  ci_results <- lapply(seq_along(boot_result$t0), function(idx) {
    coef_name <- names(boot_result$t0)[idx]
    # 计算百分位数置信区间
    ci <- quantile(boot_result$t[, idx], c(0.025, 0.975))
    data.frame(
      变量名 = coef_name,
      点估计值 = boot_result$t0[idx],
      95%置信区间下限 = ci[1],
      95%置信区间上限 = ci[2],
      row.names = NULL
    )
  })
  do.call(rbind, ci_results)
}

# 获取两个区域的结果
region1_ci <- get_coef_ci(b_list$region1)
region2_ci <- get_coef_ci(b_list$region2)

# 输出结果
cat("=== region1 系数置信区间 ===\n")
print(region1_ci)
cat("\n=== region2 系数置信区间 ===\n")
print(region2_ci)

可选：更准确的BCA置信区间 #

如果需要考虑偏差和分布偏度，可改用BCA（偏差校正加速）置信区间，修改提取函数如下：

get_coef_ci_bca <- function(boot_result) {
  ci_results <- lapply(seq_along(boot_result$t0), function(idx) {
    coef_name <- names(boot_result$t0)[idx]
    # 计算BCA置信区间
    bca_ci <- boot.ci(boot_result, index = idx, type = "bca")$bca[, 4:5]
    data.frame(
      变量名 = coef_name,
      点估计值 = boot_result$t0[idx],
      95%置信区间下限 = bca_ci[1],
      95%置信区间上限 = bca_ci[2],
      row.names = NULL
    )
  })
  do.call(rbind, ci_results)
}

# 使用BCA方法获取置信区间
region1_ci_bca <- get_coef_ci_bca(b_list$region1)
print(region1_ci_bca)

内容的提问来源于stack exchange，提问作者Joshua Smith