Skip to content

Syntax for a multivariate & multilevel MCMCglmm model

4 messages · Allan Debelle, Paul Debes

Original
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
Paul Debes
# 
Hi Allan,

This model (explicitly nesting random replicates within treatments):

fixed = var ~ 1 + treatment + cov1 + cov2,
random = ~ treatment:replicate,
rcov = ~ units

expanded to a multivariate model with unstructured covariances in G and R  
may be (you missed the "trait" interaction in the fixed statement of your  
first model):

fixed = cbind(varA, varB, varC, varD) ~ trait + trait:treatment +  
trait:cov1 + trait:cov2,
random = ~ us(trait):treatment:replicate,
rcov = ~ us(trait):units

Adding the trait interactions probably gets you back the treatment effects  
you did see in the univariate models.

How you add an additional random term, nested within replicates, will  
depend on how your random treatment:replicate effects are correlated  
between traits (e.g., "us" above may be reduced to "diag" if not  
correlated) and how the sub-replicates are correlated with the replicates  
and the traits. Plotting random effects and nested model testing may be  
required to find a useful solution.

Paul


On Fri, 04 Mar 2016 15:30:19 +0200, Allan Debelle <a.debelle at sussex.ac.uk>  
wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Allan Debelle
    

    

      
      #
    

  


  

      
Hi Paul,

Thank you very much for taking the time to have a look at this.

No need to modify the residual variance structure then? Just "rcov = ~ 
us(trait):units" ? And I don't need to change the syntax of the prior?

Thank you.

Al

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Paul Debes
    

    

      
      #
    

  


  

      
Hi Allan,

Maybe take a look at the residuals, if they show independence and a lovely  
distribution for each trait (or even across traits after standardisation)  
when plotted the "us(trait):units" works for you. You assume each trait  
has a different residual variance and residuals are correlated between  
trait pairs within individuals. You can, of course, relax this or make  
this more complicated - only model diagnostics and comparisons can tell  
you.

I don't dare giving advice on priors (being completely clueless myself  
once correlations are involved), but the "course notes" by Jarrod Hadfield  
are extremely informative on many modeling matters and also include  
information about covariance structures and priors. You can find them here:

https://cran.r-project.org/web/packages/MCMCglmm/vignettes/CourseNotes.pdf

Best,
Paul



On Fri, 04 Mar 2016 19:00:21 +0200, Allan Debelle <a.debelle at sussex.ac.uk>  
wrote: