如何为多响应变量的lmer模型实现非参数Bootstrap？

问题背景 #

需为数据中第18-1543列（示例为第5-10列）的多个响应变量执行lmer模型的非参数Bootstrap，同时要提取预测值绘制95%置信区间，现有代码在编写fit_model函数时持续报错。

核心解决思路 #

1. 动态构建模型公式：通过响应变量列名拼接公式字符串，解决多列响应变量的公式指定问题
2. 标准化模型拟合函数：编写可遍历所有响应变量的拟合函数，同时输出系数和预测值
3. 并行Bootstrap执行：结合parLapply实现并行抽样与模型拟合，提升效率

完整代码实现 #

1. 初始化并行环境 #

library(lme4)
library(parallel)

# 创建并行集群
cl <- makeCluster(3)
# 在集群节点加载lme4包
clusterEvalQ(cl, library(lme4))

2. 定义拟合函数fit_model #

fit_model <- function(boot_data, resp_vars) {
  result_list <- list()
  
  for (var in resp_vars) {
    # 动态拼接模型公式
    formula <- as.formula(paste(var, "~ Timepoint + smoking + Gender + (1 + Timepoint|ID_new)"))
    
    # 拟合lmer模型
    mod <- lmer(formula, data = boot_data, REML = TRUE, na.action = na.omit)
    
    # 提取系数表
    coef_tbl <- as.data.frame(summary(mod)$coefficients)
    
    # 生成预测值（保留随机效应）
    preds <- predict(mod, newdata = boot_data, re.form = NULL)
    
    # 存储当前响应变量的结果
    result_list[[var]] <- list(coefficients = coef_tbl, predictions = preds)
  }
  
  return(result_list)
}

3. 执行并行Bootstrap #

# 指定响应变量列名（示例为第5-10列，实际替换为18:1543）
resp_vars <- colnames(df1)[5:10]

# 向集群节点导出所需对象
clusterExport(cl, varlist = c("df1", "fit_model", "resp_vars"))

# 设置Bootstrap次数
num_boots <- 1000

# 执行并行Bootstrap
bootstrap_results <- parLapply(cl, 1:num_boots, function(boot_idx) {
  # 按行抽样（若需按聚类单元抽样，见下方注意事项）
  boot_data <- df1[sample(nrow(df1), replace = TRUE), ]
  # 拟合所有响应变量的模型
  fit_model(boot_data, resp_vars)
})

# 关闭并行集群
stopCluster(cl)

4. 提取结果并计算95%置信区间 #

# 按响应变量整理所有Bootstrap预测值
pred_list <- lapply(resp_vars, function(var) {
  do.call(rbind, lapply(bootstrap_results, function(boot_res) boot_res[[var]]$predictions))
})
names(pred_list) <- resp_vars

# 计算每个观测值的95%置信区间
ci_list <- lapply(pred_list, function(preds) {
  apply(preds, 2, function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE))
})

# 示例：绘制第一个响应变量的预测值置信区间
library(ggplot2)
var_to_plot <- resp_vars[1]
plot_data <- cbind(df1, lower = ci_list[[var_to_plot]][1, ], upper = ci_list[[var_to_plot]][2, ])

ggplot(plot_data, aes(x = Timepoint, y = .data[[var_to_plot]])) +
  geom_point(size = 1) +
  geom_line(aes(y = predict(lmer(as.formula(paste(var_to_plot, "~ Timepoint + smoking + Gender + (1 + Timepoint|ID_new)")), data = df1))), color = "blue") +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "blue") +
  facet_wrap(~ID_new) +
  labs(title = paste(var_to_plot, "的Bootstrap置信区间"), x = "时间点", y = "响应值")

关键注意事项 #

• 聚类单元抽样优化：混合模型更推荐按聚类单元（如ID_new）抽样，避免破坏组内相关性，调整抽样代码如下：

# 按个体抽样的实现
unique_subjects <- unique(df1$ID_new)
boot_subjects <- sample(unique_subjects, replace = TRUE)
boot_data <- do.call(rbind, lapply(boot_subjects, function(subj) df1[df1$ID_new == subj, ]))

• 预测值基准选择：若需基于原数据的固定效应预测，可将predict函数中的newdata参数设为原数据的设计矩阵子集。

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