[R-sig-ME] Covariance between two traits in MCMCglmm
Jarrod Hadfield
j.hadfield at ed.ac.uk
Thu Mar 24 12:40:54 CET 2011
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.
More information about the R-sig-mixed-models
mailing list