[R-sig-ME] Syntax for a multivariate & multilevel MCMCglmm model
Allan Debelle
a.debelle at sussex.ac.uk
Fri Mar 4 14:30:19 CET 2016
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]]
More information about the R-sig-mixed-models
mailing list