1
$\begingroup$

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:

enter image description here

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

$\endgroup$

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Thanks a lot, Ben, for your answer (and for the suggestion on SO about cross-posting). 1) for this specific response variable (Range_filling), the median is 6% and 3rd Qu. is 14%. But I have another response variable (Range_expansion) for which the median is 66% and 3rd Qu. is 90%, and it looks not 0-inflated. Probably for this response variable it makes more sense to use a beta-binomial, without the 0-inflation component? But will it be possible to compare the results from two models if one is 0-infl and the other is not? (continue in the next comment) $\endgroup$
    – LT17
    Commented Jun 28, 2023 at 15:22
  • $\begingroup$ However, if I try it, I run into Model convergence problem; non-positive-definite Hessian matrix and false convergence (8). Removing the random effect removes the non-positive-definite Hessian matrix error, but the false convergence persists. I checked the troubleshooting vignette but I couldn't really fix the problem. (continue) $\endgroup$
    – LT17
    Commented Jun 28, 2023 at 15:24
  • $\begingroup$ 2) I tried to log-transform and scale (should I just scale?) all my predictors before fitting the model; the resultant model has a bigger overdispersion, and Residence_time_d (it was already log-transformed and scaled, and was found to be significant by the first model), after the transformation of all the predictors it's not important anymore. Are the predictors that stay significant in both models are the most "important"? (continue) $\endgroup$
    – LT17
    Commented Jun 28, 2023 at 15:25
  • $\begingroup$ 3) I checked for 0-infl by looking at the histogram of my response variable (skewed to the 0, but I know this doesn't mean 0-infl), by performing the same glmmTMB model but without ziformula (it results in an AIC of 509, so slighlty higher), and by running DHARMa::testZeroInflation(mod). Note that this last one results in a ratioObsSim = 0.6, which from ?testZeroInflation seems to not reflect 0-infl. 4) thanks for the clarification regarding ziformula = ~1, I guess it's better to stick to that then. Overall, thanks a lot for your complete answer, and sorry for the long replies. $\endgroup$
    – LT17
    Commented Jun 28, 2023 at 15:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.