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

[Jarrod Hadfield]

Hi David,

My fault. It should be ante2(letter):day:id.

Cheers,

Jarrod



Quoting David Villegas Ríos <chirleu at gmail.com> on Fri, 9 Oct 2015  
10:24:00 +0200:

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



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.