prediction intervals from a mixed-effects models?

      How can I get prediction intervals from a mixed-effects model?
 Consider the following example:

 library(nlme)
 fm3 <- lme(distance ~ age*Sex, data = Orthodont, random = ~ 1)
 df3.1 <- with(Orthodont, data.frame(age=seq(5, 20, 5),
                    Subject=rep(Subject[1], 4),
                    Sex=rep(Sex[1], 4)))
 predict(fm3, df3.1, interval='prediction')
 #      M01      M01      M01      M01
 # 22.69012 26.61199 30.53387 34.45574

 # NOTE:  The 'interval' argument to the 'predict' function was ignored.
 # It works works for an 'lm' object, but not an 'lme' object.

      One way to do this might be via mcmcsamp of the corresponding
 'lmer' model:

 library(lme4)
 set.seed(3)
 samp3r <- mcmcsamp(fm3r, n=10000)
 samp3r[1:2,]

      Then use library(coda) to check convergence and write a function
 to simulate a single observation from each set of simulated parameters
 and compute quantile(..., c(.025, .975)) for each prediction level
 desired.

      However, before I coded that, I thought I would ask if some easier
 method might be available.

      Thanks,
      Spencer
 p.s.  RSiteSearch("lme prediction intervals") produced 3 hits including
 2 from James A Rogers over 3 years ago.  In one, he said, "I am not
 aware of any published R function that gives you prediction intervals or
 tolerance intervals for lme models."
 (http://finzi.psych.upenn.edu/R/Rhelp02a/archive/42781.html)   In the
 other, he provided sample code for prediction or tolerance intervals of
 lme variance components.
 (http://finzi.psych.upenn.edu/R/Rhelp02a/archive/44675.html)