使用sjPlot叠加ZINB模型森林图：计数与零膨胀模型显示异常

使用sjPlot包绘制全数据集与子集数据的零膨胀负二项（ZINB）模型叠加森林图时，单独调用plot_model()能正常展示ZINB模型的计数模型和零膨胀模型两部分结果，但使用plot_models()进行叠加操作时，仅能显示计数模型的对比图，无法同时呈现两个模型部分的叠加效果。

plot_models()函数默认仅提取ZINB模型的计数部分系数，要实现计数与零膨胀模型的同时叠加展示，可通过手动提取两个模型的系数数据，再用ggplot2自定义绘制森林图。以下是完整实现步骤：

1. 加载依赖包与数据 #

library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggplot2)
library(MASS)
library(pscl)
library(boot)

# 加载示例数据集
caterpillor=structure(list(id = 1:100,
               age = structure(c(1L, 1L, 2L, 1L, 
                                 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 
                                 1L, 1L, 2L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 2L, 1L, 2L, 2L, 
                                 2L, 3L, 3L, 3L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 
                                 2L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 3L, 
                                 3L, 3L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 2L, 1L, 
                                 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 3L, 3L, 3L),
                               .Label = c("a", "b", "c"), class = "factor"),
               sex = structure(c(2L, 
                                 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 
                                 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 
                                 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 
                                 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 
                                 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 
                                 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 
                                 2L, 1L, 1L),
                               .Label = c("F", "M"), class = "factor"),
               country = structure(c(1L, 
                                     1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 
                                     2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 
                                     1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 
                                     2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 
                                     3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 
                                 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 
                                 2L, 2L, 2L),
                                   .Label = c("eng", "scot", "wale"), class = "factor"), 
               edu = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 
                                 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 
                                 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 
                                 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 
                                 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 
                                 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 
                                 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L),
                               .Label = c("x", "y", "z"), class = "factor"),
               lungfunction = c(45L, 
                                23L, 25L, 45L, 70L, 69L, 90L, 50L, 62L, 45L, 23L, 25L, 45L, 
                                70L, 69L, 90L, 50L, 62L, 45L, 23L, 25L, 45L, 70L, 69L, 90L, 
                                50L, 62L, 45L, 23L, 25L, 45L, 70L, 69L, 90L, 50L, 62L, 45L, 
                                23L, 25L, 45L, 70L, 69L, 90L, 50L, 62L, 45L, 23L, 25L, 45L, 
                                70L, 69L, 90L, 50L, 62L, 45L, 23L, 25L, 45L, 70L, 69L, 90L, 
                                50L, 62L, 45L, 23L, 25L, 45L, 70L, 69L, 90L, 50L, 62L, 45L, 
                                23L, 25L, 45L, 70L, 69L, 90L, 50L, 62L, 25L, 45L, 70L, 69L, 
                                90L, 50L, 62L, 25L, 45L, 70L, 69L, 90L, 50L, 62L, 25L, 45L, 
                                70L, 69L, 90L),
               ivdays = c(15L, 26L, 36L, 34L, 2L, 4L, 5L, 
                          8L, 9L, 15L, 26L, 36L, 34L, 2L, 4L, 5L, 8L, 9L, 15L, 26L, 
                          0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
                          0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
                          0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
                          0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
                          0L, 0L, 0L, 0L, 0L, 5L, 8L, 9L, 36L, 34L, 2L, 4L, 5L, 8L, 
                          9L, 36L, 34L, 2L, 4L, 5L),
               no2_quintile = structure(c(1L, 
                                          1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                          1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                          2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                          3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                          3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 
                                          4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 
                                          5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L),
                                        .Label = c("q1", "q2", 
                                                   "q3", "q4", "q5"), class = "factor")),
          class = "data.frame", row.names = c(NA, 
                                              -100L))

2. 拟合ZINB模型 #

# 全数据集多变量模型
zinb_full_adj <- zeroinfl(ivdays~age+sex+edu, link="logit", dist = "negbin", data=caterpillor)
# 子集（eng）多变量模型
zinb_adj_sub <- zeroinfl(ivdays~age+sex+edu, link="logit", dist = "negbin", data=subset(caterpillor, country=="eng"))

3. 提取并整理系数数据 #

使用get_model_data()分别提取两个模型的计数和零膨胀部分系数，合并后统一格式：

# 提取全数据集模型的计数与零膨胀系数
df_full_count <- get_model_data(zinb_full_adj, type = "est", component = "count")
df_full_count$model <- "全数据集"
df_full_count$component <- "计数模型"

df_full_zero <- get_model_data(zinb_full_adj, type = "est", component = "zero")
df_full_zero$model <- "全数据集"
df_full_zero$component <- "零膨胀模型"

# 提取子集模型的计数与零膨胀系数
df_sub_count <- get_model_data(zinb_adj_sub, type = "est", component = "count")
df_sub_count$model <- "子集（eng）"
df_sub_count$component <- "计数模型"

df_sub_zero <- get_model_data(zinb_adj_sub, type = "est", component = "zero")
df_sub_zero$model <- "子集（eng）"
df_sub_zero$component <- "零膨胀模型"

# 合并所有数据并清理变量名
plot_data <- rbind(df_full_count, df_full_zero, df_sub_count, df_sub_zero)
plot_data$term <- gsub("count_|zero_", "", plot_data$term)

4. 绘制叠加森林图 #

通过ggplot2实现计数与零膨胀模型的分面叠加展示：

ggplot(plot_data, aes(x = estimate, y = term, color = model, shape = model)) +
  # 绘制置信区间
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2, position = position_dodge(width = 0.5)) +
  # 绘制系数点
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  # 添加系数为0的参考线
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  # 按模型组分面展示
  facet_wrap(~component, scales = "free_y") +
  # 设置标签与主题
  labs(x = "回归系数（95%置信区间）", y = "变量", color = "数据集", shape = "数据集") +
  theme_bw() +
  theme(
    panel.grid.major.y = element_blank(),
    legend