Skip to content
Prev 3486 / 20628 Next

cloglog and lmer

Hmmm.  Can you remind me what didn't work?
  Examples below give at least vaguely plausible answers ...

  cheers
    Ben Bolker

set.seed(1001)

N <- 20
n <- 100
x <- runif(N*n)
f <- factor(rep(LETTERS[1:N],each=n))
dat <- data.frame(x,f)

beta <- c(1,3)
alpha <- rnorm(N,sd=0.5)

mm <- cbind(model.matrix(~x,data=dat),
            model.matrix(~f-1,data=dat))

eta <- mm %*% c(beta,alpha)
y <- rbinom(N*n,prob=1-exp(-exp(eta)),size=1)
dat <- data.frame(dat,y)

library(lme4)
(g1 <- glmer(y~x+(1|f),data=dat,
      family=binomial(link="cloglog")))

fixef(g1)  ## matches (1,3) reasonably well
VarCorr(g1) ## matches sd=0.5 reasonably well

## now try crossed random effects
set.seed(1001)  ## reset
N1 <- 10
N2 <- 10
n <-  20
ntot <- n*N1*N2
dat <- expand.grid(f1=LETTERS[1:N1],f2=letters[1:N2],rep=1:n)
dat <- data.frame(dat,x=runif(ntot))
alpha1 <- rnorm(N1,sd=0.5)
alpha2 <- rnorm(N2,sd=1)

mm <- cbind(model.matrix(~x,data=dat),
            model.matrix(~f1-1,data=dat),
            model.matrix(~f2-1,data=dat))

eta <- mm %*% c(beta,alpha1,alpha2)
y <- rbinom(ntot,prob=1-exp(-exp(eta)),size=1)

(g2 <- glmer(y~x+(1|f1)+(1|f2),data=dat,
      family=binomial(link="cloglog")))

## warning about false convergence
fixef(g2) ## still OK.
VarCorr(g2)  ## reasonable but ??
Sasha Goodman wrote: