[R-sig-ME] predictions from a binomial GLMM

[Malcolm Fairbrother]

Dear list,

This issue has come up before, but after reading previous threads I'm still not quite clear on how to do this. After fitting a binomial model, with "cbind(successes, failures)" as the outcome and a single covariate "x", how would one generate predicted probabilities for different values of x? Below is an example, with my attempt to do this, where I am interested primarily in plotting a 95% confidence interval for predictions for a new/hypothetical unit (i.e., one not included in the data and thus for which no random effect has been estimated).

I thought that using a parametric bootstrap would be appropriate and reasonably straightforward, allowing me to get predictions with something like:

newdata %*% fixef(mod) .

However, from what I've seen (e.g., on Ben Bolker's wiki at http://glmm.wikidot.com/faq, and related threads such as https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q2/016356.html), something more is needed… though I'm not clear on whether that remains true:
(a) if one uses a bootstrap, and
(b) if the predictions are for new units (rather than ones for which the model has produced random effects).

The code below shows where I've got to. Any further assistance would be much appreciated.

Thanks,
Malcolm


library(lme4.0); library(multicore)
n <- 30
dat <- data.frame(obs=1:n, x=runif(n, -1, 1))
dat <- with(dat, data.frame(dat, y=rbinom(n, 10, prob=plogis(-1 + 2*x + rnorm(n)))))
(mod <- lmer(cbind(y,10-y) ~ x + (1 | obs), dat, family=binomial))
sims <- do.call("rbind", mclapply(1:200, function(x) fixef(refit(mod, simulate(mod))))) # parametric bootstrap
newdat <- cbind(1,seq(-1,1,0.1)) # define datapoints for which to calculate predictions
preds <- apply(sims, 1, function(x) newdat %*% x) # get predictions for each resample
loCI <- apply(preds, 1, function(x) quantile(x, 0.025)) # get, for example, the lower bound of a 95% CI
loCI <- plogis(loCI) # convert the lower bound to the data scale