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:
- Binomial model with zero-inflation
- OLRE binomial model
- 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.
- Would you agree with my conclusions that the beta-binomial model is the best of all proposed?
- Is there any other way how I could improve the fit of the models that I did not think of?
- 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?
- Is there any other way how I should test the fit of the models?
- 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.