[R-sig-ME] Error from glmer() that I do not understand.
Phillip Alday
ph||||p@@|d@y @end|ng |rom mp|@n|
Mon Oct 21 04:20:20 CEST 2019
Without taking the time to actually look in more detail about what's
going on here:
The initial nAGQ=0 / only optimize over theta and not over beta step can
do weird things for the nAGQ > 0 / optimizer over theta and beta step
when the initial fit is on the boundary. Is that what's happening here?
On 21/10/2019 12:24, 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)
>
>
>
>
>
> 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.
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list