Overlay of forest plot from ZINB model

Question

I want an overlay of a forest plot from the ZINB models of full and the subset of data using the sjPlot package. As you may know, the ZINB model produces two models: one for the count model and one for the zero-inflated model. plot_model works fine when employing the ZINB model from either full or a subset of data meaning producing a plot for both models (count and zero models), but when I overlay using plot_models then only one plot is produced for the count model. I am looking for the count and zero-inflated model plots from the full and sub-model for both the full and the subset of data. any help would be much appreciated

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


zinb_all_uni <- zeroinfl(ivdays~age,
                 link="logit",
                 dist = "negbin",
                 data=caterpillor)



summary(zinb_all_uni)
plot_model(zinb_all_uni, type="est")


zinb_full_adj <- zeroinfl(ivdays~age+sex+edu,
                 link="logit",
                 dist = "negbin",
                 data=caterpillor)

summary(zinb_full_adj)
plot_model(zinb_full_adj, type="est", terms = c("count_ageb", "count_agec", "zero_ageb", "zero_agec"))

############ second model#######


Zinb_uni_sub <- zeroinfl(ivdays~age,
                     link="logit",
                     dist = "negbin",
                     data=subset(caterpillor, country=="eng"))



summary(zinb_uni_sub)
plot_model(zinb_uni_sub, type="est")



zinb_adj_sub <- zeroinfl(ivdays~age+sex+edu,
                 link="logit",
                 dist = "negbin",
                 data=subset(caterpillor, country=="eng"))

summary(zinb_adj_sub)

plot_model(zinb_adj_sub, type="est", terms = c("count_ageb", "count_agec", "zero_ageb", "zero_agec"))



### overlying plots from both models


plot_models(zinb_all_uni, Zinb_uni_sub)

plot_models(zinb_full_adj, zinb_adj_sub)

DATA:

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, 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))

but when i overlay plots i get only one plot

Answer 1

Code below, basic points:

• when I run into trouble with automated machinery like plot_model I usually prefer to use machinery like broom::tidy() (for coefficients) or the ggeffects or emmeans packages (for predictions) and build my own ggplot — for me, it's easier than trying to figure out what the more automated tool is doing
• broom doesn't have a tidy() method for zeroinfl models, but a little googling finds one in the poissonreg package ...
• ... however, that tidy() method doesn't have machinery for constructing confidence intervals or back-transforming coefficients to a count-ratio or odds-ratio scale, so I had to implement my own below ...

library(broom)
library(poissonreg)
library(tidyverse) ## purrr::map_dfr, ggplot ...
theme_set(theme_bw())
library(colorspace)
mod_list <- list(all_uni = zinb_all_uni, uni_sub = Zinb_uni_sub,
     full_adj = zinb_full_adj, adj_sub = zinb_adj_sub)


tidy(zinb_all_uni, type = "all")

coefs <- (mod_list
    |> map_dfr(tidy, type = "all",
               .id = "model")
    ## construct CIs
    |> mutate(conf.low  = qnorm(0.025, estimate, std.error),
              conf.high = qnorm(0.975, estimate, std.error))
    |> filter(term != "(Intercept)")  ## usually don't want this
    ## cosmetic (strip results down to the components we actually need)
    |> select(model, term, type, estimate, conf.low, conf.high)
    ## back-transform
    |> mutate(across(c(estimate, conf.low, conf.high), exp))
)

ggplot(coefs, aes(x = estimate, y = term, colour = model)) +
    geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
                    position = position_dodge(width = 0.5)) +
    ## separate count-ratio and odds-ratio (conditional/zero) plots
    facet_wrap(~type, scale = "free") +
    scale_color_discrete_qualitative() ## cosmetic

If you only want to see the age-related coefficients you can add

 |> filter(stringr::str_detect(term, "^age"))

to the end of the pipeline that defines coefs.