[R-sig-ME] Mixed model predictions using sim

[Tom Wilding]

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]]