Selecting between a zero-inflated binomial, OLRE and beta-binomial model

Question

I need some help in deciding which of the following models fits best the data that I have. This was a survey where participants reported proportions of successes (defined as n/m) in condition A and B. The model predicts the proportions by the binary condition variable, and continuous x and z variables (ranging from 1 to 7), as well as random effects for each subject and 13 types of task. This is the distribution of the proportions

So the model is defined as

mod_b0 <- glmmTMB(n/m ~ x*condition + z*condition + (1|subject) + (1|task), weights = m, family = binomial)
summary(mod_b0)

     AIC      BIC   logLik deviance df.resid 
 22830.4  22883.7 -11407.2  22814.4     5781 

Random effects:

Conditional model:
 Groups  Name        Variance Std.Dev.
 task    (Intercept) 0.2094   0.4576  
 subject (Intercept) 1.5546   1.2468  
Number of obs: 5789, groups:  task, 13; subject, 225

Conditional model:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.44713    0.25706 -13.410  < 2e-16 ***
x             0.38560    0.03690  10.449  < 2e-16 ***
conditionB   -1.36826    0.20133  -6.796 1.08e-11 ***
z            -0.07328    0.02276  -3.220  0.00128 ** 
x:conditionB  0.17682    0.03807   4.644 3.41e-06 ***
conditionB:z  0.12544    0.02512   4.994 5.91e-07 ***

The residuals test by DHARMa (N = 1000 simulations) suggest that there is no overdispersion, that there is zero-inflation and that the model does not fit the data well.

I tried three solutions:

1. Binomial model with zero-inflation
2. OLRE binomial model
3. Beta-binomial model

Here are the outputs of all three of them.

Binomial model with zero-inflation

mod_bzi <- glmmTMB(n/m ~ x*condition + z*condition + (1|task) + (1|subject), 
                  data = dx, family = binomial, weights = m, ziformula = ~ 1 + condition*z)
summary(mod_bzi)
    AIC      BIC   logLik deviance df.resid 
 17949.0  18029.0  -8962.5  17925.0     5777 

Random effects:

Conditional model:
 Groups  Name        Variance Std.Dev.
 task    (Intercept) 0.09208  0.3034  
 subject (Intercept) 1.95087  1.3967  
Number of obs: 5789, groups:  task, 13; subject, 225

Conditional model:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.65838    0.29974  -8.869  < 2e-16 ***
x             0.40498    0.04874   8.309  < 2e-16 ***
conditionB   -1.31011    0.26986  -4.855 1.21e-06 ***
z            -0.01559    0.02852  -0.547   0.5847    
x:conditionB  0.14559    0.05150   2.827   0.0047 ** 
conditionB:z  0.19289    0.03291   5.861 4.59e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.393898   0.084827  -4.644 3.42e-06 ***
conditionB    0.307062   0.126750   2.423   0.0154 *  
z             0.034095   0.034146   0.999   0.3180    
conditionB:z -0.003092   0.046014  -0.067   0.9464    

Please note that the regression lines in the right plot are not significantly different from the quantile lines if the number of simulations is 250!

Now we see some slight underdispersion.

OLRE model

mod_OLRE <- glmmTMB(n/m ~ x*condition + z*condition + (1|task) + (1|subject) + (1|obs_id), 
                   data = dx, family = binomial, weights = m)

     AIC      BIC   logLik deviance df.resid 
 15588.2  15648.1  -7785.1  15570.2     5780 

Random effects:

Conditional model:
 Groups  Name        Variance Std.Dev.
 task    (Intercept) 0.4361   0.6604  
 subject (Intercept) 3.0721   1.7527  
 obs_id  (Intercept) 4.8962   2.2127  
Number of obs: 5789, groups:  task, 13; subject, 225; obs_id, 5789

Conditional model:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -4.46870    0.55951  -7.987 1.38e-15 ***
x             0.43727    0.09152   4.778 1.77e-06 ***
conditionB   -2.65037    0.53953  -4.912 9.00e-07 ***
z            -0.17483    0.06014  -2.907 0.003650 ** 
x:conditionB  0.35813    0.10186   3.516 0.000438 ***
conditionB:z  0.21831    0.06827   3.198 0.001384 ** 

Again, no zero-inflation any more, but there is some under-dispersion.

Beta-binomial model

mod_bb <- glmmTMB(n/m ~ x*condition + z*condition + (1|task) + (1|subject), 
                    data = dx, family = betabinomial(link = "logit"), weights = m)

     AIC      BIC   logLik deviance df.resid 
 15305.4  15365.4  -7643.7  15287.4     5780 

Random effects:

Conditional model:
 Groups  Name        Variance Std.Dev.
 task    (Intercept) 0.2267   0.4761  
 subject (Intercept) 0.9929   0.9965  
Number of obs: 5789, groups:  task, 13; subject, 225

Overdispersion parameter for betabinomial family (): 1.54 

Conditional model:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.51074    0.33909  -7.404 1.32e-13 ***
x             0.24238    0.05426   4.467 7.94e-06 ***
conditionB   -1.31799    0.32146  -4.100 4.13e-05 ***
z            -0.08722    0.03508  -2.486  0.01291 *  
x:conditionB  0.17975    0.06081   2.956  0.00312 ** 
conditionB:z  0.09051    0.04010   2.257  0.02400 *  

Here, there is more under-dispersion in the previous models.

My conclusions and questions

• By the looks of residual distribution, it seems to me that beta-binomial model does the best job to account the data. All models have some problems with higher levels of predictors, as there are fewer cases for these values. Thus it does not wonder that the fits are somewhat poorer in that segment of the plot.
• AIC values are lowest for the beta-binomial model. However, I am not sure whether I can compare the AIC for models with different distributions of the criterion. If yes, then that would be another argument to pick the beta-binomial model.
• The coefficients are somewhat similar in beta-binomial and binomial zero-inflated models. OLRE model has some quite different coefficients. According to Harrison (2014), beta-binomial models tend to produce more reliable estimates than OLRE. Thus, I would stick with that one.

1. Would you agree with my conclusions that the beta-binomial model is the best of all proposed?
2. Is there any other way how I could improve the fit of the models that I did not think of?
3. Can I try to tweak the zero-inflation parameter in the beta-binomial model in order to get a better fit, although no zero-inflation was diagnosed by DHARMa?
4. Is there any other way how I should test the fit of the models?
5. Is underdispersion "problematic" for beta-binomial model? According to GLMM FAQ, dispersion is a problem only for models with fixed variance like binomial or poisson ones.

Answer 1

Would you agree with my conclusions that the beta-binomial model is the best of all proposed?

Yes you seem to have done a thorough job on this analysis. Your point about whether it is OK to compare these models with AIC is a good one. I remember reading conflicting information on this point, but I quickly found a reference that supports the idea that it is OK:

Hardin, J.W. and Hilbe, J.M., 2014. Estimation and testing of binomial and beta-binomial regression models with and without zero inflation. The Stata Journal, 14(2), pp.292-303. https://journals.sagepub.com/doi/pdf/10.1177/1536867X1401400204

Is there any other way how I could improve the fit of the models that I did not think of?

You could look at predictive accuracy using a train / validate / test approach.

Can I try to tweak the zero-inflation parameter in the beta-binomial model in order to get a better fit, although no zero-inflation was diagnosed by DHARMa?

It would be worth a try but given the DHARMa output probably won't improve things.

Is there any other way how I should test the fit of the models?

Again, I would suggest looking at predictions.

Is underdispersion "problematic" for beta-binomial model? According to GLMM FAQ, dispersion is a problem only for models with fixed variance like binomial or poisson ones.

Under- and over-dispersion is "handled" by beta-binomial models, so it should not be a problem.