prediction intervals from a mixed-effects models?
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