[R-sig-ME] different overdispersion parameter for binomial GLMM in lme4, glmmADMB and glmmTMB

[Thomas Merkling]

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)
}