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

[Ned Dochtermann]

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