[R-sig-ME] MCMCglmm multi-variate response for animal model

[Shana Caro]

Hi,

I am trying to fit a bivariate response model for a phylogenetic meta-analysis. I am looking at two behaviours (y1 in offspring and y2 in adults), and trying to determine if certain life history traits (x1 and x2) influence those behaviours. Logically, I would expect y1 and y2 to covary within species, which is why I want to fit a bivariate response model instead of 2 univariate models.

I am not certain if my data are in an appropriate format for a bivariate MCMCglmm, or if the model I am trying to run is actually telling me what I think it is. Apologies for a very beginner-type question.

My data is in this form:

Animal Study y1 y2 sample size x1 x2
species1 study1 0.05 n/a 20 yes poor
species1 study2 n/a 0.20 35 yes good
species2 study3 n/a -0.10 15 no normal
species2 study3 n/a -0.50 25 no normal
species3 study4 0.75 n/a 40 no good

…etc

Y1 and y2 are Z-transformed correlations. The life history traits are factors: x1 is a categorical yes/no, and x2 is an ordered categorical poor/normal/good. Each row has either y1 or y2 (never both). Each species can have multiple effect sizes, from multiple studies. Not every species/study has data on both y1 and y2 (41 species have both, 32 have data on only y1, and 20 have data on only y2; 93 total species), but I believe MCMCglmm will impute data (missing at random). Most species have 1-3 studies, and most studies have 1-5 effect sizes. I cannot collapse my data to 1 value per species or study, because x1 and x2 vary within species and study. Not every species has an effect size for every combination of x1:x2. N is the sample size of the study, and I am using 1/(n-3) for the mev argument.

Essentially, the outputs I want are:

  *   The contrast between the levels of x1 (for y1 controlling for y2, and vice versa)
  *   Whether the slope of x2 is different than 0 (for y1 controlling for y2, and vice versa)
  *   Whether there is an interaction between x1 and x2 (for y1 controlling for y2, and vice versa)
  *   The species-level correlation between y1 and y2 (controlling for x1 and x2).

The model I’ve run is:

prior = list(R=list(V = diag(2)/3, nu = 2), G=list(G1 = list(V = diag(2)/3, nu = 2), G2 = list(V = diag(2)/3, nu = 2)))

Model <- MCMCglmm(cbind( y1, y2 ) ~ trait * x1 * x2 - 1,
                         random = ~ us(trait):animal + us(trait):study,
                         rcov = ~ us(trait):units,
                         prior = prior,
                         mev = variance,
                         pedigree = tree,
                         data = data,
                         family = c("gaussian","gaussian"), verbose=F, pr=T, slice=T, nitt=600000, burnin=100000, thin=100)

Am I correct in this interpretation of the post.means for this model?

  *   Trait y1: the intercept for y1
  *   Trait y2: the intercept for y2
  *   x1: the contrast for x1.no for y1
  *   x2 (linear): the slope for y1 from x2
  *   Trait y2:x1: the contrast for x1.no for y2
  *   Trait y2:x2 (linear): the slope for y2 from x2
  *   x1:x2 (linear): the interaction between x1 and x2 for y1
  *   Trait y2:x1:x2: the interaction between x1 and x2 for y2

Thank you!

Shana Caro

Shana.Caro at zoo.ox.ac.uk




	[[alternative HTML version deleted]]