plot_model plots conditional model of glmmTMB zero-inflated models when type = "pred"

[grantdadams]

When using plot_model(..., type = "pred") to plot zero-inflated models fit using glmmTMB, plot_model plots the prediction from the conditional model rather than both the conditional and zero-inflated model. Doing a quick test with a poisson and zero-inflated poisson seems to confirm this:

library(glmmTMB)
library(ggplot2)
library(sjPlot)
library(ggeffects)

# Fit two models. A poisson (m1) and a zero-inflated poisson (m2)
data("Salamanders")
m1 <- glmmTMB(count ~ DOP,
              family=poisson, data=Salamanders)

m2 <- glmmTMB(count ~ DOP,
               zi=~DOP,
               family=poisson, data=Salamanders)


# Create a data.frame over the range of the independent variable to predict over
pred_data <- data.frame(matrix(nrow=1000,ncol=1))
colnames(pred_data)<-c("DOP")
pred_data$DOP <- seq(from = min(Salamanders$DOP, na.rm = T), to = max(Salamanders$DOP, na.rm = T), length.out = 1000)

# Predict type = "response" using glmmTMB.predict
count_hat1 <- predict(m1, newdata=pred_data, type = "response", se.fit = TRUE)
count_hat2 <- predict(m2, newdata=pred_data, type = "response", se.fit = TRUE)

# Predict type = "conditional" using glmmTMB.predict
cond_hat1 <- predict(m1, newdata=pred_data, type = "conditional", se.fit = TRUE)
cond_hat2 <- predict(m2, newdata=pred_data, type = "conditional", se.fit = TRUE)


# Assign predictions and 95% confidence interval of response to data.frame
pred_data$fit1 <- count_hat1$fit; pred_data$lower1 = count_hat1$fit - 1.96 * count_hat1$se.fit ; pred_data$upper1 = count_hat1$fit + 1.96 * count_hat1$se.fit
pred_data$fit2 <- count_hat2$fit; pred_data$lower2 = count_hat2$fit - 1.96 * count_hat2$se.fit ; pred_data$upper2 = count_hat2$fit + 1.96 * count_hat2$se.fit

# Assign predictions and 95% confidence interval of conditional to data.frame
pred_data$cond1 <- cond_hat1$fit; pred_data$lower_cond1 = cond_hat1$fit - 1.96 * cond_hat1$se.fit ; pred_data$upper_cond1 = cond_hat1$fit + 1.96 * cond_hat1$se.fit
pred_data$cond2 <- cond_hat2$fit; pred_data$lower_cond2 = cond_hat2$fit - 1.96 * cond_hat2$se.fit ; pred_data$upper_cond2 = cond_hat2$fit + 1.96 * cond_hat2$se.fit

# Plot response using ggplot 
plot1<- pred_data%>%ggplot(aes(y = fit1,x = DOP))+
  geom_ribbon(aes(ymin = lower1,ymax=upper1),fill="grey",alpha=0.6)+
  geom_line(color="black",lwd=1)
plot2<- pred_data%>%ggplot(aes(y = fit2,x = DOP))+
  geom_ribbon(aes(ymin = lower2,ymax=upper2),fill="grey",alpha=0.6)+
  geom_line(color="black",lwd=1)

# Plot conditional using ggplot 
plot_cond1<- pred_data%>%ggplot(aes(y = cond1,x = DOP))+
  geom_ribbon(aes(ymin = lower_cond1,ymax=upper_cond1),fill="grey",alpha=0.6)+geom_line(color="black",lwd=1)
plot_cond2<- pred_data%>%ggplot(aes(y = cond2,x = DOP))+
  geom_ribbon(aes(ymin = lower_cond2,ymax=upper_cond2),fill="grey",alpha=0.6)+geom_line(color="black",lwd=1)

# Compare with plot_mode
plot1
plot_cond1
plot_model(m1, type = "pred")

plot2
plot_cond2
plot_model(m2, type = "pred")