Solved – interpreting output for glmmTMB for zero-inflated count data

generalized linear modelglmmtmblme4-nlmemultilevel-analysisr

I have been trying to read all the documentation I have, but I'm still not sure what the difference is between the "conditional" and zero-inflated models in the output of the glmmTMB. Below is some code modelling a zero-inflated count dependent variable (alcoholic drinks, standardized). The observations are nested within persons (153 persons, 90 observations each).

I am wondering if the "Conditional model" is simply the original model that doesn't consider zero-inflation, whereas the "zero-inflated" is the model with this factored in (hence I would report these stats only?). In addition, I am wondering whether the zero-inflated model includes 0 values or not, given that the estimates change in direction (and sometimes in significance). Or maybe I have this all the other way around. Any clarifications would be helpful.

Also when playing around with the output, the allEffects function only plots the conditional model. I am wondering if there's an easy way to plot the zero-inflated model as well.

Thanks in advance!

> model1 <- glmmTMB(drinks_round ~ sex + cann_used + cann_g + other_type + n_used + day + wDay 
                                   + (1|studyID), 
                   zi = ~sex + can_used + cann_g + other_type + n_used + day + wDay,
                   family=nbinom2, data=can1
                   )

> summary(model1)
 Family: nbinom2  ( log )

Formula:
   drinks_round ~ sex + cann_used + cann_g + other_type + n_used + day + wDay 
                  + (1 | studyID)
Zero inflation:
   ~sex + cann_used + cann_g + other_type + n_used + day + wDay
Data: can1

     AIC      BIC   logLik deviance df.resid 
 24795.8  24931.4 -12379.9  24759.8    13752 

Random effects:

Conditional model:
 Groups  Name        Variance Std.Dev.
 studyID (Intercept) 0.4657   0.6825  
Number of obs: 13770, groups:  studyID, 153

Overdispersion parameter for nbinom2 family (): 4.97 

Conditional model:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      1.1459925  0.2080129   5.509 3.60e-08 ***
sex             -0.1729803  0.1209996  -1.430   0.1528    
cann_used        0.0640049  0.0416621   1.536   0.1245    
cann_g           0.1866736  0.0311277   5.997 2.01e-09 ***
other_type      -0.0330603  0.0691555  -0.478   0.6326    
n_used           0.3266708  0.0660834   4.943 7.68e-07 ***
day              0.0010574  0.0005845   1.809   0.0705 .  
wDayweekday     -0.2464173  0.0314343  -7.839 4.54e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.8868254  0.1014340   8.743  < 2e-16 ***
sex              0.0638217  0.0495552   1.288   0.1978    
cann_used       -0.2525820  0.0593977  -4.252 2.12e-05 ***
cann_g           0.2174713  0.0410397   5.299 1.16e-07 ***
other_type      -0.1837736  0.0996293  -1.845   0.0651 .  
n_used          -0.3920739  0.0720337  -5.443 5.24e-08 ***
day              0.0002897  0.0008903   0.325   0.7449    
wDayweekday      0.3367418  0.0487646   6.905 5.00e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Best Answer

"The zero-inflation model estimates the probability of an extra zero such that a positive contrast indicates a higher chance of absence [...]; this is the opposite of the conditional model where a positive contrast indicates a higher abundance"

So, answering your comment, it is indeed the likelihood of zero (the alternative being "not 0" rather than 1). By the way the ziformula also allows you to add random effects, just like in your conditional formula.

Pages 382-383 explain all components of the model summary:

Brooks, M. E., Kristensen, K., van Benthem, K. J., Magnusson, A., Berg, C. W., Nielsen, A., Skaug, H. J., Machler, M., & Bolker, B. M. (2017). glmmTMB balances speed and flexibility among packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal, 9(2), 378-400. https://doi.org/10.32614/RJ-2017-066

Related Solutions

Solved – Can you use glmmTMB to simultaneously model offsets and zero-inflation

tl;dr as far as I can tell at this point,

 glmmTMB(formula=<...>+offset(log(sampling_effort)),
         ziformula = ~.,
         family=nbinom2,
         data=<...>)

should do what you want. (1) ziformula specifies zero-inflation. (2) The offset term in the conditional model (formula) adjusts for sampling effort; as I understand your context, you shouldn't need any adjustment for effort/sample depth in the zero-inflation part of the model (since it describes structural zeros, which will always be observed as zero regardless of sampling effort). (3) family=nbinom2 takes care of other sources of overdispersion. [You might want to alternately consider nbinom1, which specifies $\textrm{Var} \propto \mu$ rather than $\textrm{Var} = \mu + \mu^2/k$.]

This is an interesting question both from the statistical and the implementation point of view.

Implementation: you can add offsets to zero-inflation terms, you just can't do it via .. For example something like

 glmmTMB(y~x,
         family=nbinom2,
         zi=~x+offset(log(w)))

should work fine. It's only if you try to use zi = ~ . to match the conditional formula in a lazy way that the offset gets dropped.

Statistics: I question a couple of your premises.

First of all, it's not immediately obvious that different numbers of counts should lead to biased results; it's important to know what the source of the differences is - i.e. large natural variations in density, or variations in searching/capture effort?
Second, you have to think carefully about the form of the offset. Using log(effort) as the offset with a log-link zero-inflation term (the only choice at present) will mean that the probability of an observation being a structural zero will be proportional to the effort, i.e. $\log(p_z) = ... + \log(e) \to p_z \propto e$. In general a complementary log-log link (with log(effort) as the offset) is more appropriate for detection probabilities, as that makes the hazard of finding something proportional to effort. However, if you're really trying to model structural zeros I question whether search effort should influence this part of the model at all ...

Based on the comments, I think this question may be based on a (reasonable) misunderstanding of ?glmmTMB, which currently reads

ziformula: ... Offset terms will automatically be dropped from the conditional effects formula when using ‘~.’

This warning applies only to the zero-inflation formula: the conditional formula (formula) argument isn't modified at all, it's only the zero-inflation version of the formula (and, as discussed above, you can include an offset in the zero-inflation part of the model if you really want to, by writing out the formula explicitly). Thus if you follow the fairly standard procedure of adding + offset(log(sampling_effort)) to formula then the conditional mean number of counts will be proportional to effort (assuming a [default] log-link model for the counts).

In hopes of clarifying this I've tried editing this statement: is the following clearer?

Specifying ~. will set the right-hand side of the zero-inflation formula identical to the right-hand side of the main (conditional effects) formula; terms can also be added or subtracted. Offset terms will automatically be dropped from the conditional effects formula when using ~.

(Hmmm, now that I read this over it doesn't seem much better ...) Feel free to comment (or suggest clearer wording) here, or at https://github.com/glmmTMB/glmmTMB/issues ...

Solved – GLMM hurdle model for continuous data -Truncated negative binomial family in glmmTMB

I'm not sure why you say that glmmTMB can't handle zero-inflated Gamma responses: the glmmTMB news file says (for version 1.0.0, release 2020-02-03):

new ziGamma family (minor modification of stats::Gamma) allows zero-inflation (i.e., Gamma-hurdle models)

I'd say it's not crazy to use a truncated negative binomial, but I'd be worried as it doesn't make statistical sense (technically, the likelihood of any non-integer value is 0 ...) (If you really had count data, a zero-inflated NB rather than a hurdle would be a reasonable option ...)

Given the distribution/density functions for each distribution, parameterized in terms of the mean $\mu$ and a shape/dispersion parameter ($k$ for NB, $a$ for Gamma):

$$ \textrm{NB}: \qquad \frac{(k/(k+\mu))^k}{\Gamma(k)} \cdot \frac{\Gamma(k+x)}{x!} \cdot (\mu/(k+\mu))^x $$

$$ \textrm{Gamma}: \qquad \frac{1}{(\mu/a)^a \Gamma(a)} \cdot x^{a-1} \cdot e^{-(x/(\mu/a))} $$

I think you may? be able to show that NB converges approximately to Gamma for large $x$ (but someone better/more dedicated than I am will need to do the math ...)

An empirical demonstration (not "proof"!)

hist(rnbinom(100000,mu=100,size=2),freq=FALSE,ylim=c(0,0.008), breaks=100)
curve(dgamma(x,scale=100/2,shape=2),add=TRUE,col=2,lwd=2)

This definitely doesn't work for small mean (try it with mean=4 rather than 100 ...)

Best Answer

Related Solutions

Solved – Can you use glmmTMB to simultaneously model offsets and zero-inflation

Solved – GLMM hurdle model for continuous data -Truncated negative binomial family in glmmTMB

Related Question