[R-sig-ME] different overdispersion parameter for binomial GLMM in lme4, glmmADMB and glmmTMB
Thomas Merkling
thomasmerkling00 at gmail.com
Fri Mar 16 02:51:18 CET 2018
Hi all,
I'm trying to model proportion data using a binomial GLMM. I first
started with lme4 and saw that my data were overdispersed.
LME4Inter <- glmer(PropInter ~ Manip + Observer +(1|Date),
family=binomial, weights = NbAg, data = senior,
control=glmerControl(optimizer="bobyqa",optCtrl = list(maxfun=100000)))
Using the dispersion_glmer (provided by Ben Bolker, see below for code),
I got a dispersion parameter of 1.99 (as expected the outcome is the
same with glmmADMB).
I then tried to add an observation level random effect, and the
dispersion parameter dropped to 0.24. With such a small dispersion
parameter, my model becomes overly conservative and I'll probably not be
able to detect any (true) effect (if I understand correctly). Would you
stick to this model anyways?
I started investigating beta-binomial models to see if they could better
model overdispersion (and not be too conservative) and wanted to use
glmmTMB to do that.
I first tried to refit the binomial GLMM with glmmTMB and calculate the
dispersion parameter using the dispfun function (also provided by Ben
Bolker, see below for code), and got a very different response: the
dispersion parameter was 0.514. If this is correct, I wouldn't have to
worry about overdispersion.
TMBInter <- glmmTMB(PropInter ~ Manip + Observer +(1|Date),
family=binomial, weights = NbAg, data = senior)
The estimates and SEs were identical between the 3 packages, so I am
fitting correctly the exact same models. However, the residuals between
glmer/glmmADMB and glmmTMB models were quite different. Is that expected?
I tried to run Poisson GLMMs with the 2 packages and got identical
dispersion parameters, so the dispfun function works well for glmmTMB too.
Which dispersion parameter value should I trust for the binomial GLMM?
Is there something specific about binomial GLMMs in glmmTMB?
Thanks in advance for your help!
Thomas Merkling
dispersion_glmer <- function(model) {
## number of variance parameters in
## an n-by-n variance-covariance matrix
vpars <- function(m) {
nrow(m)*(nrow(m)+1)/2
}
model.df <- sum(sapply(VarCorr(model),vpars))+length(fixef(model))
rdf <- nrow(model.frame(model))-model.df
rp <- residuals(model,type="pearson")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
dispfun <- function(m) {
r <- residuals(m,type="pearson")
n <- df.residual(m)
dsq <- sum(r^2)
c(dsq=dsq,n=n,disp=dsq/n)
}
More information about the R-sig-mixed-models
mailing list