[R-sig-ME] MCMCglmm: predictions for posterior predictive checks

[Jarrod Hadfield]

Hi,

post.pred is equivalent to y^{rep} but you have to be careful about  
whether you want to marginlaise the random effects or not. If you  
don't marginalise the random effects then post.pred is usually the  
joint posterior predictive distribution of y (where I take y to be a  
vector, and joint refers to the fact that it is the joint distribution  
of y^{rep}_{1},y^{rep}_{2} ...y^{rep}_{n}). However, post.pred always  
marginalises the residuals (or it wouldn't be a predictive  
distribution!).  If, for example, the model was bivariate with  
observations x and y, and R-structure of the form us(trait):units then  
post.pred does not give the joint posterior predictive distribution of  
x and y: it will give the joint posterior predictive distribution of x  
marginal to the residuals of y, and vice versa.

If the random effects are marginalised then post.pred is returning the  
marginal posterior predictive distribution for each observation. This  
might not always be what you want.  Imagine you make two observations  
on the same individuals at each of two time points, and you have the  
random effect structure us(time):individual.  This fits different  
between-individual variances for observations at different times, and  
also fits a between-individual across-time covariance. If you  
marginalise the individual effects (for example if you want to get a  
posterior predictive distribution for an observation made on a random  
individual) then post.pred will not preserve the covariation structure  
of the four observations. For example, have Y_ijk as the kth  
observation at time j for individual i. We have Y_i11, Y_i12 Y_i21,  
Y_i22.  Imagine us(time):individual gives us a covariance matrix V  
with 1 and 2 along the diagonal, and sqrt(1/2) on the off-diagonals  
(ie. a correlation of 0.5). If you marginalise the individual effects  
then the they are assumed to be drawn from a distribution with  
diagonal variance structure Diag{1,1,2,2}, whereas perhaps you would  
prefer them to be drawn from kronecker(V, diag(2)).

This is a bit hard to explain! If it doesn't make sense let me know,  
and I'll try again.

Cheers,

Jarrod














Cheers,

Jarrod


Quoting Ramon Diaz-Uriarte <rdiaz02 at gmail.com> on Sat, 01 Mar 2014  
12:44:04 +0100:

> Dear List,
>
> I want to perform a simple posterior predictive check on some Poisson
> models I've fitted using MCMCglmm. I immediately though about using
> predict.MCMCglmm as
>
> predict(mymodel, type = "response", marginal = 0, interval = "prediction")
>
> However, predict returns the expectation (in fact one of its very last
> lines have pred <- matrix(colMeans(post.pred), dim(post.pred)[2], 1)).
>
>
> I can hack predict.MCMCglmm and directly return the "post.pred" object
> which, IIUC, would give me the "y^{rep}" (in Gelman et al., 1996,
> notation.).  But having to do this makes me wonder if I am understanding
> this correctly.
>
> Is directly using the "post.pred" object the right way of getting the
> y^{rep} with MCMCglmm?
>
>
> Best,
>
>
> R.
>
>
>
>
> P.S. I am using "marginal = 0" as I want what, e.g., Green et al., 2009
> ("Use of posterior predictive assessments to evaluate model fit in
> multilevel logistic regression", Vet. Res, 40) call "full": "The predictive
> distribution was conditional on all fixed effect and random effect
> parameters estimated in the final model and a replicate observation
> y_{ijk}^{full} generated from the conditional distribution (...)".
>
> --
> Ramon Diaz-Uriarte
> Department of Biochemistry, Lab B-25
> Facultad de Medicina
> Universidad Autónoma de Madrid
> Arzobispo Morcillo, 4
> 28029 Madrid
> Spain
>
> Phone: +34-91-497-2412
>
> Email: rdiaz02 at gmail.com
>        ramon.diaz at iib.uam.es
>
> http://ligarto.org/rdiaz
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.