[R-sig-ME] residual covariance structure and long format data in MCMCglmm

[David Villegas Ríos]

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]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>
>
>
> --
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
>
>
>

	[[alternative HTML version deleted]]