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

[Ben Bolker]

 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)





On 2019-10-20 8:58 p.m., Rolf Turner wrote:
> 
> Xdemo <- structure(list(Dead = c(4, 26, 71, 34, 67, 68, 32, 39, 2, 15,
> 64, 26, 29, 0, 11, 10, 89, 87, 30, 121, 2, 9, 20, 0, 4, 30, 21,
> 21, 48, 54, 0, 47, 73, 87, 33, 73, 41, 39, 0, 25, 33, 19, 53,
> 39, 20, 3, 26, 26, 80, 79, 9, 35, 39, 40, 81, 80, 103, 7, 63,
> 92, 89, 40, 76, 0, 30, 29, 44, 64, 24, 72, 49, 0, 72, 73, 78,
> 46, 3, 21, 34, 29, 25, 37, 52, 64, 49, 0, 11, 11, 0, 17, 13,
> 56, 28, 3, 19, 23, 21, 26, 33), Alive = c(50, 22, 16, 6, 6, 1,
> 1, 0, 45, 5, 3, 0, 0, 55, 70, 18, 24, 4, 0, 0, 27, 11, 1, 35,
> 23, 17, 0, 0, 0, 0, 30, 12, 4, 1, 0, 0, 0, 0, 35, 29, 11, 6,
> 5, 2, 0, 61, 21, 7, 1, 0, 13, 0, 0, 0, 0, 0, 0, 26, 0, 0, 0,
> 0, 0, 26, 30, 5, 0, 0, 0, 0, 0, 74, 24, 11, 0, 0, 18, 11, 0,
> 0, 14, 6, 6, 8, 0, 40, 17, 9, 22, 38, 11, 2, 0, 51, 10, 5, 3,
> 0, 0), Rep = structure(c(6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 12L,
> 12L, 12L, 12L, 12L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 3L, 3L,
> 3L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 15L, 15L, 15L, 15L, 4L, 4L, 4L,
> 4L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 16L, 16L, 16L, 16L, 16L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 7L, 7L, 7L, 7L, 7L, 7L, 13L, 13L,
> 13L, 13L, 13L, 13L, 13L, 13L, 5L, 5L, 5L, 5L, 5L, 11L, 11L, 11L,
> 11L, 17L, 17L, 17L, 17L, 17L, 2L, 2L, 2L, 8L, 8L, 8L, 8L, 8L,
> 14L, 14L, 14L, 14L, 14L, 14L), .Label = c("1", "2", "3", "4",
> "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15",
> "16", "17", "18"), class = "factor"), Dose = c(0, 15, 30, 45,
> 55, 65, 75, 85, 0, 30, 40, 60, 80, 0, 10, 20, 30, 40, 60, 80,
> 0, 10, 35, 0, 10, 20, 35, 40, 45, 50, 0, 25, 30, 35, 30, 45,
> 50, 55, 0, 15, 20, 25, 30, 35, 40, 0, 10, 15, 25, 35, 0, 10,
> 15, 18, 21, 24, 27, 0, 10, 18, 21, 24, 37, 0, 8, 12, 15, 18,
> 22, 20, 24, 0, 10, 14, 26, 30, 0, 4, 8, 13, 4, 8, 10, 12, 18,
> 0, 6, 9, 0, 3, 6, 13, 18, 0, 6, 8, 10, 12, 14), Trt = structure(c(6L,
> 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
> 6L, 6L, 6L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
> 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
> 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
> 5L, 5L, 5L, 5L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
> 2L, 2L), .Label = c("16hour10deg", "16hour20deg", "16hour5deg",
> "8hour10deg", "8hour20deg", "8hour5deg"), class = "factor")), class =
> "data.frame", row.names = c(1L,
> 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 11L, 12L, 14L, 15L, 16L, 17L,
> 18L, 19L, 20L, 22L, 23L, 24L, 25L, 26L, 28L, 29L, 30L, 31L, 32L,
> 33L, 34L, 35L, 38L, 39L, 40L, 43L, 45L, 46L, 47L, 48L, 49L, 50L,
> 51L, 52L, 53L, 54L, 55L, 56L, 57L, 59L, 61L, 63L, 64L, 65L, 66L,
> 67L, 68L, 69L, 71L, 72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L,
> 81L, 82L, 83L, 84L, 85L, 87L, 88L, 90L, 91L, 92L, 93L, 94L, 95L,
> 98L, 99L, 100L, 101L, 104L, 105L, 106L, 107L, 111L, 112L, 113L,
> 115L, 116L, 117L, 120L, 121L, 122L, 123L, 124L))
> 
> library(lme4)
> fit <- glmer(cbind(Dead,Alive) ~ (Trt+0)/Dose +
>              (Dose | Rep),
>               data=Xdemo,
>               family=binomial(link="logit"),nAGQ=0)
> 
> This throws the error:
> 
>> Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GHrule(0L),
>> compDev = compDev,  :   (maxstephalfit) PIRLS step-halvings failed to
>> reduce deviance in pwrssUpdate
> 
> Setting nAGQ=1 makes no difference.