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

[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:

> 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]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


-- 
Paul V. Debes
DFG Research Fellow
University of Turku
Department of Biology
Finland

Email: paul.debes at utu.fi