如何基于长格式数据集通过迭代方法快速拟合混合模型——针对多MARKER分组的lmer模型迭代拟合与诊断需求

没问题，咱们一步步来解决你这三个关于混合模型迭代拟合的问题，结合你提供的长格式数据集来展开：

问题1：长格式数据下用迭代函数拟合混合模型 #

完全可以！你的长格式数据刚好适合按MARKER分组后迭代拟合模型。我们可以用tidyverse里的purrr::map()函数（比for循环更简洁），配合dplyr::group_split()按MARKER拆分数据集，然后对每个子集调用lmer()。

具体代码示例：

# 加载所需包
library(lme4)
library(tidyverse)

# 假设你的数据集名为df（替换成你实际的数据集名称）
df <- structure(list(
 SAMPLE = c("A1", "A1", "A1", "A2", "A2", "A2", "A3", "A3", "A3", "A4"),
 GROUP = c("AA", "AA", "AA", "AA", "AA", "AA", "AA", "AA", "AA", "AA"),
 SESSION = c("L", "L", "L", "L", "L", "L", "L", "L", "L", "L"),
 CATEGORY = c("LUM-STD", "LUM-ALT", "CON-ALT", "LUM-STD", "LUM-ALT", "CON-ALT", "LUM-STD", "LUM-ALT", "CON-ALT", "LUM-STD"),
 MARKER = factor(rep("Z1Fx", 10), levels = "Z1Fx"),
 READING = c(-11.63, -11.14, -3.996, -0.314, 0.239, 5.037, -0.214, -2.96, -1.97, -2.83)
), row.names = c(NA, -10L), class = c("tbl_df", "tbl", "data.frame"))

# 按MARKER分组，迭代拟合模型
model_list <- df %>%
  group_split(MARKER) %>%
  map(~ lmer(READING ~ CATEGORY + (1 | SAMPLE), data = .x))

# 给模型列表命名（方便后续索引）
names(model_list) <- unique(df$MARKER)

# 查看某个模型的结果，比如Z1Fx的模型
summary(model_list[["Z1Fx"]])

这里的核心逻辑是：每个MARKER对应一个独立的数据集子集，每个子集都拟合以CATEGORY为固定效应、SAMPLE为随机截距的混合模型。如果需要更复杂的随机效应结构（比如随机斜率），只需要修改公式里的随机部分即可，比如(CATEGORY | SAMPLE)。

如果你更习惯用for循环，也可以这么写：

model_list <- list()
unique_markers <- unique(df$MARKER)

for (marker in unique_markers) {
  subset_data <- filter(df, MARKER == marker)
  model_list[[marker]] <- lmer(READING ~ CATEGORY + (1 | SAMPLE), data = subset_data)
}

问题2：宽格式数据集的迭代分析代码 #

如果你的数据转成宽格式（每个MARKER对应一列READING值，比如Z1Fx_READING、Z2Fx_READING等），我们可以通过迭代列名来拟合模型。假设宽格式数据集名为wide_df，其中列名包含每个MARKER的读数：

# 提取所有MARKER对应的读数列名（假设列名格式为"MARKER_READING"）
reading_cols <- grep("_READING", colnames(wide_df), value = TRUE)
# 提取MARKER名称
marker_names <- str_remove(reading_cols, "_READING")

# 迭代拟合模型
wide_model_list <- map2(reading_cols, marker_names, function(col, marker) {
  # 构造公式：读数列 ~ CATEGORY + (1 | SAMPLE)
  formula_str <- paste(col, "~ CATEGORY + (1 | SAMPLE)")
  lmer(as.formula(formula_str), data = wide_df)
})

names(wide_model_list) <- marker_names

这里的关键是动态构造公式，把每个MARKER对应的读数列作为因变量，再统一用CATEGORY和SAMPLE作为固定/随机效应。

问题3：迭代提取模型诊断图 #

当然可以！lme4模型支持基础的诊断图（比如残差拟合图、QQ图），我们可以用map()结合plot()函数，或者用ggplot2来生成更美观的诊断图。

方法1：用基础绘图函数生成并保存诊断图 #

# 遍历模型列表，生成并保存每个模型的诊断图
walk2(model_list, names(model_list), function(model, marker) {
  # 打开一个png文件
  png(paste0("model_diagnostic_", marker, ".png"), width = 800, height = 600)
  # 生成4个诊断图：残差拟合、残差QQ、尺度-位置、残差-leverage
  par(mfrow = c(2, 2))
  plot(model)
  # 关闭绘图设备
  dev.off()
})

方法2：用ggplot2生成诊断图（更灵活） #

如果你想用ggplot2，可以先提取模型的残差和拟合值，再绘图：

# 提前安装需要的包
install.packages(c("broom.mixed", "patchwork"))

# 加载包
library(broom.mixed)
library(patchwork)
library(ggplot2)

# 定义一个生成诊断图的函数
make_diagnostic_plots <- function(model, marker) {
  # 提取模型数据和残差
  model_data <- augment(model)
  
  # 残差拟合图
  p1 <- ggplot(model_data, aes(x = .fitted, y = .resid)) +
    geom_point() +
    geom_hline(yintercept = 0, color = "red") +
    labs(title = paste(marker, "残差vs拟合值"), x = "拟合值", y = "残差")
  
  # QQ图
  p2 <- ggplot(model_data, aes(sample = .resid)) +
    stat_qq() +
    stat_qq_line() +
    labs(title = paste(marker, "残差QQ图"))
  
  # 组合两个图
  wrap_plots(p1, p2, ncol = 2)
}

# 迭代生成所有模型的诊断图
diagnostic_plots <- map2(model_list, names(model_list), make_diagnostic_plots)

# 查看某个MARKER的诊断图
diagnostic_plots[["Z1Fx"]]

# 保存所有图
walk2(diagnostic_plots, names(diagnostic_plots), function(plot, marker) {
  ggsave(paste0("ggplot_diagnostic_", marker, ".png"), plot, width = 10, height = 6)
})

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