Plot a generalised mixed effects model with binomial errors

Question

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"*

How can I produce the desired plot?
While I'm at it, what other plots would be useful for this kind of model for either diagnostic, presentation or glam purposes?

Of possible interest: Using lmer for prediction. And it's probably worth to check what Doug Bates has to say in his draft on lme4: Mixed-effects Modeling with R (I didn't check it, though), or maybe LMERConvenienceFunctions on CRAN. — chl, Jan 13 '12 at 13:24
Unfortunately not much in his lme4 drafts. Although I will check his book. This came close http://glmm.wikidot.com/faq but there's no glmer plot only lmer. — Matt Albrecht, Jan 16 '12 at 05:10

score 3 · Answer 1 · answered Jan 30 '12 at 04:21

3

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

answered Jan 30 '12 at 04:21

stop.psy

31

score 2 · Accepted Answer · answered Jan 30 '12 at 05:57

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? — Matt Albrecht, Feb 03 '12 at 05:12

score 2 · Answer 3 · edited Jul 12 '12 at 20:45

2

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.

edited Jul 12 '12 at 20:45

chl

53,725

answered May 17 '12 at 03:57

Jeremy

21

Plot a generalised mixed effects model with binomial errors

3 Answers3