如何用R语言boot函数实现双中位数Bootstrapping方法？

基于你现有代码的后续实现 #

你已经完成了两个样本的独立Bootstrap抽样，接下来只需要提取抽样得到的中位数，计算差值分布并构建置信区间即可：

# 生成模拟数据
x <- rnorm(100, mean = 10, sd = 2)
y <- rnorm(100, mean = 12, sd = 3)

# 加载boot包并完成独立Bootstrap抽样
library(boot)
bx <- boot(data = x, statistic = function(x,i) median(x[i]), R = 1000)
by <- boot(data = y, statistic = function(y,i) median(y[i]), R = 1000)

# 1. 计算中位数差值的抽样分布（这里定义为x中位数 - y中位数）
diff_dist <- bx$t - by$t

# 2. 构建95%置信区间（百分位法）
ci <- quantile(diff_dist, c(0.025, 0.975))
cat("95%置信区间：", round(ci[1], 3), "到", round(ci[2], 3), "\n")

# 3. 检验原假设：若置信区间不含0则拒绝原假设
if(ci[1] > 0 || ci[2] < 0) {
  cat("拒绝原假设：两样本中位数存在显著差异\n")
} else {
  cat("不拒绝原假设：没有足够证据表明两样本中位数存在差异\n")
}

更高效的整合式实现 #

推荐直接在Bootstrap统计量中定义中位数差值，一次完成抽样和差值计算，还能使用boot.ci函数生成更准确的置信区间（如BCa法）：

# 生成模拟数据
x <- rnorm(100, mean = 10, sd = 2)
y <- rnorm(100, mean = 12, sd = 3)
combined_data <- list(x = x, y = y)

# 定义Bootstrap统计量：每次迭代同时对两个样本抽样并计算中位数差值
median_diff_stat <- function(data, .) {
  boot_x <- sample(data$x, replace = TRUE)
  boot_y <- sample(data$y, replace = TRUE)
  median(boot_x) - median(boot_y)
}

# 执行Bootstrap抽样
boot_result <- boot(data = combined_data, statistic = median_diff_stat, R = 1000)

# 提取差值的抽样分布
diff_dist <- boot_result$t

# 生成95%置信区间（同时输出百分位法和BCa法）
ci_perc <- boot.ci(boot_result, type = "perc")
ci_bca <- boot.ci(boot_result, type = "bca")

cat("百分位法95%置信区间：", round(ci_perc$percent[4],3), "到", round(ci_perc$percent[5],3), "\n")
cat("BCa法95%置信区间：", round(ci_bca$bca[4],3), "到", round(ci_bca$bca[5],3), "\n")

# 检验原假设（以BCa法结果为例）
if(ci_bca$bca[4] > 0 || ci_bca$bca[5] < 0) {
  cat("拒绝原假设：两样本中位数存在显著差异\n")
} else {
  cat("不拒绝原假设：没有足够证据表明两样本中位数存在差异\n")
}

关键注意事项 #

• 差值方向：明确你定义的差值是x中位数 - y中位数还是y中位数 - x中位数，置信区间的解读要对应这个方向。
• Bootstrap次数：建议设置R=1000以上，样本量较小时可增加到R=5000，提升置信区间的稳定性。
• 置信区间方法：BCa法（偏置校正加速法）比简单百分位法更准确，能修正抽样分布的偏态问题，优先使用。

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