[R-sig-ME] Prior for multi-response model MCMCglmm

[Mirjam Amcoff]

Hi all,

I am trying to do a phylogenetic test of correlated evolution between  
two traits using MCMCglmm but am having trouble with setting the priors.

My response variable consists of one multinomial variable (cbind(yes,  
no)) and one ordinal variable with three levels (low, some, high). I  
have a phylogeny with 61 species and the raw data are:

table(data$Pred, data$Yes)

       0  1
High  1  8
Low  24  1
Some  8 19

The full model I want to run looks like this:

MCMCglmm(cbind(cbind(Yes, No), Pred ~ trait-1, random =  
~us(trait):animal, rcov = ~us(trait):units, family = c("multinomial2",  
"ordinal"), pedigree = tree, prior = prior1, data = data, nitt =  
1100000, burnin = 100000, thin = 100)

I want to compare this model with one where the covariance is  
constrained to estimate correlated evolution between the two response  
traits:

MCMCglmm(cbind(cbind(Yes, No), Pred ~ trait-1, random =  
~idh(trait):animal, rcov = ~us(trait):units, family =  
c("multinomial2", "ordinal"), pedigree = tree, prior = prior1, data =  
data, nitt = 1100000, burnin = 100000, thin = 100)

As far as I have understood the following prior should be appropriate  
with the residual variance fixed to 1 and the random effect prior  
being uninformative for the covariance and parameter expanded to  
improve mixing.

prior1 <- list(R = list(V = diag(2), nu = 0.002, fix=1), G = list(G1 =  
list(V = diag(2)*1e-6 , nu = 3, alpha.mu=c(0,0), alpha.V=diag(2)*25^2)))

However, when I try with this prior the mixing is exceedingly poor and  
when I try to increase the number of iterations I get the following  
error message:

Error in MCMCglmm(cbind(cbind(Yes, No), Pred) ~ trait - 1, random =  
~us(trait):animal,  :
   ill-conditioned G/R structure: use proper priors if you haven't or  
rescale data if you have

I find this strange as the only difference is that I increased the  
number of iterations?

When I try with the following prior things look much better:

prior2 <- list(R = list(V = diag(2), nu = 0.002, fix=1), G = list(G1 =  
list(V = diag(2) , nu = 3, alpha.mu=c(0,0), alpha.V=diag(2)*25^2)))

The only difference is that I used V = diag(2) instead of V =  
diag(2)*1e-6. The mixing in the ordinal variable is still quite poor  
but I think I could get away from that by running the chain longer  
than I have tried so far.

I do not understand enough about how the priors work in these models  
so I would be very grateful if someone could help me and tell me if  
the second prior is OK or not and, if it's not, what can I do to  
improve the mixing using the first model?


Here is an example of the output I get from the test models I have run  
using the second prior, note that the effective sample size is still  
very low for the ordinal (Pred) trait:

  Iterations = 6599801
  Thinning interval  = 600001
  Sample size  = 30000

  DIC: 393.0197

  G-structure:  ~us(trait):animal

                    post.mean   l-95% CI u-95% CI eff.samp
Yes:Yes.animal        0.4458  3.053e-08    1.410    29465
Pred:Yes.animal      1.0145 -8.075e-01    3.653    22451
Yes:Pred.animal      1.0145 -8.075e-01    3.653    22451
Pred:Pred.animal   15.5067  5.446e-01   41.764     5068

  R-structure:  ~us(trait):units

                   post.mean l-95% CI u-95% CI eff.samp
Yes:Yes.units             1        1        1        0
Pred:Yes.units           0        0        0        0
Yes:Pred.units           0        0        0        0
Pred:Pred.units         1        1        1        0

  Location effects: cbind(cbind(Yes, No), Pred) ~ trait - 1

            post.mean l-95% CI u-95% CI eff.samp pMCMC
traitYes     -0.1474  -1.0387   0.6735    27585 0.706
traitPred   -0.8047  -4.7517   2.9530    10050 0.650

  Cutpoints:
                       post.mean l-95% CI u-95% CI eff.samp
cutpoint.traitPred.1    0.2244   0.2244   0.2244        0


Thank you very much in advance,

//Mirjam