[R] MCMCglmm priors including phylogeny

[Stott, Iain]

Hi all,

I'm hoping I might be able to get some help with some issues specifying priors for MCMCglmm.

I'm trying to fit a gaussian glmm using MCMCglmm to a data set with two (correlated) response variables.  The response variables are both logit-transformed proportions (there are a few reasons why I've chosen these with gaussian error over binomal glmm, which I won't  go into).  Each probability represents a response for a single plant species (c.150 species), and I want to incorporate the phylogeny (almost fully resolved with only 2 polytomies, including branch lengths) in the error structure.  There is only 1 explanatory variable, a 5-level factor.  I'm using uninformative priors for fixed parameters (i.e. not specifying any priors).

Here's my code:

#prior model1
prior1<-list(R = list(V = diag(2)*x, n=2),
             G = list(G1=list(V = diag(2)*x, n=2)))

#model1
m1<-MCMCglmm(cbind(y1.logit,y2.logit) ~ trait:Ecotype - 1, 
             random=~us(trait):animal, rcov=~us(trait):units, family=c("gaussian","gaussian"),
             prior=prior1, data=meandata, pedigree=meantree, nodes="TIPS",
             thin=100, nitt=150000, burnin=30000, verbose=F)

#prior model2
prior2<-list(R = list(V = diag(2)*x, n=2))

#model2
m2<-MCMCglmm(cbind(reac.prob.logit,inert.prob.logit) ~ trait:Ecotype - 1, 
             rcov=~us(trait):units, family=c("gaussian","gaussian"),
             prior=prior2, data=meandata, 
             thin=100, nitt=150000, burnin=30000, verbose=F)

I've tried prior variance values of 5, 10, 50, 100, 500, 1000, 5000 and 10000 (by varying x in the code above), and I've used a large thinning interval (100) to deal with autocorrelation arising from the correlated responses. Chains seem to mix well as long as I use enough iterations (over 100000 seems to be best).  Posterior distributions from plot(model) are normal-distributed for fixed effects and residuals in every case.  However, posteriors for errors don't behave as nicely:
- Lower prior values (superficially <100) result in positively-skewed distributions for trait:animal whereas higher prior values (superficially >100) are much more normal.  With the logit data running from ~-3.5 to ~10, it seems that the prior values that result in well-fitting models are higher than the variance of the data itself.  
- Fixed parameter values vary (but not considerably) with choice of prior, but error parameters seem to be more sensitive to prior choice.  
- Posterior means for errors and residuals are larger and more variable for y1:y1 and y2:y2 combinations (in the hundreds) but smaller for y1:y2 and y2:y1 combinations (in the tens).

Considering this, my first question is: what are appropriate priors for G and R?

I'm aware that my priors assume a priori independence between the two responses, and I know that they will covary, but my understanding is that specifying us(trait) in my errors and residuals in my residuals should estimate covariance anyway.  But, my second question is this: what are the merits in using a multi-response model?  The phylogenetic models are taking a fair amount of time to run, and I hope that without the need to estimate covariances, they might be a bit speedier.

At some point, I will want to incorporate replication within species in the analysis (populations nested within species), but for now I'm really just interested in looking at mean response per species.  I'm assuming that priors will probably be of similar value (but more of them) once I include within-species replication... or will I have to go through this again??


Any help appreciated, thanks in advance

Iain


--------------------------------------------------------------------------------
Iain Stott
Centre for Ecology and Conservation
University of Exeter, Cornwall Campus
Tremough
Treliever Road
Penryn
Cornwall
TR10 9EZ
Tel: 01326 371852
http://biosciences.exeter.ac.uk/cec/staff/postgradresearch/iainstott/
--------------------------------------------------------------------------------