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,rcov=~us(trait):units,data=yr2,family=c("poisson","gaussian"),verbose=FALSE,prior=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@gmail.com










	[[alternative HTML version deleted]]


