lme4 / lme4

Mixed-effects models in R using S4 classes and methods with RcppEigen
Other
619 stars 145 forks source link

another edge case: cloglog/binomial -> pwrssUpdate failure #68

Closed bbolker closed 10 years ago

bbolker commented 11 years ago
## testing identity, cloglog links, gaussian family
set.seed(1002)
nblock <- 10
nperblock <- 50
sd.u <- 1
beta <- c(0.2,0.5)
ntot <- nblock*nperblock
d <- data.frame(x=runif(ntot),f=factor(rep(LETTERS[1:nblock],each=nperblock)))
r <- rnorm(nblock,mean=0,sd=sd.u)
gshape <- 1.5
d$offset <- rgamma(ntot,1,1)
d <- within(d,
            {
              eta0 <- beta[1]+beta[2]*x+offset
              eta <- eta0+r[f]
            })

## cloglog:
cc <- binomial(link="cloglog")
d$mu0 <- cc$linkinv(d$eta0)
d$mu <- cc$linkinv(d$eta)

d$y0 <- rbinom(ntot,prob=d$mu0,size=1)
d$y <- rbinom(ntot,prob=d$mu,size=1)

## 

library(lme4)
g1A <- glmer(y~x+(1|f)+offset(offset),data=d,
             family=binomial(link="cloglog"),
             optimizer="bobyqa",
             control=glmerControl(tolPwrss=1e-3))

if (require(glmmADMB)) {
g1 <- glmmadmb(y~x+(1|f)+offset(offset),data=d,
               family="binomial",link="cloglog")

dput(coef(g1))
dput(sqrt(VarCorr(g1)[[1]]))
}

glmmadmb_ests <- list(theta=1.0866922287382,
                             fixef=c(0.964558165725369, 0.039712606911609))

g1B <- glmer(y~x+(1|f)+offset(offset),data=d,
             family=binomial(link="cloglog"),
             start=glmmadmb_ests,
             control=glmerControl(tolPwrss=1e-3),
             verbose=20)

(Still gets nan: something still isn't clamped!)

mmaechler commented 11 years ago

maybe I need to be told that I need to install non-CRAN packages.... but I don't have glmmADMB on my laptop and the code above uses coef(g1) Can you just dput(coef(g1)) and show the output of that... so I could add the above to an example collection ... I've been musing about.

bbolker commented 11 years ago

See edited version on Github (since you're probably reading in e-mail, the fixed effect coefficients are: c(0.964558165725369, 0.039712606911609)