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

[Jarrod Hadfield]

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.