7
$\begingroup$

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?
$\endgroup$
2
$\begingroup$

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?

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Matt Albrecht Feb 3 '12 at 5:12
3
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$
2
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.