如何解决glmmTMB Hurdle模型二元与计数组分绘图重复问题？

问题背景 #

使用glmmTMB拟合Hurdle混合模型，分别通过二元组分（zi，对应感染流行率的零/非零数据）和零截断负二项组分（cond，对应感染强度的计数数据）分析响应变量对环境变化的影响。模型拟合、固定效应显著性检验（car::Anova）、组间成对比较（emmeans）均能得到两组分的独立结果，但绘制预测图时，即使指定了component参数，两组分的图完全一致，与成对比较的差异结果矛盾。

问题原因 #

查看示例代码后，核心错误出在二元组分（zi）绘图的误差棒数据引用：

• 绘制zi组分预测图时，geom_errorbar的data参数错误地使用了predict_nest.nb（cond组分的预测结果），而非predict_nest.zi，导致误差棒与点的来源数据不匹配，最终两张图呈现出相同的视觉效果。

解决办法 #

1. 修正数据引用错误：确保绘图时每个组分的所有图层都使用对应组分的预测数据。
2. 优化绘图代码结构：将公共美学映射（aes）统一放在ggplot()函数中，减少重复代码，提升可读性。
3. 添加组分标识：为每张图添加标题或坐标轴标签，明确区分二元组分（流行率）和计数组分（感染强度）。
4. 可选：合并两组分图：使用分面（facet）将两组分的预测图合并展示，更便于对比差异。

修正后的代码示例 #

1. 单独绘制两组分预测图 #

library("glmmTMB") 
library("emmeans")
library("ggeffects")
library("ggplot2")

data("Owls")

# 拟合模型
mod <- glmmTMB(SiblingNegotiation ~ FoodTreatment + Nest + offset(log(BroodSize)), 
               zi= ~ FoodTreatment + Nest, 
               family=truncated_nbinom2, data=Owls)

# 生成两组分的预测数据
predict_nest.nb <- ggeffect(mod, "Nest", component="cond")
predict_nest.zi <- ggeffect(mod, "Nest", component="zi")

# 绘制计数组分（cond）预测图
ggplot(predict_nest.nb, aes(x=x, y=predicted, color=x)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1, size=1) +
  geom_point(size=4) +
  labs(title="计数组分（感染强度）预测值", x="巢（Nest）", y="预测谈判次数") +
  theme(axis.text.x=element_text(size=10, angle=90, vjust=1, hjust=1))

# 绘制二元组分（zi）预测图
ggplot(predict_nest.zi, aes(x=x, y=predicted, color=x)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1, size=1) +
  geom_point(size=4) +
  labs(title="二元组分（感染流行率）预测值", x="巢（Nest）", y="非零概率") +
  theme(axis.text.x=element_text(size=10, angle=90, vjust=1, hjust=1))

2. 合并两组分预测图（分面展示） #

# 合并两组分数据，添加组分标识
predict_combined <- rbind(
  cbind(predict_nest.nb, component="计数组分（感染强度）"),
  cbind(predict_nest.zi, component="二元组分（感染流行率）")
)

# 分面绘图
ggplot(predict_combined, aes(x=x, y=predicted, color=x)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1, size=1) +
  geom_point(size=4) +
  facet_wrap(~component, scales="free_y") +
  labs(x="巢（Nest）", y="预测值") +
  theme(axis.text.x=element_text(size=10, angle=90, vjust=1, hjust=1))

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