如何绘图展示Negative Binomial GLMM分类变量与计数结果及解读

一、cat_plot绘图方式的问题 #

你当前的绘图核心问题是混淆了尺度：cat_plot默认用link尺度（对数尺度）展示固定效应原始系数，但你想呈现的是发生率比（IRR，即exp(系数)），这就导致数值看起来“过高”——本质是把对数尺度的系数直接当成IRR展示，属于尺度错误。

正确的绘图方式需要明确指定转换为指数尺度：

• 在cat_plot中添加transform = "exp"参数，Y轴就会显示真实的IRR值，而非对数尺度的原始系数。
• 关于“无显著性”：cat_plot默认不会自动标注显著性，需添加show.p = TRUE参数来显示固定效应的p值；如果要展示组间两两比较的显著性，得先通过emmeans包做事后检验，再手动把显著性标记加到图中。

二、模型结果解读（负二项GLMM） #

你的模型针对过度离散的有蹄类计数数据选择负二项分布，加入位点随机截距(1|Location)控制位点间的非独立性，核心解读逻辑如下：

1. 固定效应（Aspect、Predator）
   • 所有结果要基于**IRR（exp(模型系数)）**解读：IRR=1代表该水平与参考水平的计数发生率无差异；IRR>1代表该水平计数发生率是参考水平的对应倍数，IRR<1则表示发生率更低。
   • Aspect某水平无显著性：说明该水平和参考水平的计数发生率差异在统计上不显著（p>0.05），现有数据无法证明坡向会影响有蹄类计数。
   • Predator的Yes水平无显著性：说明捕食者存在与否对有蹄类计数的影响没有统计学意义，没有足够证据支持“捕食者影响数量”的结论。
2. 随机效应(1|Location)
   查看模型输出的Random effects部分，截距的方差越大，说明不同位点之间的有蹄类计数差异越大——这个随机效应的作用是把位点间的变异从残差中分离出来，让固定效应的估计更准确。
3. 模型拟合有效性
   • 对比带固定效应的模型和零模型（仅含随机效应）的AIC/BIC值：如果目标模型的AIC/BIC更低，说明加入Aspect、Predator确实改善了模型拟合。
   • 检查过度离散参数：负二项模型的离散参数若接近1，说明泊松模型也适用；若远大于1，说明选择负二项是合理的。

三、修正后的R代码示例 #

假设你的原始代码如下：

library(lme4)
library(sjPlot)

# 拟合负二项GLMM
nb_model <- glmer.nb(Count ~ Aspect + Predator + (1|Location), data = your_data)

# 原始绘图（存在尺度错误）
cat_plot(nb_model, type = "fe", terms = c("Aspect", "Predator"))

调整为正确的IRR绘图代码：

# 1. 绘制带IRR和显著性的图
cat_plot(nb_model, 
         type = "fe", 
         terms = c("Aspect", "Predator"),
         transform = "exp",  # 转换为IRR尺度
         show.values = TRUE, # 显示IRR数值
         show.p = TRUE,      # 显示固定效应的p值
         vline.color = "red")# 添加IRR=1的参考线（方便对比）

# 2. 做组间两两事后检验（若需细致的显著性对比）
library(emmeans)
# Aspect各水平的两两比较（输出IRR和p值）
aspect_posthoc <- pairs(emmeans(nb_model, ~ Aspect, type = "response"))
print(aspect_posthoc)

# Predator两组的比较
pred_posthoc <- pairs(emmeans(nb_model, ~ Predator, type = "response"))
print(pred_posthoc)

# 3. 把事后检验的显著性标注到图中（以Aspect为例）
library(ggplot2)
library(ggsignif)
p <- cat_plot(nb_model, type = "fe", terms = "Aspect", transform = "exp") +
  geom_signif(comparisons = list(c("你的Aspect水平1", "你的Aspect水平2")),
              annotations = as.character(round(aspect_posthoc$p.value[1], 3)),
              y_position = 你需要的Y轴位置)
print(p)

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