[R-sig-ME] Quadrature appears to incorrectly calculate deviance in binomial GLMM

[Russell Millar]

Dear All,

I'm fitting to some historical grouped binomial data. There are 16 
binomial obs arranged in 8 blocks of 2 (trmt and control). I know what 
the model deviance calculated via the Laplace approx and quadrature 
should be since I've fitted the model using SAS PROC GLIMMIX, ADMB and 
via simulated likelihood (Millar 2004, ANZJS). However, lmer is only 
giving the correct deviance under the Laplace approximation.  Also, the 
deviance calculated by lmer using quadrature depends on the choice of 
quadrature points???

By the way, I'm not directly interested in the deviance. The calculated 
value of 35.8 on 13 dof would suggest lack of fit, which is nonsense. 
I've done a parametric bootstrap and this observed deviance is typical. 
I'm really interested in the accuracy of the calculated log-likelihood.

Thanks for any help,

Russell Millar
Dept of Stats
U. Akld


#Binomial data in 8 blocks of 2
clinic=rep(1:8,rep(2,8))
trmt=as.factor(rep(c("Drug","Control"),8))
y=c(11,10,16,22,14,7,2,1,6,0,1,0,1,1,4,6) #Successes
n=c(36,37,20,32,19,19,16,17,17,12,11,10,5,9,6,7) #Number of trials

#log-likelihood of saturated model
sum(log(dbinom(y,n,y/n)))
#[1] -19.19493

#Using 3 forms of alternative software, it has been calculated that
#the true log-likelihood (with consts) of this model is -37.0313
#or -37.1006 if calculated using the Laplace approximation.
#These correspond to deviances of
#2*(-19.19493+37.0313)=35.67274 (exact)
#2*(-19.19493+37.1006)=35.81134 (from Laplace approx)


library(lme4)
#Laplace approximation
deviance( lmer(y/n~trmt+(1|clinic),family="binomial",weight=n) )
#[1] 35.81135 #Looks good

#Quadrature with 2 quadrature points
deviance( lmer(y/n~trmt+(1|clinic),family="binomial",weight=n,nAGQ=2) )
#[1] 38.90244 #???

#Quadrature with 3 quadrature points
deviance( lmer(y/n~trmt+(1|clinic),family="binomial",weight=n,nAGQ=3) )
#[1] 38.90244 #At least this is the same as above

#Quadrature with 9 quadrature points
deviance( lmer(y/n~trmt+(1|clinic),family="binomial",weight=n,nAGQ=9) )
#[1] 39.05139 #Losing numerical accuracy perhaps?

#Quadrature with lots of quadrature points
deviance( lmer(y/n~trmt+(1|clinic),family="binomial",weight=n,nAGQ=99) )
#[1] 31.59593 #Yikes!