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

[Jarrod Hadfield]

Hi,

i) The syntax for your response is incorrect it should be cbind(Yes,  
No, Pred) (but I guess you had this right when you ran it).  If the  
multinomial response is Bernoulli it would be slightly quicker to have  
it as a single column of Yes/No's and have the family argument as  
"categorical".


ii) your output from table(data$Pred, data$Yes) looks like the levels  
are not ordered properly as it suggestes High>Some>Low rather than   
High>Low>Some (but perhaps this intentional?)

iii) although the residual variances are not identified for either  
trait, the residual correlation is. I would advise estimating it:   
rcov=~cor(trait):units.

iv) I would suggest saving the latent variables (pl=TRUE) and checking  
that they are not outside of the range -7/7 or -20/20 for the ordinal  
variable or categorical/multinomial variables respectively. If they  
are the model cannot be trusted because the probit/logit functions are  
under/overflowing.

v) prior2 should behave relatively OK (although having nu=2 (us) or  
nu=1 (idh) rather than nu=3 is more usual). However, it does put  
weight on a phylogenetic H^2 of 1 which could be a problem if there is  
high support in this area.


vi) the ordinal part is not mixing at all well: cutpoint.traitPred.1  
has zero effective sample size.  I would try and diagnose the problem  
with a univariate analysis of the ordinal response - perhaps your  
phylogeny has two clades; one clade containing species with all  
"High"/"Some" and the other clade containing species with all "Low"?

Cheers,

Jarrod





Quoting Mirjam Amcoff <mirjam.amcoff at ebc.uu.se> on Sun, 15 Apr 2012  
19:32:12 +0200:

> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.