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

[Vincent Dorie]

Hi Tom,

If I'm following you correctly, the difference is due to the inclusion of
uncertainty in the slope term which is not present in the term for open !=
0. The quantiles are actually identical if you omit the signif call.

A more accurate representation of uncertain for new observations might be
to use the posterior predictive distribution. This can be obtained by
adding some rnorm using the sample of sigma (assuming that you were
interested in observations in a new group).

Vince

On Wed, Dec 9, 2015 at 5:18 AM, Tom Wilding <Tom.Wilding at sams.ac.uk> wrote:

> 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]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]