Thanks Jarrod, this seems exactly what I need.
I'm trying to implement it now, but for some reason I get an error
prior$R<-list(V=diag(3), nu=0, beta.mu=rep(0,3), beta.V=diag(3)*100)
prior$R$beta.V[1,1]<-1e-8
mod1=MCMCglmm(value~(letter-1), random=~us(letter):id,
rcov=~ante2(letter):day,family="gaussian", data=ALL2,
nitt=510000,burnin=10000, thin=500,verbose=FALSE, prior=prior,pr=TRUE)
Error in MCMCglmm(value ~ (letter - 1), random = ~us(letter):id, rcov =
~ante2(letter):day, :
R-structure miss-specified: each residual must be unique to a data point
An example of the data for two individuals:
id value letter day
1 379 -0.52909133 mirror 2
2 379 -0.36287842 mirror 4
3 379 -0.43505872 exploration 1
4 379 -0.31524632 exploration 4
5 379 -1.12034977 boldness 5
6 379 -0.15906279 exploration 5
7 379 -0.50567659 boldness 3
8 379 -0.51900960 exploration 3
9 379 0.03790821 exploration 2
10 380 -0.18712819 mirror 2
11 380 -0.05894661 exploration 5
12 380 -0.77426854 mirror 4
13 380 -0.61255224 exploration 4
14 380 -0.65062409 exploration 3
15 380 -0.06765472 exploration 2
16 380 -0.21452781 exploration 1
17 380 -0.72320719 boldness 3
18 380 -0.41338280 boldness 5
Can you see where the error comes from?
Thanks
2015-10-08 17:27 GMT+02:00 Jarrod Hadfield <j.hadfield at ed.ac.uk>:
Hi,
The appropriate residual syntax would be:
rcov=~idh(letter):day
although I would estimate the residual covariances, given that you assume
the individual-level random effects for the three traits can be correlated.
However, this raises an issue because trait 2 and trait 3 have never been
measured on the same individual on the same day and so you cannot estimate
their residual covariance. The covariance matrix looks like:
x x x
x x 0
x 0 x
where I have x's for things that can be estimated and 0 for the things
that cannot. Older versions of MCMCglmm didn't allow you to fit this type
of constraint (i.e. in the absence of information fix the unidentifiable
covariance to zero). However, you can trick it into doing this using an
antedependence model. Reorder your traits so trait 1 is last, and so the
constraint you want is:
x 0 x
0 x x
x x x
Then fit a second order autoregressive model:
rcov=~ante2(letter):day
and fix the regression of trait 3 on trait 2 to be zero (with the new
ordering this would be 2 on 1, i.e. the first regression coefficient).
prior$R<-list(V=diag(3), nu=0, beta.mu=rep(0,3), beta.V=diag(3)*100)
prior$R$beta.V[1,1]<-1e-8
To see why this works, define the matrix B where B[i,j] is the regression
of j on i. In this case B has the form:
. . .
0 . .
x x .
where 0 is the regression coefficient set to zero, and . are regression
coefficients that are zero by design in an antependence model. The
residual covariance matrix has the same sparse structure as
solve(I-B)%*%t(solve(I-B)) which is
x 0 x
0 x x
x x x
as desired. I believe (could be wrong, so it would be good to get
confirmation) that the only other constraint imposed on the x's is that
they result in a positive definite matrix irrespective of the true value of
the covariance between trait 2 and trait 3.
Note that MCMCglmm will return the residual covariance matrix. If you want
the matrix in terms of regression coefficients (and innovation variances)
you can use posterior.ante(mod1$VCV[,9:18], k=2).
Note that you may want to increase/decrease the prior variance on the
identifiable regression coefficients depending on scale, and you may not
want flat improper priors on the innovation variances.
Cheers,
Jarrod
Quoting David Villegas R?os <chirleu at gmail.com> on Thu, 8 Oct 2015
11:47:12 +0200:
Hi all,
I'r trying to run a multivariate MCMCglmm with 3 traits. I was suggested
to
use the long format since I have unequal number of replicates per trait.
Trait 1 and 2 were replicated twice, trait 3 was replicated five times.
Traits are gaussian.
The way I measured the trait for each individual is as follows:
day 1: trait1
day 2: trait1 and trait 2
day 3: trait1 and trait 3
day 4: trait1 and trait 2
day 5: trait1 and trait 3
From the model, I'm interested in extracting the between-individual
variances/covariances and if possible, the within-individual
variances/covariances.
This is my attemp so far. Letter identifies the trait.
mod1=MCMCglmm(value~(letter-1), random=~us(letter):id,
rcov=~idh(letter):xxxx, family=c("gaussian"), data=ALL)
My questions are about the rcov bit, the residual variances/covariances...
- First I don't know if with my experimental design, it makes sense to
estimate residual covariances ("us" structure) or constrain them as in the
model above ("idh" structure)
- Second, I don't know how to define the xxxx variable according to my
experimental design.
I guess both questions are related.
Any advise will be appreciated.
Thanks
David
[[alternative HTML version deleted]]