Prediction on mixed effect models: what to do with random effects?

Let’s consider this hypothetical dataset:

set.seed(12345)

num.subjects <- 10

dose <- rep(c(1,10,50,100), num.subjects)
subject <- rep(1:num.subjects, each=4)
group <- rep(1:2, each=num.subjects/2*4)

response <- dose*dose/10 * group + rnorm(length(dose), 50, 30)

df <- data.frame(dose=dose, response=response, 
                 subject=subject, group=group)

we can use lme to model the response with a random effect model:

require(nlme)
model <- lme(response ~ dose + group + dose*group, 
             random = ~1|subject, df)

I would like to use predict on the result of this model to get, for instance, the response of a generic subject of group 1 to a dose of 10:

pred <- predict(model, newdata=list(dose=10, group=1))

However, with this code I get the following error:

Error in predict.lme(model, newdata = list(dose = 10, group = 1)) : 
cannot evaluate groups for desired levels on 'newdata'

To get rid of it I need to do, for instance

pred <- predict(model, newdata=list(dose=10, group=1, subject=5))

This, however, does not really make much sense to me… the subject is a nuisance factor in my model, so what sense does it have to include it in predict? If I put a subject number not present in the dataset, predict returns NA.

Is this the wanted behaviour for predict in this situation? Am I missing something really obvious?

Answer

If you look at the help for predict.lme you will see that it has a level argument that determines which level to make the predictions at. The default is the highest or innermost which means that if you don’t specify the level then it is trying to predict at the subject level. If you specify level=0 as part of your first predict call (without subject) then it will give the prediction at the population level and not need a subject number.

Attribution
Source : Link , Question Author : nico , Answer Author : Greg Snow

Leave a Comment