[R-sig-ME] Covariance matrix specification in MCMCglmm

Alexandre Courtiol alexandre.courtiol at gmail.com
Sun Mar 13 12:22:43 CET 2016

Dear all,

I am trying to figure out the syntax required to constrain a covariance
matrix in a particular way with MCMCglmm.

I have a multivariate problem, so let's say my response variables are "y1"
and "y2".
I have a fixed effect coding a treatment called "treatment" with "A" or "B".
I want the covariance matrix characterizing individual variation for each

So I wrote something along:

model <- MCMCglmm(
    fixed = cbind(y1, y2) ~  trait:treatment,
    random = ~ us(trait:treatment):individuals,
    rcov = ~ us(trait):units, ...)

and in the VCV object of the model, I can find something like:

|sigmaA1A1, sigmaA1A2, sigmaA1B1, sigmaA1B2|
|sigmaA2A1, sigmaA2A2, sigmaA2B1, sigmaA2B2|
|sigmaB1A1, sigmaB1A2, sigmaB1B1, sigmaB1B2|
|sigmaB2A1, sigmaB2A2, sigmaB2B1, sigmaB2B2|

which is almost what I want.

The only problem is that my individuals are never in both treatments, so I
would like to constrain the previous matrix to be:

|sigmaA1A1, sigmaA1A2,         0,         0|
|sigmaA2A1, sigmaA2A2,         0,         0|
|        0,         0, sigmaB1B1, sigmaB1B2|
|        0,         0, sigmaB2B1, sigmaB2B2|


In practice the covariances that I would rather not to estimate are close
to zero, which makes sense, so I could stick to my solution... but I could
gain a bit of extra power if I could constrain those to be exactly zero and
it would look a bit cleaner.

Any idea would be appreciated.



Alexandre Courtiol


*"Science is the belief in the ignorance of experts"*, R. Feynman

	[[alternative HTML version deleted]]

More information about the R-sig-mixed-models mailing list