# Finding a prediction interval for an lmer model via mcmcsamp versus simulate

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?

-
Have you looked at the residuals to see if the ones associated with values below the intercept tend to be positive and those above the intercept tend to be negative? –  Michelle Feb 3 '12 at 10:05
@Michelle: Good call. I'll look at it this evening when I get home. Doing an eyeball of the actual (red) versus fitted (green), it doesn't look like this is the case, but it's easy to think "it looks like it evens out" when it doesn't. –  Wayne Feb 3 '12 at 16:01
This is a difficult problem, the prediction intervals don't look like they centre on either the fitted values or the actual data. With 5000 predictions, you would think the intervals would be quite nice. –  Michelle Feb 3 '12 at 18:30
@Michelle: I think the issue in the first approach/graph is WHAT is being simulated, which is tricky with mixed models. In the first approach, I think we are seeing a simulation that applies to the fixed effects, perhaps accounting for uncertainty in the random effects. Perhaps something like that, but more subtle, is happening in the second approach. Perhaps it's as simple as my doB3mc function being wrong -- I'll have to check that again. Any insights that come to you are appreciated. –  Wayne Feb 3 '12 at 19:33
I'm not familiar with mcmcsamp, I found this: rss.acs.unt.edu/Rdoc/library/lme4/html/mcmcsamp.html , have you tried deviance=TRUE in the command? I'm sort of clutching at straws here. –  Michelle Feb 3 '12 at 19:52