Interpretation of glmmTMB output for zero-inflated negative binomial regression

Question

I fitted, using glmmTMB R package, a zero-inflated negative binomial GLMM, with offset and a random factor, to investigate which variables could explain animal species' range filling. The response variable (Range_filling) is expressed as cell counts, and to account for the total area available I have put its log as an offset with offset(log(POAR)). I also accounted for variation among different species' orders (with (1|Order)). I log-transformed and scaled some predictors because they ranged on a huge scale; I transformed adult_mass_g, Residence_time_d, and Native_range.

This is my model call:

mod <- glmmTMB(Range_filling ~ adult_mass_g + Residence_time_d + Pathways_tot + Native_range + Points_succ_intro +
                             offset(log(POAR)) + (1|Order),
                           data = sp_tbl_scaled,
                           ziformula = ~ 1,
                           family = nbinom1)

And this is my model output:

summary(mod)

Family: nbinom1  ( log )
Formula:          Range_filling ~ adult_mass_g + Residence_time_d + Pathways_tot +  
    Native_range + Points_succ_intro + offset(log(POAR)) + (1 |      Order)
Zero inflation:                 ~1
Data: sp_tbl_scaled

     AIC      BIC   logLik deviance df.resid 
   507.0    523.6   -244.5    489.0       38 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Order  (Intercept) 0.4897   0.6998  
Number of obs: 47, groups:  Order, 5

Dispersion parameter for nbinom1 family (): 85.2 

Conditional model:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -1.91215    0.45947  -4.162 3.16e-05 ***
adult_mass_g      -0.01812    0.23620  -0.077  0.93886    
Residence_time_d  -0.29223    0.13907  -2.101  0.03561 *  
Pathways_tot      -0.64082    0.12099  -5.297 1.18e-07 ***
Native_range      -0.33345    0.11940  -2.793  0.00523 ** 
Points_succ_intro  0.02085    0.00306   6.812 9.64e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.6250     0.7017  -3.741 0.000183 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Now I would like to interpret my results, and specifically, I would like to find which explanatory variables are more important in predicting my response variable. To do so, I am looking at the model output and I also produced some response plots using sjPlot, like this:

Looking specifically at those 3 predictors, I have a couple of questions:

1. from the slope of the red line, it seems that the negative relationship between Range_filling and Residence_time_d is stronger compared to the one of Range_filling with Pathways_tot. If I look at the Estimate, Residence_time_d has a -0.29, while Pathways_tot has a -0.64. As Residence_time_d is log-transformed and scaled, I wonder if this big difference depends on that, and how to successfully "compare" those two estimates (for instance, if I want to say which of the two predictors has a stronger influence on the response variable). I guess I shall look at the p-value, as for Pathways_tot is highly significant, while it has a 0.03 value for Residence_time_d.
2. I specified ziformula = ~1, to have a constant term (intercept) used for the zero-inflation model, assuming that the probability of excess zeros in the response variable is the same across all levels of the predictor variables, regardless of their values. This leaves me with only one line in my zero-inflation model summary, which appears to be significant based on the p-value. I'm not really sure how to interpret this. Moreover, what is zero-inflated is actually my response variable (Range_filling), not my predictors. Should I then specify the model using ziformula = ~ Range_filling? If I do that, the significance of the predictors doesn't change, and I get a p-value close to 1. As in page 383 of glmmTMB documentation, "The zero-inflation model estimates the probability of an extra zero such that a positive contrast indicates a higher chance of absence (e.g. minedno <0 means fewer absences in sites unaffected by mining); this is the opposite of the conditional model where a positive contrast indicates a higher abundance (e.g., minedno >0 means higher abundances in sites unaffected by mining)". If I understand this correctly, then, a p-value close to 1 should be better, but I'm not really sure.

Thanks in advance! Cross-posted on SO

Answer 1

There's a bunch here.

• If the proportion of filled cells in a range is sometimes large (e.g. greater than 30 or 40%) you might want to consider a binomial-type model (e.g. a zero-inflated beta-binomial model) rather than a Poisson-plus-offset model, which becomes unrealistic when the outcome is not rare.
• The simplest way to quantify 'effect size' (in my opinion) is to fit a version of your model with all of the predictor variables scaled to have a standard deviation of 1 (see Schielzeth 2010; there is an lm.beta package that handles this standardization after fitting for linear models, but it won't work for glmmTMB objects - it can be done after model fitting but might be easiest to use scale() on your predictors before you fit the model). I would not compare the importance of the parameters by looking at p-values, as this can be influenced by the uncertainty in the parameter as well as its magnitude.
• The p-value in the zero-inflation term is not really useful (and should probably be eliminated from the model summary); it compares the coefficient (log-odds of zero-inflation) to a reference value of zero, which would correspond to 50% zero-inflation — not really a meaningful value. Testing for zero-inflation is tricky; the simplest (although crude and conservative) way to test the importance of zero-inflation would be to compare AIC values or do a likelihood ratio test (anova()) between a model with and without the zero-inflation component.
• Adding ziformula = ~ Range_filling makes no sense at all (sorry). The existing model (ziformula = ~1) already models zero-inflation in the response variable. The p-value and the zero-inflation probability are completely different; in your case the estimated zero-inflation probability is plogis(-2.6) (approx. 7%)

Schielzeth, Holger. 2010. “Simple Means to Improve the Interpretability of Regression Coefficients: Interpretation of Regression Coefficients.” Methods in Ecology and Evolution 1 (2): 103–13. https://doi.org/10.1111/j.2041-210X.2010.00012.x.