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]]