使用R实现非参数bootstrap计算线性模型预测值95%枢轴置信区间

问题定位与修正方法 #

你的代码存在3个会导致结果错误的核心问题，逐一对应如下：

• 未指定预测目标x值，预测结果完全偏离目标场景
  代码中predict(fit)默认返回重抽样后模型对原数据集所有样本的预测向量，你将这个长向量和单个残差相加后赋值给pred[ii]（标量），R会自动循环运算并只取第一个元素的结果，根本不是你需要的x=-6.10529场景的预测值。
• 错误引入残差导致区间人为放大
  你要构建的是预测值（条件均值E(Y|X=x0)）的置信区间，不需要额外加残差。残差是单个观测的随机噪声，只有在构建单个新观测的预测区间时才需要加入，这就是你觉得区间过宽的最主要原因。
• 混淆置信区间计算逻辑
  你最后直接取bootstrap结果的2.5%和97.5%分位数，得到的是百分位置信区间，不是要求的枢轴（基本法）置信区间。枢轴法的计算逻辑需要基于原始预测值反推，不是直接取分位数。

正确实现代码 #

首先明确前提：模型自变量是log(Metabolic.by.mass)，目标预测值对应的x0=-6.10529是该变量的取值，采用和你原本一致的案例重抽样非参数bootstrap逻辑。

# 定义目标预测点
x0 <- -6.10529

# 计算原始模型的目标预测值，可验证是否为17.34775
original_fit <- lm(lifespan ~ log(Metabolic.by.mass), data = anage)
y_hat_original <- predict(original_fit, newdata = data.frame(Metabolic.by.mass = exp(x0)))

B <- 1000
boot_pred <- numeric(B)
n <- nrow(anage)

for (ii in 1:B) {
  # 重抽样本索引
  boot_idx <- sample(seq_len(n), size = n, replace = TRUE)
  boot_data <- anage[boot_idx, ]
  # 拟合重抽样模型
  boot_fit <- lm(lifespan ~ log(Metabolic.by.mass), data = boot_data)
  # 对固定x0预测，不添加残差
  boot_pred[ii] <- predict(boot_fit, newdata = data.frame(Metabolic.by.mass = exp(x0)))
}

# 枢轴法置信区间计算：公式为 [2*原始预测值 - 97.5%分位数, 2*原始预测值 - 2.5%分位数]
ci_lower <- 2 * y_hat_original - quantile(boot_pred, 0.975)
ci_upper <- 2 * y_hat_original - quantile(boot_pred, 0.025)

# 输出结果
cat("95%枢轴置信区间为：", ci_lower, " ~ ", ci_upper, "\n")

补充说明 #

• 如果你确实需要构建单个新观测的预测区间，再在predict(boot_fit, ...)后加上+ sample(resid(boot_fit), size = 1)即可，注意区分置信区间（针对条件均值）和预测区间（针对单个观测）的差异。
• 为了避免newdata赋值出错，可提前对数据集做预处理，把log后的变量单独存为一列，降低出错概率：

  anage$log_met <- log(anage$Metabolic.by.mass)
  original_fit <- lm(lifespan ~ log_met, data = anage)
  y_hat_original <- predict(original_fit, newdata = data.frame(log_met = x0))
  # 后续bootstrap的newdata直接传log_met=x0即可，无需计算指数
  

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