[R-sig-ME] MCMCglmm variance query

[Thomas Houslay]

Dear all,

I am using MCMCglmm to investigate genotype x environment interactions. I have a number of genetic lines, and two environments of interest; several individuals from each line were assigned to a single environment, and various characters were measured.

At the very simplest level, I just want to estimate the genetic correlation of a single character expressed in two different environments. The character I am using to get things going ('PE') fits a normal distribution, but I also have measurements which will require poisson / ZIpoisson once I’m a little more confident. I feel like this should be very straightforward, but I'm still having some problems.

I have used the ‘us’ variance function for the random effect of Line because I want to measure the covariance between measurements in L and H environments for each Line; ‘idh’ is used for the rcov term because no individual is measured in both environments. Envt  has two levels, 'L' and 'H'.

Here is my original model, with one of the various priors I tried using:

prior.3 <- list(R=list(V=diag(2),nu=0.002), G=list(G1=list(V=diag(2)*0.02,nu=3)))
model.1 <- MCMCglmm(PE~Envt-1, random=~us(Envt):Line, rcov=~idh(Envt):units, data=mydata, prior=prior.3, nitt=3550000, thin=500, burnin=50000, verbose=FALSE)

> summary(model.1)

 Iterations = 50001:3549501
 Thinning interval  = 500
 Sample size  = 7000 

 DIC: 1294.103 

 G-structure:  ~us(Envt):Line

         post.mean  l-95% CI u-95% CI eff.samp
H:H.Line     1.823  0.004461    8.528     7000
L:H.Line     1.687 -2.169747   11.880     7000
H:L.Line     1.687 -2.169747   11.880     7000
L:L.Line     2.738  0.003271   15.105     7000

 R-structure:  ~idh(Envt):units

        post.mean l-95% CI u-95% CI eff.samp
H.units     320.9    221.2    433.1     7000
L.units     288.4    201.9    383.3     7000

 Location effects: PE ~ Envt - 1 

           post.mean l-95% CI u-95% CI eff.samp  pMCMC    
EnvtH     47.02    42.74    51.36     7000 <1e-04 ***
EnvtL     46.59    42.52    50.50     7000 <1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

> posterior.mode(model.1$VCV)
     H:H.Line      L:H.Line      H:L.Line      L:L.Line       H.units       L.units 
-2.640172e-04  1.785168e-01  1.785168e-01 -1.828135e-03  3.201142e+02  2.699961e+02 


I think that this indicates that the variances are getting stuck at values close to zero...? I received similar results from all of the different priors which I used. 

Assuming this to be the case, I used the ever-excellent course notes and moved on to try parameter-expanded priors. I have tried various different forms (a couple of which I have given below), all of which give me similar results to that shown below:

## Parameter-expanded priors ##
## inverse-gamma prior for the variances with shape=scale=0.001
prior.parex.4 <- list(R=list(V=diag(2),nu=1.002), G=list(G1=list(V=diag(2),nu=2,alpha.mu=c(0,0),alpha.V=diag(2)*1000)))
## flat for the correlation from -1 to 1
prior.parex.5 <- list(R=list(V=diag(2)*1e-6, nu=3), G=list(G1=list(V=diag(2),nu=2,alpha.mu=c(0,0),alpha.V=diag(2)*1000)))

## Model using prior.parex.5
model.px <- MCMCglmm(PE~Envt-1, random=~us(Envt):Line, rcov=~idh(Envt):units, data=mydata, prior=prior.parex.5, nitt=5050000, thin=500, burnin=5000, verbose=FALSE)

> summary(model.px)

 Iterations = 50001:5049501
 Thinning interval  = 500
 Sample size  = 10000 

 DIC: 1285.159 

 G-structure:  ~us(Envt):Line

         post.mean   l-95% CI u-95% CI eff.samp
H:H.Line     56.43  1.114e-06    182.4    10000
L:H.Line     26.43 -3.325e+01    108.0    10000
H:L.Line     26.43 -3.325e+01    108.0    10000
L:L.Line     88.10  1.459e-04    237.7     9586

 R-structure:  ~idh(Envt):units

        post.mean l-95% CI u-95% CI eff.samp
H.units     286.2    192.0    386.2    10000
L.units     237.3    164.1    316.9    11141

 Location effects: PE ~ Envt - 1 

           post.mean l-95% CI u-95% CI eff.samp  pMCMC    
EnvtH     46.33    39.91    52.92     9222 <1e-04 ***
EnvtL     44.15    36.73    51.31    10000 <1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

The means of the variance components look more reasonable from the summary than they did in previous models, but the posterior modes seem a bit strange:

> posterior.mode(model.px$VCV)
   H:H.Line    L:H.Line    H:L.Line    L:L.Line     H.units     L.units 
  0.9244938   1.4572214   1.4572214  31.2085852 257.8094014 239.5498357 

This is also having a large effect on calculating heritability of the character in each environment, as L heritability is around 100 times that of H.

I feel as if the parameter-expanded terms are headed more in the right direction, but I don't really know where to go next, and maybe there’s something obvious I’m missing here. If anyone has any ideas, I’d be extremely grateful! 

If there is anything else I should include in order to get some help then please let me know...

Thanks,

Tom


-- 
The University of Stirling is ranked in the top 50 in the world in The Times Higher Education 100 Under 50 table, which ranks the world's best 100 universities under 50 years old.
The University of Stirling is a charity registered in Scotland, 
 number SC 011159.