[R-sig-ME] Covariance between two traits in MCMCglmm

[Jarrod Hadfield]

Hi Ned,

You should probably interact your fixed effects with trait, because at  
the moment you are assuming things like Sex have the same effect on  
P.TimeFront and M.TimeFront  - I'm not sure if this in intentional?

Regarding priors - the ones you have used are quite informative with  
respect to some parameters. For example, specifying  nu=2 as a prior  
for a 2x2 covariance matrix means that the marginal prior for each of  
the variances is inverse-Wishart with nu=1. This can have a strong  
effect if there are not much data.

 From experience I find

list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*a))

where a is something large (e.g. 1000, depending on the scale of the  
data) works well for the two standard deviations and the correlation,  
in terms of informativeness.

You can't use parameter expanded priors for the residual term yet, so  
you will have to stick with the standard inverse-Wishart (or use  
another program). Generally, data are highly informative for the  
residual part so often the posterior is not very sensitive to the  
prior specification. Nevertheless, you should check alternatives:

V=diag(2), nu=1.002 gives the inverse-gamma prior for the variances  
with shape=scale=0.001
V=diag(2)*1e-6, nu=3 is flat for the correlation from -1 to 1

Cheers,

Jarrod




On 23 Mar 2011, at 19:59, Ned Dochtermann wrote:

> List members,
>
> I'm currently trying to estimate the covariance and correlation  
> between two
> variables for which I have repeated measurements at the level of
> individuals. I'm getting output but would appreciate feedback  
> regarding my
> coding, if I've learned anything it is that getting output isn't the  
> same as
> getting output because you've done things correctly!
>
> BACKGROUND:
> I had conducted univariate analyses using glmer and the binomial  
> family, the
> variables are the amount of time spent essentially in one area versus
> another (based on some other validation work these variables equate to
> tolerance of predation risk by individuals and aggressiveness towards
> conspecifics). I then extracted the BLUPS from glmer, recognizing  
> that they
> were neither unbiased nor linear in this case, and calculated  
> Spearman's
> correlation (rs ~ 0.4). For well known reasons this is not a  
> preferred way
> to go so I wanted to calculate the correlation directly from a
> multi-response model. I know you can do this in glmer by including  
> one of
> the two variables as a covariate but I wanted to do it from a single  
> model.
>
> MCMCglmm MODEL:
> I chose a generally uninformative prior because I don't have any a  
> priori
> expectations for the variances but used a lengthy burnin. However,  
> to be
> honest, I don't have a good grasp of priors despite having read all  
> the
> MCMCglmm documentation and some other readings besides those.
>
> "mass.bar" and "mass.dev" represent mass, obviously, but split to
> distinguish between within versus among individual effects. I'm not  
> actually
> interested in the fixed effects from a multi-response model though,  
> just the
> correlation. Autocorrelation and other diagnostics look ok. Anyway,  
> the
> code:
>
> ####
> multi.prior<- 
> list(R=list(V=diag(2),nu=2),G=list(G1=(list(V=diag(2),nu=2))))
>
> corr.mcmc<- 
> MCMCglmm(cbind(P.TimeFront,P.TimeBack,M.TimeFront,M.TimeBack)
> ~(trait-1)+Site+Sex+mass.bar+mass.dev, random=~us(trait):Sub_ID,
> rcov=~us(trait):units,
> family=c("multinomial2","multinomial2"),
> data=Compiled.3,
> nitt=1080000,thin=480,burnin=120000,prior=multi.prior, verbose=FALSE)
>
> xx<-matrix(c(posterior.mode(corr.mcmc$VCV) 
> [1],posterior.mode(corr.mcmc$VCV)[
> 2],
> posterior.mode(corr.mcmc$VCV)[3],posterior.mode(corr.mcmc$VCV)[4]),2)
> cov2cor(xx)
> ####
>
> This produces an estimate (0.9) of the correlation between  
> P.TimeFront and
> M.TimeFront (TimeBack's are linearly dependent on TimeFronts) which  
> is the
> value I want. Unfortunately I can't find an example of a multi- 
> response
> binomial model in the documentation or in the archives of this  
> listserv so
> I'm a bit concerned that my specification of the response variables in
> particular, or other model terms, is incorrect. Any feedback on this  
> would
> be appreciated.
>
> Thanks a lot,
> Ned
>
> --
> Ned Dochtermann
> Department of Biology
> University of Nevada, Reno
>
> ned.dochtermann at gmail.com
> http://wolfweb.unr.edu/homepage/mpeacock/Ned.Dochtermann/
> http://www.researcherid.com/rid/A-7146-2010
> --
>
> _______________________________________________
> 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.