[R-sig-ME] R-sig-mixed-models Digest, Vol 81, Issue 51

[Rafael Sauter]

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:
> Send R-sig-mixed-models mailing list submissions to
> 	r-sig-mixed-models at r-project.org
> 
> To subscribe or unsubscribe via the World Wide Web, visit
> 	https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> or, via email, send a message with subject or body 'help' to
> 	r-sig-mixed-models-request at r-project.org
> 
> You can reach the person managing the list at
> 	r-sig-mixed-models-owner at r-project.org
> 
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of R-sig-mixed-models digest..."
> 
> 
> Today's Topics:
> 
>    1. Re: Maximum nAGQ=25? (Ben Bolker)
> 
> 
> ----------------------------------------------------------------------
> 
> Message: 1
> Date: Thu, 26 Sep 2013 21:21:33 +0000 (UTC)
> From: Ben Bolker <bbolker at gmail.com>
> To: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Maximum nAGQ=25?
> Message-ID: <loom.20130926T231632-517 at post.gmane.org>
> Content-Type: text/plain; charset=us-ascii
> 
> Ross Boylan <ross at ...> writes:
> 
> > On Thu, Sep 26, 2013 at 04:23:47PM +0000, Ben Bolker wrote:
> > > Rafael Sauter <rafael.sauter <at> ...> writes:
> 
>  [snip]
>  
> > > > As I did not find any discussion about this change in the new
> > > > lme4-version let me allow to ask:
> > > > 
> > > > 1) Why is 25 a reasonable upper bound for nAGQ? What were the reasons to
> > > > implement this upper bound? Is the increasing complexity as mentioned in
> > > > the details of '?glmer' the the main reason for this?
> > > > 
> 
>  [snip]
> 
> > >    I will only speak for myself: other lme4-authors (especially Doug
> > > Bates) may chime in on this one.  I believe there isn't a rigorous
> > > argument for why >25 quadrature points is too many: ?glmer says
> > > " A model with a single, scalar random-effects term could
> > > reasonably use up to 25 quadrature points per scalar integral."
> 
> [snip]
> 
> > I think we would certainly
> > > be willing to reconsider this limit if you can show that there is
> > > some sensible case where it matters ...
> > 
> > If the limit is hard-coded to 25, it will be hard to discover if using
> > >25 matters.  That seems to me an argument for not hard coding it.  I
> > suppose if the results had not stabilized by 25 that would be an
> > indication.
> > 
> > OTOH, 25 is a lot of quadrature points.
> > 
> > One problem I've encountered with high number of quadrature
> > points--not in lmer, but I think it's a general issue--is that as the
> > number of quadrature points goes up the extreme x values go up, and
> > numerical problems are more likely.  Usually one can compensate by
> > coding the likelihood defensively.
> > 
> > Ross Boylan
> 
>   As may have been reflected in discussions on the list, the lme4
> authors have been having lots of internal discussions about how much
> flexibility to allow, when to give users helpful advice in the form
> of warnings, etc etc etc.. Not completely tongue-in-cheek, I could say
> that if you're capable of compiling a package from source, it's not 
> very complicated to search for "nAGQ <= 25L" in R/modular.R and modify
> or remove this limitation for yourself ... 
> 
>   Ben Bolker
> 
> 
> 
> ------------------------------
> 
> _______________________________________________
> R-sig-mixed-models mailing list
> R-sig-mixed-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 
> End of R-sig-mixed-models Digest, Vol 81, Issue 51
> **************************************************