如何用glmmTMB计算负二项混合效应模型的预测值？

关键要点 #

• 访花序列设为因子型完全可行，predict()函数对因子变量支持良好，无需担心兼容性
• 加随机效应后预测报错，核心原因是预测数据集格式不匹配，或未明确随机效应的处理方式
• 同图展示多体型蜜蜂结果，只需构造包含所有变量组合的预测数据集，再用ggplot2直接绘图

分步实现代码 #

1. 构建正确的nbinom1混合模型 #

假设你的数据集为bee_df，包含字段：pollen（花粉粒数）、visit_seq（访花序列，已转因子）、bee_type（蜜蜂体型，如"large"/"small"）、trial（实验轮次，随机效应）：

library(glmmTMB)
# 模型包含体型×序列的交互项，轮次作为随机截距
mod <- glmmTMB(pollen ~ bee_type * visit_seq + (1|trial),
               data = bee_df,
               family = nbinom1)

2. 生成预测数据集 #

构造覆盖所有体型和访花序列组合的数据集，同时指定随机效应的处理（这里选择忽略随机效应，聚焦固定效应的趋势；若要保留随机效应，需指定具体轮次）：

# 生成所有变量组合的预测框架
pred_df <- expand.grid(
  bee_type = unique(bee_df$bee_type),
  visit_seq = levels(bee_df$visit_seq),
  trial = unique(bee_df$trial)[1]  # 选单个轮次，或用re.form=~0忽略随机效应
)

# 生成响应尺度的预测值
pred_df$pred_pollen <- predict(mod, newdata = pred_df, re.form = ~0, type = "response")

3. 同图绘制双体型结果 #

用ggplot2直接实现曲线对比，无需依赖ggeffects/ggpredict：

library(ggplot2)
ggplot(pred_df, aes(x = visit_seq, y = pred_pollen, color = bee_type, group = bee_type)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  labs(x = "访花序列", y = "预测花粉粒数量", color = "蜜蜂体型") +
  theme_classic()

4. 常见问题排查 #

• 随机效应报错：检查trial是否为因子型，预测数据中的trial必须是原模型中存在的水平
• 曲线异常：确认visit_seq已转为因子（用bee_df$visit_seq <- as.factor(bee_df$visit_seq)），且模型包含了体型与序列的交互项（bee_type * visit_seq），否则无法捕捉不同体型的序列趋势差异

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