Plot effects of glmmTMB separately for zero-truncated count component and for binary component

119 views Asked by At

I fit a Hurdle mixed model (glmmTMB function in glmmTMB package) to simultaneously explore how infection prevalence (binary part of the model, zeros and non-zeros data) and infection intensity (zero-truncated negative binomial model, count data) respond to environmental changes. I fitted the model, used Anova for significant fixed effects, and then emmeans and contrasts for pairwise comparisons. For each step, I was able to get separate outputs for the binary component and for the count component. How do I get separate plots, based on the model, for these two components?

As an example, I use the Owls data.

library("glmmTMB") 
library("emmeans")

data("Owls")

# fit model
mod <- glmmTMB(SiblingNegotiation ~ FoodTreatment + Nest + offset(log(BroodSize)), 
                       zi= ~ FoodTreatment + Nest, 
                       family=truncated_nbinom2, data=Owls)

# Anova for significant fixed effects
#count part
car::Anova(mod, type="II", component="cond") 
#binary part
car::Anova(mod, type="II", component="zi")

# pairwise comparisons
#binary part
emm.nest.zi <- emmeans::emmeans(mod, c("Nest"), component="zi")
contrast(emm.nest.zi, alpha=0.05, method="pairwise", adjust="bh")
#count part
emm.nest.nb <- emmeans::emmeans(mod, c("Nest"), component="cond")
contrast(emm.nest.nb, alpha=0.05, method="pairwise", adjust="bh")

Now, I would like to get two model-based plots. One for the binary part of the model and another for the count part of the model. This is my attempt.

# plot the count part of the glmmTMB model
predict_nest.nb <- ggeffect(mod,"Nest", component="cond")
ggplot(predict_nest.nb,aes(x=x, y=predicted), fill=group, group=group, color=group)+
  geom_errorbar(data=predict_nest.nb, 
                mapping=aes(x=x, ymin=conf.low, ymax=conf.high), 
                width=0.1, size=1) +
  geom_point(size=4, aes()) +
  theme(axis.text.x=element_text(size=10, angle=90, vjust=1, hjust=1, face="plain"))

# plot the binary part of the glmmTMB model
predict_nest.zi <- ggeffect(mod,"Nest", component="zi")
ggplot(predict_nest.zi,aes(x=x, y=predicted), fill=group, group=group, color=group)+
  geom_errorbar(data=predict_nest.nb, 
                mapping=aes(x=x, ymin=conf.low, ymax=conf.high), 
                width=0.1, size=1) +
  geom_point(size=4, aes()) +
  theme(axis.text.x=element_text(size=10, angle=90, vjust=1, hjust=1, face="plain"))

But I get two identical plots, despite having specified the component="zi" or "cond". In my original df, pairwise contrasts show some significant differences in the zero-part, that are not the same as in the count-part of the model. So, the two plots should look different.

Why it doesn't work? How can I plot separately the two components of the glmmTMB model?

Thank you for any input =)

1

There are 1 answers

2
Russ Lenth On

The plotting functions in the emmeans package itself do this kind of thing pretty well...

require("emmeans")

plot(pairs(emm.nest.zi))
plot(pairs(emm.nest.nb))