[R-sig-ME] Multiresponse MCMCglmm model specification and questions

[Jarrod Hadfield]

Hi Rafael,

1/ The reason that nu=1.002 is sometimes used for a 2x2 covariance  
matrix is because the marginal prior for one of the variances is  
equivalent to an inverse-Wishart with nu=0.002. You would also have to  
rescale V to get an inverse-gamma distribution for the marginal  
variance, as given in the CourseNotes footnote.

2/ This is the correlation in the phyloegentic effects, after  
accounting for the fixed effects. It may also make sense to have the  
residual structure as us(trait):units to model any residual  
correlation between the traits. The chain can get `stuck' at extreme  
correlations (-1 or 1) as it can for variances close to zero. The best  
way to see if this is a problem is to plot the MCMC time-series for  
the correlation (not the covariance) and see if it looks horrible.  
Parameter expansion can alleviate this problem, but the resulting  
priors on the covariance matrix are not inverse-Wishart but scaled  
Multivariate-F.

3/ Tricky  - and demanding of the data.

us(at.level(X1,1):trait):spp+us(at.level(X1,2):trait):spp

Fits two phylogentic covariance matrices, one for each of the two  
levels of X1. However, it assumes that there are no phylogenetic  
correlations in the response
between species that have different levels of X1.

us(X1:trait):spp

fits the phylogenetic covariances between the different levels of X1 too.

However, unless you have a very large phylogeny I would consider  
sticking to the simpler model. In addition, MCMCglmm cannot fit block  
diagonal structures for the residual (co)variances, and I think you  
would want to do this to ensure that the phylogenetic covariances  
weren't actually due to a mis-specified residual structure.  A  
saturated residual variance structure would be 2 2x2 covariance  
matrices, one for each level of X1. A single 4X4 covariance matrix  
(us(trait:X1):units) can be fitted, but the 2 2x2 off-diagonal blocks  
(the residual covariances between observations made on different  
levels of X1) cannot be estimated (no observation is associated with  
*both* levels of X1). As discussed elsewhere (I think a post by Celine  
Teplitsky) this may cause problems.

Cheers,

Jarrod





Quoting Rafael Maia <queirozrafaelmv at yahoo.com.br> on Fri, 17 Jan 2014  
12:05:16 -0500:

> Dear list members,
>
> I am attempting to run a multiresponse model to investigate factors  
> influencing the evolution of a trait measured in both males and  
> females of multiple species. Therefore I chose a multiresponse  
> approach using MCMCglmm, but my model has a series of little caveats  
> (don’t they all?) that I’d like to make sure I’m addressing  
> correctly. Here is a brief description of the model I’m trying to fit:
>
> Ymale, Yfemale: the same response variable, measured in males and females
>
> X1: a variable (factor, 2 levels) that I expect to affect both sexes  
> differently
> X2: a variable (factor, 2 levels) that I expect to affect both sexes  
> equally (or estimate the joint effect across sexes)
> X3: a continuous variable measured only for males (which I therefore  
> expect to only have an effect in male measurements)
>
> so I am currently specifying the model as below:
>
> mod <- cbind(Ymale, Yfemale) ~ trait*X1 + X2 + at.level(trait, 1):X3
>
> ainv <- inverseA(phylo)$Ainv
>
> Prior <- list(R=list(V=diag(2), nu=0.002), G=list(G1=list(V=diag(2),  
> nu=0.002)))
>
> result <- MCMCglmm(mod,
>     random = ~us(trait):spp, rcov = ~idh(trait):units,
>     data=dat, prior = Prior, family = rep(‘gaussian’, 2),
>     nitt=2100000, burnin=100000, thin=1000)
>
> with that, I have a couple questions (besides: does this model  
> specification look reasonable?):
>
> 1. Prior: I wanna make sure I’m specifying an inverse-gamma with  
> scale=shape=0.001 for the variance components. Based on footnote 1  
> of the Course Notes (p.102) I think I am based on the changes, but  
> I’ve seen several times in mailing lists responses and other papers  
> nu=1.002 being used for a matrix of same dimensions. My results are  
> robust to either specification but I’d like to make sure I’m doing  
> what I think I’m doing, given that there have been changes on this  
> regard.
>
> 2. I am calculating the correlation between ymale and female as:
>     result$VCV[,2]/sqrt(result$VCV[,1]*result$VCV[,4])
> as described in the tutorial vignette. My values, however, are  
> extremely high. The two traits are indeed strongly correlated among  
> species, but I’m getting credible intervals of 0.96 - 0.99 (for  
> comparison, a raw non-phylogenetic correlation of the traits is  
> 0.75). Am I forgetting to include any variance components to  
> calculate the intersexual correlation? Is this increase in the  
> correlation expected given the inclusion of the fixed and random  
> effects in the model?
>
> 3. I have a hypothesis that the correlation between male and female  
> traits should be stronger for one level of X1 than for the other. Is  
> it possible to specify random terms in order to calculate  
> (co)variances of the traits conditional on X1?
>
> I deeply appreciate any help you can provide. Cheers!
>
>
> Abraços,
> Rafael Maia
> ---
> http://www.rafaelmaia.net/
> PhD Candidate, Integrated Bioscience
> University of Akron
> "A little learning is a dangerous thing; drink deep, or taste not  
> the Pierian spring." (A. Pope)
>
> _______________________________________________
> 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.