[R-sig-ME] MCMCglmm bivariate random regression

[Jarrod Hadfield]

Dear Adam,

This is a complex model and I'm not sure MCMCglmm is going to be able  
to do all that you wish of it.

The key thing that you point out is that you want to estimate the  
covariance between the individual effect of a repeated measured trait  
and the residual of a non-repeated measure trait. In ASReml you can  
achieve this by having random=~us(individual) and  
rcov=~idh(trait):units, and then fix the residual (units) variance in  
the second (non-repeated measure trait) to zero.

In MCMCglmm this is not possible because the chain will not mix. I  
have suggested in the past that for these types of model the `units'  
variance could be left unfixed and the residual variance in the  
non-repeated trait estimated as the sum of the two variances, which is  
identifiable.  Alternatively the units variance could be fixed at  
something which is less than the minimum value of V_2-(C^2)/V_1  that  
has support in the posterior. Here V_2 is the residual variance in  
LBS, C is the covariance between residual LBS and the individual  
effect on body size, and V_1 is the variance of the individual effects  
on body size. By ensuring that the units variance is fixed at value  
less than this you will not have the problem that the correlation  
estimate hits the boundary conditions of -1/1.  I have not tried  
either of these solutions so you need to do some work to make sure  
what I say is sensible.

For the second problem I think you can switch syntax to  
~idh(EQ:trait):units and the LBS variances will be in position 5:8 and  
can be fixed. I keep meaning to allow direct products in the  
R-structure which should make this type of model easier to fit.

The model is complex so be careful!

Jarrod







Quoting Adam Hayward <a.hayward at sheffield.ac.uk> on Wed, 29 May 2013  
13:37:07 +0100:

> Dear list users,
>
> I wonder if anyone would be able to help me with a problem which I've been
> trying to solve for quite some time. I want to run a bivariate random
> regression, where the two dependent variables are body weight and lifetime
> breeding success. Body weight is measured multiple times per individual,
> while LBS is measured once. I then want to fit two random slopes in the
> individual-level variance component for body weight- I x A and I x E. The
> key is then to estimate the individual-level covariances involving LBS to
> look for evidence of selection on body weight and individual slopes of
> weight on A and E. Finally, I would like to fix the residual variance in
> LBS to close to zero (I appreciate that it is impossible to fix it to
> exactly zero), but I would like the residual variance in body weight to
> vary across four quartiles of E (EQ).
>
> I've been having two issues. Firstly, how to fix the residual variance in
> LBS to zero, while allowing that for body weight to vary freely. The
> problem is that when I specify "rcov = ~idh(trait:EQ):units", the resulting
> 8 x 8 VCV matrix (where only the variances are estimated), rather than
> havig elements 1-4 being the four residual variances in weight and elements
> 5-8 being those for LBS, has the weight variances at the diagonal locations
> 1,3,5, and 7 and those for LBS at 2,4,6, and 8. In the former scenario, it
> would be easy to fix the residual variances in LBS by using V = 0.001, fix
> = 5 in the prior specification, but since the weight and LBS residual
> variances alternate, this is not possible. I am yet to find a way to fix
> specific variances in the residual structure. I have been able to do this
> in ASReml, but the distribution of LBS is far from Gaussian, hence checking
> it in MCMCglmm.
>
> The second issue is that I have found it very difficult (even without
> fixing residual variances in any way) to obtain satisfactory diagnostics-
> in particular, the autocorrelation is very high for the individual VCV
> components involving LBS and the random slopes, even with a run as long as
> 10 million iterations. I though that treating LBS as zero-inflated could
> help, but have been unable it work out how to incorporate the necessary
> zipoisson prior specifications into a bivariate model.
>
> If anyone has any thoughts on fixing individual elements of the R
> structure, or on sue of zipossion errors in bivariate models, I'd be very
> grateful and pleased to hear from you.
>
> Many thanks and best wishes,
> Adam
>
> A "basic" model is shown below. Nothing is fixed, but it illustrates the
> structure of the model:
>
> prior.x0 = list (R = list (V = diag(8), n = 0.002),
>               G = list (G1 = list (V = diag(4), n = 0.002),
>                         G2 = list (V = 1, n = 0.002),
>                         G3 = list (V = diag(2), n = 0.002),
>                 G4 = list (V = 1, n = 0.002)))
>
> model.x0<-MCMCglmm(cbind(WT, LBS) ~ trait-1 + some fixed effects,
>       random = ~ us(trait + at.level(trait,1):E + at.level(trait,1):Age):ID
> + us(at.level(trait,1)):Year + us(trait):MumID + us(at.level(trait,2)):BY,
>       rcov = ~idh(trait:EQ):units, data = data, family = c("gaussian",
> "poisson"), nitt = , burnin = , thin = ,
>       DIC = TRUE, prior = prior.x0, verbose = TRUE)
>
>
>
>
>
>
>
>
>
>
> --
> Adam Hayward
> Post-Doctoral Research Associate
> Department of Animal and Plant Sciences
> Alfred Denny Building
> University of Sheffield
> Western Bank
> Sheffield S10 2TN
> UK
> http://www.huli.group.shef.ac.uk/adam-personal.html
> http://adhayward.wordpress.com/
> https://twitter.com/adhayward18
>
> 	[[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.