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

[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 <mailto:paul.debes at utu.fi>
> 4 March 2016 14:59
> 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 <mailto:a.debelle at sussex.ac.uk>
> 4 March 2016 13:30
> 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

-- 
Allan Debelle - PhD
Research Fellow at University of Sussex <http://www.sussex.ac.uk/>
Phone: +44 (0)748 095 3874 <callto:415-526-2339>
Web: sites.google.com/site/allandebelle
LinkedIn: fr.linkedin.com/in/allandebelle