# [R-sig-ME] prediction intervals from a mixed-effects models?

Gorjanc Gregor Gregor.Gorjanc at bfro.uni-lj.si
Wed Apr 16 09:20:56 CEST 2008

Spencer,

>      Thanks for the suggestion.  Scanning "help(pac=Zelig)", I didn't
> see anything that sounded to me like "prediction intervals from a
> mixed-effects model".  Without a more focused suggestions, it seems
> simpler to me to transform somehow the MCMC results into the desired
> prediction intervals.

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