Skip to content
Prev 853 / 20628 Next

prediction intervals from a mixed-effects models?

Spencer,
library(lme4)
library(MEMSS)
set.seed(3)
data(Orthodont)
fm3r <- lmer(distance ~ age * Sex + (1 | Subject), Orthodont)
fm3r

## Call Zelig wrapper function
fmZ <- zelig(distance ~ age * Sex + tag(1 | Subject), data=Orthodont, model="ls.mixed")
fmZ

xF <- setx(fmZ, Sex = "Female")
xM <- setx(fmZ, Sex = "Male")

## First differences for sex
s.out <- sim(fmZ, x=xF, x1=xM)
summary(s.out)
plot(s.out)

## Only for females
s.outF <- sim(fmZ, x=xF)
summary(s.outF)
plot(s.outF)

## Only for Males
s.outM <- sim(fmZ, x=xM)
summary(s.outM)
plot(s.outM)

The document at

http://gking.harvard.edu/zelig/docs/ls.mixed.pdf

says the following about prediction intervals: "The predicted values (qi$pr) are draws from the normal distribution
defined by mean \mu_{ij} and variance \sigma^2

\mu_{ij} = X_{ij}\beta + Z_{ij}b_{i}

given X_{ij} and Z_{ij} and simulations of \beta and b_{i} from their posterior distributions. The
estimated variance covariance matrices are taken as correct and are themselves not simulated."

I would like to warn that functions in Zelig are a bit error prone as it pulls various bits from "lmer" objects. As we
all witness, lme4 is a fast moving target and Zelig may lag behind! I am also not sure how Zelig handles models
in lmer() that have special structure of random effects ...

Of course, you can always sample from posterior distributions of parameters obtained with mcmcsamp()
function. Perhaps rv package by Jouni Kerman can be used. Anyway, predictions are very easily
done in BUGS, say

model
{
  ## --- Likelihood ---

  for(i in 1:N) {

    y[i] ~ dnorm(mu[i], tauE)
    mu[i] <- int + ...

    y.pred[i] ~ dnorm(mu[i], tauE)
  }

  ## --- Priors ---
  ...


}

Regards, Gregor