[R-sig-ME] Syntax for a multivariate & multilevel MCMCglmm model

[Allan Debelle]

Hello everyone,

I am trying to analyse a dataset using a multivariate MCMCglmm model, 
but I cannot solve a syntax problem related to the structure of my dataset.

Basically, 50 individuals were exposed to either an A or a B treatment, 
and then 4 traits were measured on each individual, as well as 2 
covariates. This experiment was replicated 4 times, meaning that I have 
4 replicates (rep1, rep2, rep3 and rep4) for each of the two treatments 
(A and B).

The dataset looks something like this, with 'treat_rep' just being a 
unique identifier for each combination of treatment x replicate (I 
simplified it and generated random numbers, just FYI):

varA    varB    varC    varD    cov1    cov2    treatment    replicate   
  treat_rep

0.109488619    0.675591081    0.940580782    0.366980736    
0.911079042    0.468285642    A    rep1    A1
0.400887809    0.43162503    0.823351715    0.76152362    0.003221553    
0.622398329    A    rep2    A2
0.176948608    0.074676865    0.91156597    0.747663021    
0.028740661    0.856395358    A    rep3    A3
0.16468968    0.776755434    0.367321135    0.886352101    
0.083174005    0.246884218    A    rep4    A4
0.090596435    0.785823554    0.061103075    0.894436144    
0.989249957    0.50533554    B    rep1    B1
0.839349822    0.639784216    0.293986243    0.55812818    
0.307394708    0.036590982    B    rep2    B2
0.711702334    0.174145468    0.634296876    0.442112773    
0.480754889    0.863199351    B    rep3    B3
0.052352675    0.492622311    0.668647151    0.683243455    
0.958397833    0.857768988    B    rep4    B4
...

##############################

What I want to do is to test for the effect of the treatment on the 
traits I measured, while taking into account the effects of my 
covariates, so something like this:

varA + varB +varC +varD ~ treatment + cov1 + cov2

The problem is that I have to deal with a nested data structure. Each 
individual belongs to a replicate of a treatment (e.g. rep1 of treatment 
A, or rep3 of treatment B, etc.), and I am not sure how to specify this 
correctly in MCMCglmm. I tried a few things, but something is wrong 
either in the prior specification or/and in my specification of the 
variance structure of the random effects (and potentially of the 
residuals too).

##############################

Among other things, I tried to use the treat_rep variable as a grouping 
factor:

prior<- list(R=list(V=diag(4),nu=4),G=list(G1=list(V=diag(4),nu=4)))

   model <- MCMCglmm(

   fixed = cbind(varA, varB, varC, varD) ~ -1 + treatment + cov1 + cov2,

   random = ~ us(trait):treat_rep,

   rcov = ~ us(trait):units,

   prior = prior,

   family = c("gaussian", "gaussian", "gaussian", "gaussian"), nitt = 
100000, burnin = 5000,

   thin = 25, data = dataset)

However, and unless I am wrong, it is not nesting 'replicate' within 
'treatment'. Which means I end up with no effect of 'treatment', whereas 
I do have this effect when I run independent models using lme4 for each 
of the 4 traits (and in that case the nesting syntax is more 
straightforward).

##############################

I also tried this:

   prior<- 
list(R=list(V=diag(16),nu=16),G=list(G1=list(V=diag(8),nu=8))) # 4 
traits x 2 treatments x 4 replicates // 4 traits x 2 treatments

   model <- MCMCglmm(

   fixed = cbind(varA, varB, varC, varD) ~ -1 + treatment + cov1 + cov2,

   random = ~ us(trait:treatment):replicate,

   rcov = ~ us(trait:treatment:replicate):units,

   prior = prior,
   family = c("gaussian", "gaussian", "gaussian", "gaussian"), nitt = 
100000, burnin = 5000,
   thin = 25, data = dataset)

but I get the following error message:

Error in `[.data.frame`(data, , components[[1]]) :
   undefined columns selected

##############################

So, my first question is: What would be the syntax for this kind of 
nested structure (replicate nested within treatment) in MCMCglmm?

My second question is: What would be the syntax if there was another 
level of variance related to this nested structure? I'm thinking of the 
situation where there would be a sort of treatment nesting within each 
replicate (e.g. if each of the four pairs of replicates was kept in a 
different incubator, or if each pair had been done at a different moment 
in time, etc.). Would you just add a 'replicate' random effect in the 
model? How?

Any help with this would be fantastic.

Thanks a lot in advance.

Al

	[[alternative HTML version deleted]]