R中自定义Bootstrap与Wald置信区间函数的实现及报错修复

报错根因拆解 #

• 命名不一致：自定义函数定义为waldconfidenceinterval，但循环内调用时写为waldCI，R无法匹配对应函数对象。
• 输入维度错配：循环内取bootdata[i]作为函数输入，仅为单个向量元素，和函数设计的“输入完整数据集按列计算统计量”的逻辑完全不匹配。
• 赋值长度不匹配：预定义的bootresults矩阵为200行99列，但单轮Wald计算会返回99组上下限共198个数值，直接赋值时长度差为2倍，触发“替换项长度不是原长度整数倍”的报错。
• 重采样逻辑偏差：原代码sample(B, replace=T)仅从1-200的行号中有放回抽样，未覆盖原始数据集全部1500行样本，不符合bootstrap有放回重采样的核心逻辑。
• 计算逻辑偏差：原Wald公式内的分母n*200没有对应实际样本量的计算逻辑，SE计算的基础统计量维度和样本规模不匹配。

目标功能对齐 #

本次实现需要覆盖的核心能力：

• 输出结构化的bootstrap结果，包含4类字段：真实总体p值、Wald置信区间下限、Wald置信区间上限、真实p值是否被区间覆盖
• 自动计算每个p值场景下的经验区间覆盖率
• 生成Wald区间覆盖率随真实p值变化的折线图

可直接运行的修正代码 #

set.seed(42)
# 1. 生成模拟基础数据：每个列对应一个真实p值，共1500个观测，每个观测为10次二项试验的成功数
true_ps <- seq(0.01, 0.99, 0.01)
samples10 <- lapply(true_ps, function(p) rbinom(n = 1500, size = 10, prob = p))
sample10 <- data.frame(samples10)
colnames(sample10) <- true_ps

# 2. 定义修正后的Wald置信区间函数
wald_ci <- function(boot_data, alpha = 0.05) {
  # 计算列合计成功数
  total_success <- colSums(boot_data)
  # 每个bootstrap样本的总试验次数：重采样行数 * 单观测的试验次数(size=10)
  total_trials <- nrow(boot_data) * 10
  # 估计成功概率p_hat
  p_hat <- total_success / total_trials
  # 计算Wald标准误
  se <- sqrt(p_hat * (1 - p_hat) / total_trials)
  # 计算95%CI上下限
  z_val <- qnorm(c(alpha/2, 1 - alpha/2))
  lower <- p_hat + z_val[1] * se
  upper <- p_hat + z_val[2] * se
  # 返回结构化结果
  return(data.frame(lower = lower, upper = upper))
}

# 3. 初始化存储结构
B <- 200 # bootstrap重复次数
n_col <- length(true_ps)
# 分别存储每轮bootstrap的上下限
ci_lower <- matrix(NA, nrow = B, ncol = n_col)
ci_upper <- matrix(NA, nrow = B, ncol = n_col)

# 4. 执行bootstrap循环
set.seed(42)
for(i in seq_len(B)) {
  # 有放回抽取样本索引，得到本轮bootstrap重采样数据集
  boot_idx <- sample(nrow(sample10), size = 200, replace = TRUE)
  boot_data <- sample10[boot_idx, ]
  # 计算本轮的Wald区间
  res <- wald_ci(boot_data)
  # 存储上下限
  ci_lower[i, ] <- res$lower
  ci_upper[i, ] <- res$upper
}

# 5. 计算覆盖率与覆盖标识
# 逐列判断每个真实p值是否被对应轮次的区间覆盖
cover_matrix <- mapply(
  function(true_p, lower_vec, upper_vec) {
    as.integer(lower_vec < true_p & true_p < upper_vec)
  },
  true_ps,
  as.data.frame(ci_lower),
  as.data.frame(ci_upper)
)
# 每个p值对应的经验覆盖率
coverage <- colMeans(cover_matrix)

# 6. 生成结构化结果表
boot_results <- data.frame(
  true_p = rep(true_ps, each = B),
  ci_lower = as.vector(ci_lower),
  ci_upper = as.vector(ci_upper),
  is_covered = as.vector(cover_matrix)
)

# 7. 绘制Wald区间覆盖图
plot(
  x = true_ps,
  y = coverage,
  type = "l",
  lwd = 2,
  xlab = "真实成功概率p",
  ylab = "经验覆盖率",
  ylim = c(0, 1),
  main = "Wald置信区间覆盖情况"
)
# 添加95%名义置信水平参考线
abline(h = 0.95, col = "red", lty = 2)

修正说明 #

• 全程统一函数命名，消除函数找不到的报错
• 调整输入逻辑，每轮将完整重采样数据集传入Wald函数，自动完成99个p值场景的统计量计算，彻底解决长度不匹配的赋值报错
• 修正重采样逻辑，从全部1500行原始数据中有放回抽取样本，符合bootstrap核心设计原则
• 修正Wald公式的样本量计算逻辑，严格匹配二项试验的总试验次数，避免SE计算偏差
• 结构化存储覆盖标识，可直接导出包含4类目标字段的结果表
• 绘图环节添加名义覆盖率参考线，可直观对比不同p值下的区间覆盖表现

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