[R-sig-ME] MCMCglmm and prior specification

[Jarrod Hadfield]

Dear Cristina,

I am currently writing a better guide to MCMCglmm which deals more  
extensively with prior specifications - it should be ready soon.

In your first prior you have fixed the residual variance to be V=1, so  
it's not surprising that this gives very different answers from models  
2 and 3.

Often variance components are sensitive to prior information, in part  
because the data may not contain much information but also because the  
inverse-Wishart prior used is not really flexible enough (again, I am  
working on this).

The various improper priors often lead to numerical problems which  
precludes their use, but the resulting posterior does have some useful  
properties:

when V=anything & nu=0 the joint posterior modes  should be equal to  
the ML estimate
when V=0 & nu=-2 the joint posterior modes should be equal to the REML  
estimate (in practice you will have to use V=1e-6, or something small)

A commonly used proper prior is V=1 & nu=0.002.  This is equivalent to  
the inverse gamma distribution with shape=scale=0.001. This was the  
"default" in WinBUGS for a while before Andrew Gelman showed why it  
shouldn't be - it can be informative if the variances are close to zero.

Hope this helps,

If anyone wants an unfinished copy of the MCMCglmm guide just email me.

Jarod



On 24 Nov 2009, at 16:42, ledonret at email.unc.edu wrote:

> Dear all,
>
> I am trying to use the MCMCglmm package to create credibility  
> intervals for random variables in my data. I'm having a bit of  
> trouble though determining what the best prior to use for each model  
> is, since the results seem to differ tremendously depending on which  
> prior I am using, for instance, I've tried these three types of  
> priors,
>
>> halfFam<-var(data$Family)/2
>> prior1 
>> = 
>> list(R=list(V=1,n=1,fix=1),G=list(G1=list(V=1,n=1),G2=list(V=1,n=1)))
>> prior2 
>> =list(R=list(V=1,n=1),G=list(G1=list(V=1,n=1),G2=list(V=1,n=1)))
>> prior3 
>> = 
>> list 
>> (R 
>> = 
>> list 
>> (V 
>> =halfFam,n=1),G=list(G1=list(V=halfFam,n=1),G2=list(V=halfFam,n=1)))
>
> For the model,
>> model<-MCMCglmm(Length~1,random=~Family 
>> +Rep,data=data,verbose=FALSE,prior=prior,burnin=10000,nitt=75000)
>
> Where the random factors are Family and Replicate.
> From these priors, I get intervals for my Family effect,
>
>> HPDinterval(model1$VCV[,"Family"])
>         lower     upper
> var1 0.09660338 0.8888039
>
>> HPDinterval(model2$VCV[,"Family"])
>        lower    upper
> var1 0.1944570 2.120540
>
>> HPDinterval(model3$VCV[,"Family"])
>        lower    upper
> var1 0.2099238 1.529794
>
>
> I feel bad that I don't understand better how to specify the  
> components of these priors, but from what I understand, the model  
> should return similar values even if the priors are very different.  
> I've looked through the vignette thoroughly, but didn't get a sense  
> of what I was supposed to do if alternate priors returned different  
> answers. I'm not sure whether this is telling me that all the  
> information is coming from my priors (and there is, in fact, no  
> information in the data), or I am just incorrectly specifying my  
> priors.
>
> Any insight would be very much appreciated! Happy holidays,
>
> Cristina Ledon-Rettig
> UNC-Chapel Hill
>
> *I am using lme4 version 0.99375-28 with Mac OS X version 10.5
>
> _______________________________________________
> 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.