混合效应模型拟合数据的Bootstrap置信区间偏移问题

问题背景 #

我在Stack Exchange统计学板块提问了如何在R中使用bootMer为混合效应模型的预测结果生成有效的Bootstrap置信区间与预测区间，涵盖两种场景：1）模型拟合时已见过的数据；2）未见过的新数据。我采纳了回答中的方法，基于sleepstudy数据集构建了一个简单混合效应模型，其中Days为固定效应，Subject为随机效应。可复现的R代码如下：

# Load necessary libraries
library(lme4)
library(ggplot2)

# Load the sleepstudy dataset
data("sleepstudy")

# Fit the mixed effect model
mem <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)

# Define the subject for the plots
subject <- 337

# Seen subject
seen_data <- sleepstudy[sleepstudy$Subject == subject, ]

# New subject
new_data <- data.frame(Days = seq(min(sleepstudy$Days), max(sleepstudy$Days), length.out = 10),
                       Subject = factor(rep("999", 10)))

# Define prediction functions for bootstrapping
predfn <- function(.) { predict(., newdata = seen_data, re.form = NULL) }
cifn <- function(.) { predict(., newdata = seen_data, re.form = NA) }
sfun <- function(.) { simulate(., newdata = new_data, re.form = NA, allow.new.levels = TRUE)[[1]] }
cifn_new <- function(.) { predict(., newdata = new_data, re.form = NA) }

# Perform bootstrapping for seen subject
boot_pred <- bootMer(mem, FUN = predfn, nsim = 500, re.form = NULL, type = "parametric", seed = 2804)
boot_ci <- bootMer(mem, FUN = cifn, nsim = 500, re.form = NA, type = "parametric", seed = 2804)

# Perform bootstrapping for new subject
boot_pred_new <- bootMer(mem, FUN = sfun, nsim = 500, re.form = NA, type = "parametric", seed = 2804)
boot_ci_new <- bootMer(mem, FUN = cifn_new, nsim = 500, re.form = NA, type = "parametric", seed = 2804)

# Helper function to calculate confidence and prediction intervals
calc_intervals <- function(boot_obj, level = 0.95) {
  alpha <- 1 - level
  lower <- apply(boot_obj$t, 2, quantile, probs = alpha / 2)
  upper <- apply(boot_obj$t, 2, quantile, probs = 1 - alpha / 2)
  list(lower = lower, upper = upper)
}

# Calculate intervals
ci_seen <- calc_intervals(boot_ci)
pi_seen <- calc_intervals(boot_pred)
ci_new <- calc_intervals(boot_ci_new)
pi_new <- calc_intervals(boot_pred_new)

# Predictions
predictions_seen <- predict(mem, newdata = seen_data, re.form = NULL)
predictions_new <- predict(mem, newdata = new_data, re.form = NA, allow.new.levels = TRUE)

# Plotting
plot_intervals <- function(days, predictions, ci, pi, title) {
  data_plot <- data.frame(Days = days, 
                          Prediction = predictions, 
                          CI_Lower = ci$lower, 
                          CI_Upper = ci$upper, 
                          PI_Lower = pi$lower, 
                          PI_Upper = pi$upper)
  
  ggplot(data_plot, aes(x = Days)) +
    geom_line(aes(y = Prediction, color = "Prediction")) +
    geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper, fill = "Confidence Interval"), alpha = 0.2) +
    geom_ribbon(aes(ymin = PI_Lower, ymax = PI_Upper, fill = "Prediction Interval"), alpha = 0.1) +
    scale_color_manual(name = "Legend", values = c("Prediction" = "blue")) +
    scale_fill_manual(name = "Legend", values = c("Confidence Interval" = "green", "Prediction Interval" = "red")) +
    labs(title = title, y = "Reaction Time", x = "Days") +
    theme_minimal() +
    theme(legend.position = "bottom")
}

# Create plots
plot_seen_ci <- plot_intervals(seen_data$Days, predictions_seen, ci_seen, pi_seen, 
                               "Seen Subject - Confidence and Prediction Intervals")
plot_new_ci <- plot_intervals(new_data$Days, predictions_new, ci_new, pi_new, 
                              "New Subject - Confidence and Prediction Intervals")

# Print plots
print(plot_seen_ci)
print(plot_new_ci)

新主体的可视化结果（）表现正常，但针对模型拟合时已见过的主体，其置信区间相对于模型预测值和预测区间出现偏移（）。需注意的是，cifn与cifn_new完全相同，boot_ci与boot_ci_new也完全一致（遵循此前回答的建议）。核心问题是：针对已见过的主体，使用re.form=NA得到的置信区间并未正确以使用re.form=NULL（包含模型拟合时学习到的随机效应）得到的预测值为中心。

正式问题 #

1. 对混合效应模型拟合过的数据，计算其预测结果的Bootstrap置信区间是否具备统计学意义？
2. 若具备意义，如何让置信区间正确围绕模型预测值（与预测区间一致）？
3. 当前脚本中的函数参数选择是否存在错误？

解答 #

问题1：统计学意义 #

有意义，但需要明确置信区间的定义：

• 针对已见过的主体计算置信区间时，存在两种不同的统计目标：
  1. 该主体的条件均值置信区间：反映对该主体真实平均反应时间的估计不确定性，包含固定效应估计误差+该主体随机效应的估计误差，对应re.form=NULL的预测逻辑。
  2. 总体平均置信区间：反映整个群体在对应Days下的平均反应时间的估计不确定性，与单个主体的随机效应无关，对应re.form=NA的预测逻辑。

你当前代码中用re.form=NA计算的是总体平均的置信区间，而非单个已见主体条件均值的置信区间，所以它不会以单个主体的条件预测值为中心，这是逻辑上的正常结果，而非错误。

问题2：让CI围绕已见主体的预测值 #

要让置信区间围绕re.form=NULL得到的条件预测值，需要重新定义Bootstrap函数，计算条件均值的置信区间：

• 修改针对已见主体的Bootstrap函数，确保每次抽样后预测时包含当前样本拟合出的该主体随机效应：

cifn_seen <- function(.) {
  predict(., newdata = seen_data, re.form = NULL)
}

• 重新运行Bootstrap：

boot_ci_seen <- bootMer(mem, FUN = cifn_seen, nsim = 500, type = "parametric", seed = 2804)

• 用这个新的boot_ci_seen计算的区间，会以predict(mem, newdata=seen_data, re.form=NULL)的结果为中心，因为它捕获了条件均值的全部估计不确定性。

另外，若需计算已见主体的预测区间（包含残差随机波动），应替换当前的predfn为simulate逻辑（和新主体的处理一致）：

sfun_seen <- function(.) {
  simulate(., newdata = seen_data, re.form = NULL)[[1]]
}
boot_pred_seen <- bootMer(mem, FUN = sfun_seen, nsim = 500, type = "parametric", seed = 2804)

问题3：当前脚本的参数错误 #

存在核心逻辑错误：

• 针对已见主体的置信区间，错误使用re.form=NA，导致计算的是总体平均的置信区间，而非该主体条件均值的置信区间，这是区间偏移的直接原因。
• 针对已见主体的预测区间，当前predfn用predict而非simulate，predict默认仅输出条件均值的点估计，不包含残差变异，得到的区间实际是条件均值的置信区间，而非真正的预测区间（预测区间需包含残差的随机波动）。

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