I've created a model with lme4's lmer and wanted to create a prediction interval around my model fit. I figured I could do it Bayesian-style and simulate from the model to do this. My model looks like:
b3 <- lmer (dollars ~ 1 + I(refund/days) + (1 | month) + (1 | regime), data=elect5)
and my first attempt to graph this was with:
plot (elect5$dollars, type="o", ylim=c(200, 1000), ylab="Y", xlab="X")
for (i in 1:5000) lines (simulate (b3), col=rgb (0, 0, 0, 0.01))
lines (fitted (b3), col="green")
points (elect5$dollars, col="red")
abline (h=fixef (b3)[1], lty=3, col="blue")
which results in the graph (red points actual data, green line fitted model, gray "prediction interval"):
which doesn't look so good. It looks like it's simply a variance around the fixed effect intercept of the model. OK, so then I do things the "hard way" and code up a function to implement the model:
doB3mc <- function (x, f, r) f[1] + (f[2] * (x$refund / x$days)) + (r[x$month]) + (r[as.numeric (x$regime) + 12])
And then fill a matrix with simulations based on mcmcsamp, which samples the coefficients (not done the most efficient way, I'm sure):
a <- mcmcsamp (b3, 5000, saveb=T)
res <- matrix (NA, nrow=5000, ncol=49)
for (i in 1:5000) for (j in 1:49) res[i, j] <- doB3mc (elect5[j,], a@fixef[,i], a@ranef[,i])
And also graph:
plot (elect5$dollars, type="o", ylim=c(200, 1000), ylab="Y", xlab="X")
for (i in 1:5000) lines (res[i,], col=rgb (0, 0, 0, 0.01))
lines (fitted (b3), col="green")
points (elect5$dollars, col="red")
abline (h=fixef (b3)[1], lty=3, col="blue")
Resulting in the graph:
Which looks more like what I'd expect, except it's biased: at values below the fixed effect intercept, the model is low in the prediction interval, and vice versa.
So my questions are: 1) why does the first option, with simulate work strangely, 2) why does the second option look biased (fitted line not centered in the interval), and most importantly 3) is either of these the proper method for doing what I set out to do?
