CODEWITHSUNDEEP
lme4
Ewen
Ewen

Reputation: 1381

Simulation from a nested glmer model

I'm having a problem generating simulations from a 3 level glmer model when conditioning on the random effects (I'm actually using predict via bootMer but the problem is the same).

This works:

library(lme4)
fit1 =   glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                                            data = cbpp, family = binomial)
simulate(fit1, re.form=NULL)

This fails:

cbpp$bigherd = rep(1:7, 8)
fit2 =   glmer(cbind(incidence, size - incidence) ~ period + (1 | bigherd / herd),
                             data = cbpp, family = binomial)
simulate(fit2, re.form=NULL)    
Error: No random effects terms specified in formula

Many thanks for any ideas.

Update

Ben, many thanks for your help below, really appreciate it. I wonder if I can impose on you again.

What I want to do is simulate predictions on the response scale and I'm not sure if I can use your work around? Or if there is an alternative to what I'm doing. Thank you!

This works as expected, but is not conditional on random effects:

FUN = function(.){
    predict(., type="response")
}
bootMer(fit2, FUN, nsim=3)$t

This doesn't work, as would be expected given above problem:

bootMer(fit2, FUN, nsim=3, use.u=TRUE)$t

As far as I can see, I can't pass re.form to bootMer.

Does the alternative below result in simulated predictions conditional on random effects without passing use.u to bootMer?

FUN = function(.){
    predict(., type="response", re.form=~(1|herd:bigherd) + (1|bigherd))
}
bootMer(fit2, FUN, nsim=10)$t

Upvotes: 2

Views: 648

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226342

I'm not sure what's going on yet, but here are two workarounds that do work:

simulate(fit2, re.form=lme4:::reOnly(formula(fit2)))
simulate(fit2, re.form=~(1|herd:bigherd) + (1|bigherd))

There must be something going wrong with the expansion of the "slash" term, because this doesn't work:

simulate(fit2, re.form=~(1|bigherd/herd))

I've posted this as an lme4 issue

These workarounds don't work for bootMer (which only takes the use.u argument, not re.form) in the current CRAN release (1.1-9).

It is fixed in the development version on Github (1.1-10): devtools::install_github("lme4/lme4") will install it, if you have compilation tools installed.

In the meantime you could just go ahead and implement your own parametric bootstrap (for parametric bootstrapping, bootMer is actually a very thin wrapper around simulate()/[refit()orupdate()]/FUN`). Much of the complication has to do with parallel computation (you'd have to add some of it back in if you want parallel computation in your own PB implementation).

This is the outline of a hand-rolled parametric bootstrap:

nboot <- 10
nresp <- length(FUN(orig_fit))
res <- matrix(NA,nboot,nresp)
for (i in 1:nboot) {
  res[i,] <- FUN(update(orig_fit,data=simulate(orig_fit,...)))
             ## or use refit() for LMMs
             ## ... are options to simulate()
}
t(apply(res,2,quantile,c(0.025,0.975)))

Upvotes: 2

Related Questions