# Plot a generalised mixed effects model with binomial errors

Plotting a glm binomial model is reasonably simple with the predict function. I'm having trouble creating a similar plot for a glmer model; predict doesn't work:

id    <- factor(rep(1:20, 3))
age   <- rep(sample(20:50, 20, replace=T), 3)
age   <- age + c(rep(0, 20), rep(3, 20), rep(6, 20))
score <- rbinom(60, 15, 1-age/max(age))
dfx   <- data.frame(id, age, score)

library(lme4)
glmerb  <- glmer(cbind(score, 15-score) ~ age + (1|id), dfx, family=binomial)
ndf     <- expand.grid(age=10:60) #for extensibility, usually also have factors
ndf$fit <- predict(glmerb, ndf, type="response") *Error in UseMethod("predict") : no applicable method for 'predict' applied to an object of class "mer"*  1. How can I produce the desired plot? 2. While I'm at it, what other plots would be useful for this kind of model for either diagnostic, presentation or glam purposes? ## 3 Answers The issue is that the mer class in R, and the lmer etc commands in lme4 all produce mer objects, and these are not compatible with some "normal" R commands. You can get the fits out by using fitted but amending your code to ndf$fit <- fitted(glmerb, ndf, type="response")


gave me the error

Error in \$<-.data.frame(*tmp*, "fit", value = c(0.213527879025905,  :
replacement has 60 rows, data has 51


But this worked:

> fit <- fitted(glmerb, ndf, type="response")
> str(fit)
num [1:60] 0.214 0.282 0.335 0.154 0.335 ...


Is that what you were after?

• Close. What I wanted was perhaps simpler, and that is the straight fixed effects curve. Not so much individual predictions. However this is useful and satisfies Q2. Perhaps a manual command to "curve" will do the trick? Feb 3 '12 at 5:12

Take a look at the ez package, particularly at ezPredict.

ps. if you would like to use 'to_predict' parameter, you'll need the dev version, see the instructions here: https://github.com/mike-lawrence/ez

To plot curves of fixed effects, I typically use code like this:

model.coefs <- fixef(model)
curve( invlogit( cbind(1, x) %*% model.coefs ), add=TRUE )


Note that invlogit is in the arm package.