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