如何在R中绘制混合效应模型的预测结果（含95%置信区间）

步骤1：加载必要的R包 #

我们需要用到lme4拟合混合效应模型，ggplot2绘图，dplyr处理数据：

library(lme4)
library(ggplot2)
library(dplyr)

步骤2：拟合混合效应模型 #

根据数据生成逻辑，将ID设为随机截距，固定效应包含TRTYEAR、N、Myc的三阶交互项：

# 拟合模型
mixed_model <- lmer(r ~ TRTYEAR * N * Myc + (1|ID), data = tempEf)

# 查看模型摘要
summary(mixed_model)

步骤3：生成预测数据及95%置信区间 #

创建包含所有自变量组合的新数据集，基于模型生成群体水平的预测值和置信区间：

# 创建预测用的新数据集：覆盖所有N、Myc组合，TRTYEAR取连续序列值
new_data <- expand.grid(
  TRTYEAR = seq(min(tempEf$TRTYEAR), max(tempEf$TRTYEAR), length.out = 100),
  N = unique(tempEf$N),
  Myc = unique(tempEf$Myc),
  ID = unique(tempEf$ID)[1]  # 固定ID，仅计算固定效应的预测结果
)

# 生成预测值及标准误，忽略随机效应
predictions <- predict(mixed_model, newdata = new_data, re.form = ~0, se.fit = TRUE)
new_data$pred_r <- predictions$fit
# 计算95%置信区间
new_data$lower <- new_data$pred_r - 1.96 * predictions$se.fit
new_data$upper <- new_data$pred_r + 1.96 * predictions$se.fit

步骤4：绘制拟合结果图 #

用ggplot2绘制预测线，并添加95%置信区间的填充区域：

ggplot(new_data, aes(x = TRTYEAR, y = pred_r, color = N, linetype = Myc)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = N), alpha = 0.2, color = NA) +
  labs(
    x = "处理年份",
    y = "响应变量r",
    color = "氮水平",
    linetype = "菌根类型",
    fill = "氮水平"
  ) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    panel.grid = element_blank()
  )

补充：添加原始数据散点（可选） #

如果需要结合原始数据展示拟合效果，可以在图中加入散点：

ggplot() +
  geom_point(data = tempEf, aes(x = TRTYEAR, y = r, color = N), alpha = 0.3) +
  geom_line(data = new_data, aes(x = TRTYEAR, y = pred_r, color = N, linetype = Myc), linewidth = 1) +
  geom_ribbon(data = new_data, aes(x = TRTYEAR, ymin = lower, ymax = upper, fill = N), alpha = 0.2, color = NA) +
  labs(
    x = "处理年份",
    y = "响应变量r",
    color = "氮水平",
    linetype = "菌根类型",
    fill = "氮水平"
  ) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    panel.grid = element_blank()
  )

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