[R-sig-ME] MCMCglmm - detecting covariance patterns?

[James Gilbert]

Dear list,

I'm a relative newcomer to MCMCglmm and have a correspondingly newbie question.

I am trying to conduct a multivariate GLMM in order to ask whether a
continuous predictor affects the covariance between two response
variables, while at the same time incorporating non-independence due
to phylogeny.

I have 62 species of syrphid flies for each of which I have two
response variables, PCA1 and PCA2 -  calculated using phylogenetic PCA
(Liam Revell's method). n.b. these come from a wider dataset of 17
morphological characters from 259 species.

The 62 species are a subset of these for which I also have information
on total chromosome length (TCL, i.e. the sum of the lengths of the
chromosomes from the karyotype).

As well as asking whether TCL is related directly to PCA1 or PCA2 as a
fixed effect, I also want to ask whether TCL affects the strength of
correlation between these two principal components (which by
definition are uncorrelated in the wider sample, but may be related
among subsets of the species).

I would like to know whether it is appropriate for me to ask this
question by incorporating TCL as part of the random term in the model
or part of the residual term?

After reading the course notes and a few relevant posts on this list,
here are what I believe are my options:

EITHER incorporate the covariate into the random term, as follows:

str(pmdat)
'data.frame':	62 obs. of  5 variables:
 $ PCA1: num  30.75 -2.93 17.2 -31.58 27.81 ...
 $ PCA2: num  1.92 11.82 2.84 -1.96 6.35 ...
 $ TCL   : num  53.7 58.6 54.9 51.5 46.2 52.2 72.3 53.6 57 63.6 ...
 $ animal: chr  "Allograpta_obliqua" "Asemosyrphus_mexicanus"
"Cheilosia_burkei" "Cheilosia_illustrata" ...

phen.var<-matrix(c(var(pmdat$PCA1,na.rm=TRUE),0,0,
var(pmdat$PCA2,na.rm=TRUE)),2,2)

pr    <-list(G=list(G1=list(V=phen.var/3,nu=2, alpha.mu=c(0,0),
alpha.V=diag(2)*1000),
                    G2=list(V=diag(2)*2,nu=2, alpha.mu=c(0,0),
alpha.V=diag(2)*1000)),
                    R=list(V=diag(2)*2,n=3))

mcmod <- MCMCglmm(cbind(PCA1,PCA2) ~ (trait-1) * TCL,
                  random=~us(trait):animal + us(trait):TCL,
                  rcov=~us(trait):units, family=c("gaussian","gaussian"),
                  pedigree=tree.dat, nodes='TIPS',
                  data=pmdat, prior=pr,verbose=T, thin=100,
nitt=330000, burnin=30000)

TCL_effect_on_covariance <- posterior.mode(mcmod$VCV[,'PCA2:PCA1.TCL'])
TCL_effect_on_covariance_CI <- HPDinterval(mcmod$VCV[,'PCA2:PCA1.TCL'])


OR incorporate the covariate into the residual term, as follows -
although I am doubtful of my syntax for the rcov term:

pr    <-list(G=list(G1=list(V=phen.var/3,nu=2, alpha.mu=c(0,0),
alpha.V=diag(2)*1000)),
             R=list(V=diag(2)*2,n=3))

mcmod1<- MCMCglmm(cbind(PCA1,PCA2) ~ (trait-1) * TCL,
                  random=~us(trait):animal,
                  rcov=~us(trait:TCL):units,
                  family=c("gaussian","gaussian"),
                  pedigree=tree.dat, nodes='TIPS',
                  data=pmdat, prior=pr,
                  verbose=T, thin=100, nitt=330000, burnin=30000)

TCL_effect_on_covariance <-
posterior.mode(mcmod1$VCV[,'PCA2:TCL:PCA1:TCL.units'])
TCL_effect_on_covariance_CI <-
HPDinterval(mcmod1$VCV[,'PCA2:TCL:PCA1:TCL.units'])

Both models appear to run fine without any obvious issues, and I'm
just wondering which (if either) is more appropriate - and, if
neither, what would be the correct way to specify the model and
extract the effect I am interested in.

Many thanks in advance!

James Gilbert


--
t:+612 9351 4661  m:+614 5804 4038  w: http://jdjgilbert.wordpress.com
Heydon-Lawrence Bdg A08, Univ. of Sydney, NSW 2006, Australia
"I am dying by inches, for lack of anyone to talk to about insects.."
Darwin, diary entry