负二项GLMM显著系数与对数优势比/发生率比绘图方法问询

问题核心 #

运行负二项广义线性混合模型（glmer.nb）后，需绘制类似目标样式的对数尺度系数图，但遇到以下问题：

• 选择单个变量（如Predator、Aspect）绘图时，缺失参考类别（如Predator的No、Aspect的NE）
• 部分绘图工具输出的Y轴为响应变量Number，而非对数系数
• 不想显示截距，因其包含其他变量的综合信息

模拟数据与模型代码 #

Number <- c(25,16,16,13,12,9,15,25,4,5,230,259,224,259,588,604,583,576,599,664)
Location <- c("Borakolalo","Borakolalo","Borakolalo","Borakolalo","Bloemhof","Bloemhof","Bloemhof",   
              "Bloemhof","Boskop","Boskop","Boskop","Boskop","Kgaswane","Kgaswane","Kgaswane",   
              "Kgaswane","Mafikeng","Mafikeng","Mafikeng","Mafikeng")
Nitrogen<-c(0.0889,0.0406,0.0835,0.0737,0.0578,0.0806,0.0914,0.09630,0.0718,0.08955,1.0211,1.9489,
            1.9808,1.0053,1.9682,1.9794,1.0959,1.0028,1.9281,1.9887)
Dist_water<- c(2156.0,3783.8,3285.8,2574.7,2242.3,2729.5,2418.1,2874.9,2869.0,2563.0,257.1,660.4,
               440.4,417.7,562.6,528.5,426.5,591.2,435.9,306.5)
Predator<-c("Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","No","No","No","No",
            "No","No","No","No","No","No")
Aspect<-c("SE","S","S","NE","NW","S","SW","SE","S","S","NE",
          "S","S","SE","SE","SE","SE","SE","SE","SE")
bles<-data.frame(Number,Location,Nitrogen,Dist_water,Predator,Aspect)

bles_sc <- transform(bles,
                     Nitrogen = drop(scale(Nitrogen)),
                     Dist_water = drop(scale(Dist_water)))

library(lme4)
mod<-glmer.nb(Number~Nitrogen + Dist_water + Predator+Aspect+(1|Location), data=bles_sc)

summary(mod)

解决方案：用broom.mixed提取系数+ggplot2自定义绘图 #

通过手动提取系数并补充参考类别信息，可完全控制绘图内容，匹配目标样式：

1. 加载所需包：

library(broom.mixed)
library(ggplot2)
library(dplyr)

2. 提取模型系数及置信区间，整理数据：

# 提取固定效应系数
coef_df <- tidy(mod, effects = "fixed", conf.int = TRUE) %>%
  # 移除截距
  filter(term != "(Intercept)") %>%
  # 拆分变量名和类别
  separate(term, into = c("variable", "level"), sep = "(?<=[a-zA-Z])(?=[A-Z])", fill = "right") %>%
  # 处理连续变量的level值
  mutate(level = ifelse(is.na(level), variable, level))

# 手动添加参考类别（系数为0，置信区间也为0）
ref_cats <- tibble(
  variable = c("Predator", "Aspect"),
  level = c("No", "NE"),
  estimate = 0,
  conf.low = 0,
  conf.high = 0,
  std.error = NA,
  statistic = NA,
  p.value = NA
)

# 合并系数数据和参考类别
full_coef_df <- bind_rows(coef_df, ref_cats) %>%
  # 按变量和类别排序，确保参考类别在前面
  arrange(variable, level)

3. 绘制目标样式的系数图：

ggplot(full_coef_df, aes(x = estimate, y = level, color = variable)) +
  # 绘制置信区间
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2) +
  # 绘制系数点
  geom_point(size = 3) +
  # 添加x轴参考线（系数为0）
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  # 按变量分面，每个变量单独一列
  facet_wrap(~variable, scales = "free_y", ncol = 1) +
  # 设置坐标轴标签
  labs(x = "对数系数", y = "", color = "变量") +
  # 调整主题，匹配目标图简洁风格
  theme_bw() +
  theme(
    panel.grid.major.y = element_blank(),
    legend.position = "none"
  )

关键说明 #

• 补充参考类别：R线性模型默认不输出参考类别系数，手动添加后可完整展示所有类别对比
• 自定义绘图：ggplot2可完全控制绘图元素，避免其他工具的默认限制（如Y轴显示错误、类别缺失）
• 对数尺度：直接使用模型输出的estimate（对数尺度系数，对应负二项模型的对数均值）

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