R语言lme4混合模型ggplot分面绘制分组拟合线置信区间

使用lme4包拟合线性混合效应模型后，需要在ggplot分面图中为每个Direction组的拟合线添加匹配组内数据范围的置信区间。

数据结构 #

使用的数据框包含以下字段：

• Plot_label：字符型，样方标识
• PD_avg：数值型，响应变量
• Year：因子型，年份随机效应
• GS_Prec：数值型，生长季降水量预测变量
• Direction：因子型，分组变量

原有代码问题 #

初始代码设置se=T后置信区间未正常显示，使用emmeans方法绘制的置信区间范围与分组散点范围不匹配：

# 初始拟合与绘图代码
mixed.lm <- lmer(PD_avg ~ log(GS_Prec) * Direction + (1|Plot_label) + (1|Year), data = data, REML=TRUE)
pred1 <- predict(mixed.lm, newdata = data, re.form = NA) 

ggplot(data, aes(log(GS_Prec), PD_avg, colour = Direction)) +
  geom_point(alpha = .2) +
  facet_wrap(~Direction) +
  geom_smooth(aes(y = pred1, colour = Direction), method = "lm", size = 1.5, se = T)

# emmeans尝试的无效代码
gr <- ref_grid(mixed.lm, cov.keep = c("GS_Prec", "Direction"))
emm <- emmeans(gr, spec = c("GS_Prec","Direction"), level = 0.95)

ggplot(data, aes(log(GS_Prec), PD_avg, colour = Direction)) +
  geom_point(alpha = .2) +
  facet_wrap(~Direction) +
  geom_smooth(aes(y = pred1, colour = Direction), method = "lm", size = 1.5) +
  geom_ribbon(data = data.frame(emm), aes(ymin = lower.CL, ymax = upper.CL, y = NULL, fill = Direction), alpha = 0.1)+
  geom_smooth(aes(y = pred1, colour = Direction), method = "lm", size = 1.5)

测试用子集数据：

data.1 <- data.frame(Plot_label = c("BT 1-1-3", "BT 1-1-3", "BT 1-2-1", "BT 1-2-1",
                                    "GW 1-1-1", "GW 1-1-1", "GW 1-5-2", "GW 1-5-2",
                                    "SP 1-5-2", "SP 1-5-2", "SP 2-8-2", "SP 2-8-2"),
                     PD_avg = c("1196.61", "1323.15", "1172.17", "757.18",
                                "1516.02", "801.87", "1422.93", "1062.10",
                                "1580.51", "1520.30", "1326.25", "1321.89"),
                     Year = c("2016", "2017", "2016", 2017,
                              "2016", "2017", "2016", "2017",
                              "2016", "2017", "2016", "2017"),
                     Direction = c("BT-BT", "BT-BT", "BT-BT", "BT-BT",
                                   "GW-BT", "GW-BT", "GW-BT", "GW-BT",
                                   "SP-SP", "SP-SP", "SP-SP", "SP-SP"),
                     GS_Prec = c("130.5", "190.5", "130.5", "190.5",
                                 "130.5", "190.5", "130.5", "190.5",
                                 "593.26", "480.29", "593.26", "593.26"))

错误原因 #

1. geom_smooth(se=T)不生效：传入提前预测的pred1后，geom_smooth会对x和传入的y重新拟合普通线性模型，不会自动基于混合效应模型计算置信区间，重复调用还会导致拟合线叠加。
2. emmeans结果错位：一是没有按分组生成对应x轴取值范围内的预测网格，导致预测点覆盖范围不对；二是x轴绘图用的是log(GS_Prec)，但emmeans输出的预测值基于原始尺度的GS_Prec，尺度不匹配导致位置偏移。

正确实现代码 #

核心思路是手动生成分组预测网格，基于模型输出的标准误计算置信区间，直接用geom_line和geom_ribbon绘制拟合线和置信区间，不需要geom_smooth做二次拟合。

library(tidyverse)
library(lme4)

# 先清洗数据，将误存为字符型的变量转为对应格式
data <- data.1 %>% 
  mutate(
    PD_avg = as.numeric(PD_avg),
    GS_Prec = as.numeric(GS_Prec),
    Year = as.factor(Year)
  )

# 拟合混合效应模型
mixed.lm <- lmer(PD_avg ~ log(GS_Prec) * Direction + (1|Plot_label) + (1|Year), 
                 data = data, REML=TRUE)

# 按分组构建预测网格：每个Direction组内取GS_Prec从最小值到最大值的100个等距点
pred_grid <- data %>% 
  group_by(Direction) %>% 
  reframe(GS_Prec = seq(min(GS_Prec), max(GS_Prec), length.out = 100))

# 预测固定效应拟合值与标准误，计算95%置信区间
pred_res <- predict(mixed.lm, newdata = pred_grid, re.form = NA, se.fit = TRUE)
pred_grid <- pred_grid %>% 
  mutate(
    fit = pred_res$fit,
    lwr = fit - 1.96 * pred_res$se.fit,
    upr = fit + 1.96 * pred_res$se.fit
  )

# 绘图
ggplot(data, aes(x = log(GS_Prec))) +
  geom_point(aes(y = PD_avg, colour = Direction), alpha = 0.2) +
  facet_wrap(~Direction) +
  # 绘制置信区间
  geom_ribbon(data = pred_grid, aes(ymin = lwr, ymax = upr, fill = Direction), alpha = 0.1) +
  # 绘制拟合线
  geom_line(data = pred_grid, aes(y = fit, colour = Direction), linewidth = 1.5) +
  labs(y = "PD_avg")

补充说明：大样本下用1.96作为95%置信区间的临界值足够准确，小样本可以替换为对应自由度的t分布临界值。如果需要考虑随机效应的不确定性，可以用lme4::bootMer做bootstrap抽样计算置信区间，结果会更保守。

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