[R-sig-ME] Error from glmer() that I do not understand.

[Rolf Turner]

On 21/10/19 2:54 PM, Ben Bolker wrote:
> 
>   Great question (as usual). It is possible to fit with glmer - by
> skipping the "nAGQ0" step completely, which usually helps but sometimes
> hurts - but you also have to tweak starting values (see below). Not
> immediately obvious to me why this is a particularly tough problem;
> glmmTMB and GLMMadaptive can do it ...
> 
> library(lme4)
> 
> ## 'works', but crazy answer
> fit <- glmer(cbind(Dead,Alive) ~ (Trt+0)/Dose +
>               (Dose | Rep),
>                data=Xdemo,
>               family=binomial(link="logit"),
>               control=glmerControl(nAGQ0initStep = FALSE))
> 
> 
> ## in this case starting with GLM fixed-effect estimates
> ##  seems to work OK (switched optimizers too)
> fit0 <- glm(cbind(Dead,Alive) ~ (Trt+0)/Dose,
>                data=Xdemo,
>                family=binomial(link="logit"))
> summary(coef(fit0))
> fit5 <- update(fit,
>                 start=list(fixef=coef(fit0),theta=c(1,0,1)),
>                 control=glmerControl(nAGQ0initStep = FALSE,
>                                      optimizer="nloptwrap"))
> 
> ## cross-check with glmmTMB
> 
> library(glmmTMB)
> fit2 <- glmmTMB(cbind(Dead,Alive) ~ (Trt+0)/Dose +
>               (Dose | Rep),
>                data=Xdemo,
>               family=binomial(link="logit"))
> 
> nrow(Xdemo) ## 99
> length(fixef(fit2)$cond + getME(fit2,"theta")) ## 12
> 
> library(GLMMadaptive)
> ## tried to do this with GLMMadaptive,
> ##
> https://hypatia.math.ethz.ch/pipermail/r-sig-mixed-models/2019q3/028057.html
> fit3 <- mixed_model(cbind(Dead,Alive) ~ (Trt+0)/Dose,
>                      random = ~ Dose | Rep,
>                      data=Xdemo,
>                      family=binomial(link="logit"))
> 
> fit4 <- MASS::glmmPQL(cbind(Dead,Alive) ~ (Trt+0)/Dose,
>                      random = ~ Dose | Rep,
>                      data=Xdemo,
>                      family=binomial(link="logit"))
> 
> ## comparison.  first glmer fit is terrible.
> ## glmmPQL is a little farther off from others (not surprising ...)
> 
> cbind(glmer=fixef(fit),glmmTMB=fixef(fit2)$cond,
>        glmer2=fixef(fit5),
>        GLMMadaptive=fixef(fit3),
>        glmmPQL=fixef(fit4))
> 
> Xdemo <- transform(Xdemo,
>                      n= Alive+Dead,
>                      prop=Dead/(Alive+Dead))
> 
> library(ggplot2)
> ggplot(Xdemo,aes(Dose,prop,colour=Trt)) +
>      geom_point(aes(size=n),alpha=0.5) +
>      geom_line(aes(group=Rep),alpha=0.5)

Thanks Ben for the very thorough answer.  (How on earth do you find the 
time to respond the way you do?)

The error that I described actually arose in the context of a fairly big
(by my standards) simulation study.  Such errors seemed to arise very 
frequently and the randomness of the data make it difficult to keep the 
misbehaviour under control.  I need the simulations to run reliably and
not crash in the middle of a test sequence.

I may have to set things up to catch errors and if there is one to 
switch the fitting method from glmer() to mixed_model() --- or vice 
versa.  (I'm sure that GLMMadaptive will crash unexpectedly with certain 
simulated data sets, just as lme4 does.)

Life never does go smoothly! :-)

Thanks again.

cheers,

Rolf

-- 
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276