[R-sig-ME] Multi-response model to estimate Correlation

[Ned Dochtermann]

Hi Xavier,

I'm far from an expert (to put it graciously) so it might be best to wait on
some of the heavy hitters but for the correlation you might consider
specifically estimating the posterior mode for the correlation. What you are
doing now is estimating the correlation from the posterior modes of the
covariance and two variances, which is slightly different. Hopefully someone
can weigh in on whether one is actually "better" than the other (as an
aside, it would be nice if you reported the unstandardized covariance and
variance posterior modes in addition to the correlation).

to return the correlation:
       posterior.mode(posterior.cor(multi3prior2yr $VCV[,1:4])[,2])

credibility interval:
       HPDinterval (posterior.cor(multi3prior2yr $VCV[,1:4])[,2])

Presumably the posterior.cor function effectively converts the covariance to
a correlation for each MCMC sample, from which you can then draw the
appropriate values. If your burn-in was sufficient and thinning appropriate
I doubt this makes too much of a difference but it may make a small one.

On the topic of priors, I don't know much and have had a very difficult time
finding approachable material on the topic. The MCMCglmm program seems to
also conflate the shape of the distribution with the "confidence" placed on
the distribution within a single parameter which makes it even harder for me
to follow. I started a thread awhile back listing out some prior
specifications I had gleaned from the listerv:
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q2/006088.html
As you can see I had to back off the topic in the second post due to my
ignorance and unfortunately the topic wasn't picked up by others.

Good luck!
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
--



------------------------------

Message: 3
Date: Tue, 21 Jun 2011 14:19:36 +0100
From: Xavier Harrison <xav.harrison at gmail.com>
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] Multi-response model to estimate Correlation
	between	Traits using MCMCglmm
Message-ID: <255B0F83-79E0-4338-897D-6BB5956FC653 at gmail.com>
Content-Type: text/plain

Dear List

I'm studying Heterozygosity-Fitness correlations in a migratory bird, and
have used MCMCglmm to fit Poisson models of the form "Reproductive
Success~Heterozygosity". These models run well, and show a positive effect
of heterozygosity on number of offspring produced. However, one reviewer has
asked that I fit a multi-response model where no. of juveniles and
heterozygosity are run as a single response, rather than outcome and
predictor, so that I may estimate the correlation between them. I list my
analysis below because I would greatly appreciate some advice to ensure I've
gone about this the correct way, as I would like this to be watertight when
this goes back out to review.

To do so, I have largely followed the guidelines from the course notes (~ pg
90) and from this post here:
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/005694.html

My model is as follows:

multi3prior2yr<-MCMCglmm(cbind(juv,ir)~trait-1,random=~us(trait):tempcode,rc
ov=~us(trait):units,data=yr2,family=c("poisson","gaussian"),verbose=FALSE,pr
ior=multiprior3,burnin=20000,nitt=120000,thin=100)

"tempcode" is the indicator variable referencing members of breeding pairs,
as sometimes both the male and female have been sampled, and will have the
same reproductive success but different heterzygosities. "juv" is number of
offspring and "ir" is heterozygosity, calculated as internal relatedness
following Amos (2001).

The priors I have used for these models are as follows. As far as I can see
the choice of prior doesn't seem to affect the estimates too much, which I
take as a good sign, although when using sets 2 and 3 the estimate
attenuates slightly, which I expected as I read something about how set1 of
priors might actually be fairly informative for this kind of model. I'm not
certain, however, that I've specified the "uninformative" priors correctly.
For these models, autocorr stats all look good (<0.1 between stored
iterations, as suggested by the course notes).

multi.prior<-list(R=list(V=diag(2),nu=2),G=list(G1=list(V=diag(2),nu=2)))

multiprior2<-list(R=list(V=diag(2), nu=2),G=list(G1=list(V=diag(2), nu=2,
alpha.mu=c(0,0), alpha.V=diag(2)*1000)))

multiprior3<-list(R=list(V=diag(2)*1e-6,
nu=1.002),G=list(G1=list(V=diag(2)*1e-6, nu=1.002, alpha.mu=c(0,0),
alpha.V=diag(2)*1000)))


To assess the correlation between outcome variables, I follow the example in
the link above where the vcov matrix is extracted from the model object and
then converted to a correlation using the "cov2cor" function. I've written a
basic function to do this for any model object from this set of models

mc.corr<-function(mcobject){
xx<-matrix(c(posterior.mode(mcobject$VCV)[1],posterior.mode(mcobject$VCV)[2]
,posterior.mode(mcobject$VCV)[3],posterior.mode(mcobject$VCV)[4]),2)
return(cov2cor(xx))
}

This returns a matrix of the form:

mc.corr(multi3prior2yr)
           [,1]       [,2]
[1,]  1.0000000 -0.2760332
[2,] -0.2760332  1.0000000

which I believe to be the correlation between heterozygosity and
reproductive success at the group level (breeding pair) rather than at the
unit level.

This is in agreement with the mixed effects regression where increasing
homozygosity (IR increases) causes a reduction in number of juveniles, so at
least superficially appears correct. My concern is that the nature of my
model parameterisation may be overestimating the degree of correlation.
Taking a step back from that, is the "cov2cor" function an acceptable way of
estimating this correlation?

Any advice regarding model or prior specification is greatly appreciated.

Many thanks

Xavier Harrison


---------------------------------------------------------
Dr Xavier Harrison
BBSRC Postdoctoral Research Associate
Centre for Ecology and Conservation
University of Exeter
Tremough Campus
Penryn
Cornwall
TR10 9EZ
Tel: (01326) 371872
xav.harrison at gmail.com