Mixed model predictions using sim

Dear Mailing list

I'm using lme4 to conduct mixed models and then, in accordance with
several authors I'm using the library 'arm' to simulate from that model to
generate confidence intervals for my parameter estimates (see
Korner-Nievergelt, F. et al ., 2015. Bayesian data analysis in ecology
using linear models with R, BUGS and Stan. Elsevier).

I would like to know more about the relationship between the Credible
Intervals (CrI) and the Confidence Intervals (2.5,50 and 97.5% quantiles)
generated from model predictions.  I'm particularly interested to know why
the Credible Intervals (given in 'CrI.MLexamp.9' below), for values of the
intercept (i.e. when the predictor is zero) are reflected very accurately
in the quantiles from the Predictions ('MLexamp.9.Pred' as below) yet there
is a discrepancy between values predicted by the CrI and the 95% CI
quantiles from the Predictions at, e.g., 60 units of 'open' in the example
below.  How does one move between the CrI to the model predictions?  Why
the difference?  (In my real example, the discrepancy is more meaningful).

The data (thanks) I've used in this example is from:


http://jaredknowles.com/journal/2013/11/25/getting-started-with-mixed-effect-models-in-r

Any pointers would be much appreciated.

Thanks

Tom

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

The following script might take 30 seconds to run...

library(lme4)
library(arm)
lmm.data=read.table("
http://www.unt.edu/rss/class/Jon/R_SC/Module9/lmm.data.txt",
                       header = TRUE, sep = ",", na.strings = "NA", dec =
".", strip.white = TRUE)
MLexamp.9=lmer(extro ~ open  + (1 + open | school/class),
                  data = lmm.data)
#summary(MLexamp.9)
nsim=20000
MLexamp.9.sim=sim (MLexamp.9,n.sim=nsim); #
colnames(MLexamp.9.sim at fixef)=names(fixef(MLexamp.9))<mailto:
MLexamp.9.sim at fixef)=names(fixef(MLexamp.9))>

#generate CrIs.
CrI.MLexamp.9=as.data.frame(signif(apply(MLexamp.9.sim at fixef
,2,quantile,prob=c(0.025,0.5,0.975)),5))

MLexamp.9.Pred=data.frame(open=c(0,20,40,60))#realistic values and zero.
MLexamp.9.Pred.Matrix=model.matrix(~open,data=MLexamp.9.Pred)
MLexamp.9.Pred$fit=(MLexamp.9.Pred.Matrix%*%fixef(MLexamp.9))#

fitmat=matrix(nrow=nrow(MLexamp.9.Pred),ncol=nsim)#
for(i in 1:nsim) fitmat[,i]=(MLexamp.9.Pred.Matrix%*%MLexamp.9.sim at fixef
[i,])#

MLexamp.9.Pred$LowerCI=apply(fitmat,1,quantile,prob=0.025)
MLexamp.9.Pred$Middle=apply(fitmat,1,quantile,prob=0.500)
MLexamp.9.Pred$UpperCI=apply(fitmat,1,quantile,prob=0.975)
CrI.MLexamp.9
MLexamp.9.Pred
The Scottish Association for Marine Science (SAMS) is registered in
Scotland as a Company Limited by Guarantee (SC009292) and is a registered
charity (9206). SAMS has two actively trading wholly owned subsidiary
companies: SAMS Research Services Ltd (SC224404) and SAMS Ltd (SC306912).
All Companies in the group are registered in Scotland and share a
registered office at Scottish Marine Institute, Oban Argyll PA37 1QA. The
content of this message may contain personal views which are not the views
of SAMS unless specifically stated. Please note that all email traffic is
monitored for purposes of security and spam filtering. As such individual
emails may be examined in more detail.

        [[alternative HTML version deleted]]