Skip to content

Maximum nAGQ=25?

1 message · Rafael Sauter

#
Sorry for double posting.
But last Friday I forgot to change the subject line.
This is still about the nAGQ maximum.

Thanks for your replies.
I am aware that 25 is already a large number of quadrature points for a
scalar random effect.
I had one particular example in mind when asking you about this 25
quadrature points. It is the toenail data discussed for example in
section 14.8 (p. 278) of "Models for discrete longitudinal
data" (Molenberghs, Verbeke, 2005).This data is available on
http://onlinelibrary.wiley.com/journal/10.1111/%28ISSN%
291467-9876/homepage/50_3.htm or in the glmmAK package.

As illustrated with the code below, when using just a scalar random
intercept the differences in the estimates is not very large. Though it
looks different when using a random intercept and a random slope:

library(lme4)
packageVersion("lme4")
# [1] ?0.999999.2?
library(glmmAK)
data(toenail)

RI1 <- glmer(infect~trt+time+trt:time+(1|idnr), data=toenail,
family=binomial, nAGQ=25)
RI2 <- glmer(infect~trt+time+trt:time+(1|idnr), data=toenail,
family=binomial, nAGQ=50)
fixef(RI2)-fixef(RI1)

RS1 <- glmer(infect~trt+time+trt:time+(1+time|idnr), data=toenail,
family=binomial, nAGQ=25)
RS2 <- glmer(infect~trt+time+trt:time+(1+time|idnr), data=toenail,
family=binomial, nAGQ=50)
fixef(RS2)-fixef(RS1)


The fixed effect intercept differs by 0.139 when using nAGQ=25 or
nAGQ=50. For the model with a random intercept only the differences
reduce to the third digit. Comparable differences can be found when
using NLMMIXED in SAS.

I don't know if this is a large difference or not but the point is that
I assume from this example that with even more complicated random
effects the upper limit of 25 might be too low. Or is it just generally
a bad idea to fit GLMM models which still have varying estimates with
nAGQ=25?
So my point is that the maximum at 25 should perhaps be reconsidered
when it comes to implementing non-scalar random effects.

Rafael Sauter



On Fri, 2013-09-27 at 12:00 +0200,
r-sig-mixed-models-request at r-project.org wrote: