如何用ggplot2绘制glmer模型的拟合值与观测值对比图

1. 数据准备与模型构建 #

先构建与你的数据结构匹配的示例数据，再完成二项式混合效应模型的搭建：

# 加载依赖包
library(lme4)
library(ggplot2)
library(dplyr)

# 模拟符合需求的数据集：处理组、存活/死亡标记、培养皿ID
set.seed(123)
dat <- expand.grid(Treatment = c("Control", "Low", "Medium", "High"),
                   Dish_ID = paste0("Dish_", 1:10)) %>%
  mutate(Ind_Count = sample(10:20, nrow(.), replace = TRUE)) %>%
  uncount(Ind_Count) %>%
  mutate(Survived = rbinom(nrow(.), 1, plogis(-0.5 + as.integer(Treatment)*0.3)))

# 构建混合效应模型
model <- glmer(Survived ~ Treatment + (1 | Dish_ID), data = dat, family = binomial)

2. 计算观测死亡率 #

按处理组汇总，计算每组实际观测到的死亡率（你的数据中0代表死亡，因此死亡率为1 - 存活占比）：

obs_mortality <- dat %>%
  group_by(Treatment) %>%
  summarise(
    mortality = 1 - mean(Survived),
    .groups = "drop"
  )

3. 生成模型拟合预测数据 #

创建覆盖所有处理组的预测数据集，提取模型的固定效应拟合值（群体水平）并转换为死亡率，同时计算置信区间：

# 生成预测用基础数据集
pred_dat <- expand.grid(Treatment = unique(dat$Treatment))

# 获取固定效应拟合的存活概率，转换为死亡率
preds <- predict(model, newdata = pred_dat, type = "response", re.form = NA)
pred_dat$mortality_fit <- 1 - preds

# 用正态近似计算95%置信区间
se <- predict(model, newdata = pred_dat, type = "response", se.fit = TRUE, re.form = NA)$se.fit
pred_dat$lower <- 1 - (preds + 1.96*se)
pred_dat$upper <- 1 - (preds - 1.96*se)

注：re.form = NA表示忽略随机效应，只计算群体水平的拟合线；如果需要展示每个培养皿的个体拟合，可去掉该参数。

4. ggplot2绘图组合 #

将观测死亡率红点、模型拟合线及置信区间整合到同一张图中：

ggplot() +
  # 绘制模型拟合线与置信区间
  geom_line(data = pred_dat, aes(x = Treatment, y = mortality_fit, group = 1), 
            color = "darkblue", linewidth = 1) +
  geom_ribbon(data = pred_dat, aes(x = Treatment, ymin = lower, ymax = upper), 
              alpha = 0.2, fill = "lightblue") +
  # 绘制观测死亡率红点
  geom_point(data = obs_mortality, aes(x = Treatment, y = mortality), 
             color = "red", size = 3, shape = 19) +
  # 调整图表标签与主题
  labs(x = "处理组", y = "死亡率", title = "死亡率 vs 处理组：模型拟合与观测值") +
  theme_bw()

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