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

[Jarrod Hadfield]

Hi,

Sorry, kronecker(V, diag(2)) should have been kronecker(V, matrix(1,2,2)).

Cheers,

Jarrod



Quoting Jarrod Hadfield <j.hadfield at ed.ac.uk> on Sat, 01 Mar 2014  
18:32:36 +0000:

> 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.
>
> _______________________________________________
> 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.