[R-sig-ME] MCMCglmm bivariate model

Sardell, Rebecca Joanne r.sardell at med.miami.edu
Thu Sep 24 17:47:36 CEST 2015


Hi all,
I'm using MCMCglmm to estimate the heritability of 2 traits and their genetic correlation in a repeated measures animal model. My data consists of 500 individuals each measured twice for each trait with a few missing values. y1 is Gaussian, y2 is ordered discrete categories 0-5. Separate models for each trait run ok but I have 2 issues with a bivariate model.


1)      the bivariate model gives a bimodal posterior distribution for the trait y2:ID term with one peak at ~4 (similar to the univariate model) and a second peak at zero. I'm using parameter expanded priors and both mixing and convergence look ok according to trace plots, Geweke plots and Heidelberg tests. Autocorrelation is <0.13. The ID term does have some support at zero in a univariate model of y2 but has a clear mode at 4, whereas in the bivariate model there is equal support for 0 and 4.



2)      The fixed effect, location (a 3 level factor), is only relevant for trait y1. In the univariate model for y1 the 3 levels are estimated ok with no error but in the bivariate model using trait-1 + at.level(trait, 1):location gives the error "In MCMCglmm(cbind(y1, y2) ~ trait - 1 + at.level(trait,  :  some fixed effects are not estimable and have been removed. Use singular.ok=TRUE to sample these effects, but use an informative prior!" and gives estimates in the summary table for locations 1 and 2, suggesting location 3 could not be estimated? Using trait-1+trait:location does not give an error and gives estimates for locations 2 and 3 in the summary table for both traits, which presumably are contrasts from the "trait y1" that corresponds to location 1.


Can anyone advise on these 2 issues? Am I using at.level(trait, 1) correctly? I've copied my code and output below:
Thanks in advance,
Rebecca
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
univariate models:
P1<- list(G = list(G1 = list(V = 1, n = 0.002), G2 = list(V = 1,n = 0.002)), R = list(V = 1, n = 0.002))
M1<-MCMCglmm(y1~age+location+size,
random=~animal+ID,
pedigree=ped,
data=pheno2,prior=P1,nitt=1010000,thin=1000,burnin=10000,verbose=FALSE)

> summary(M1)
Iterations = 10001:1009001
Thinning interval  = 1000
Sample size  = 1000

DIC: 10909.28

G-structure:  ~animal

       post.mean l-95% CI u-95% CI eff.samp
animal      1144    444.3     1885     1162

               ~ID

   post.mean l-95% CI u-95% CI eff.samp
ID      1750     1184     2462     1000

R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units     793.1    700.6    896.2     1000

Location effects: y1 ~ age + location + size

            post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)   469.474  434.740  502.941     1000 <0.001 ***
age            -3.409   -3.950   -2.958     1000 <0.001 ***
location2    23.355   10.235   37.260     1000  0.002 **
location3   -14.113  -29.935    3.468     1000  0.114
size      4.782    2.300    7.316     1000 <0.001 ***
---

> posterior.heritability<-M1$VCV[,"animal"]/(M1$VCV[,"animal"]+M1$VCV[,"ID"]+M1$VCV[,"units"])
> mean(posterior.heritability)
[1] 0.31

~~~~~~~~~~~~~~~~~~~~~~~~~~
P2 <- list(R=list(V = 1 , fix=1), G = list(G1=list(V=1, nu=1,alpha.mu=0, alpha.V=1000),G2=list(V=1, nu=1, alpha.mu=0, alpha.V=1000)))
M2<-MCMCglmm(y2~age,random=~animal+ID,pedigree=ped,family="threshold",data=pheno2,prior=P2,verbose=FALSE, nitt=10500000,thin=2000,burnin=500000)
> summary(M2)

Iterations = 500001:10498001
Thinning interval  = 2000
Sample size  = 5000

DIC: 1270.806

G-structure:  ~animal

       post.mean l-95% CI u-95% CI eff.samp
animal     10.38    4.473    16.46     4162

               ~ID

   post.mean l-95% CI u-95% CI eff.samp
ID     4.291 2.09e-05     8.01     3910

R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

Location effects: y2 ~ age

            post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)  -11.6004 -14.5259  -8.9838     5000 <2e-04 ***
age            0.1794   0.1404   0.2218     5000 <2e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cutpoints:
                      post.mean l-95% CI u-95% CI eff.samp
cutpoint.traity2.1     2.895    2.461    3.314     5351
cutpoint.trait y2.2     4.610    4.069    5.268     4793
cutpoint.trait y2.3     7.398    6.527    8.348     5000
cutpoint.trait y2.4    11.917   10.013   13.981     5141

> posterior.heritability<-M2$VCV[,"animal"]/(M2$VCV[,"animal"]+M2$VCV[,"ID"]+M2$VCV[,"units"])
> mean(posterior.heritability)
[1] 0.66


~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Bivariate model:
P3 = list(R = list(V = diag(2), nu = 1.002, fix = 2), G = list(G1 = list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000),G2 = list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000)))
M3<-MCMCglmm(cbind(y1, y2) ~ trait - 1 + at.level(trait, 1):location + trait:age + at.level(trait, 1):size,
pedigree=ped, data=pheno2,
random = ~us(trait):animal+us(trait):ID,
rcov = ~us(trait):units,  family = c("gaussian", "threshold"),
prior = P3, verbose = FALSE, nitt=10500000,thin=2000,burnin=500000)

> summary(M3)

Iterations = 500001:10498001
Thinning interval  = 2000
Sample size  = 5000

DIC: 12151.2

G-structure:  ~us(trait):animal

                           post.mean l-95% CI u-95% CI eff.samp
y1:y1.animal   1156.88  435.518 1870.316     5000
y2:y1.animal        -39.70  -83.129    5.598     4666
y1:y2.animal        -39.70  -83.129    5.598     4666
y2:y2.animal             11.04    4.659   17.753     3657

               ~us(trait):ID

                       post.mean   l-95% CI u-95% CI eff.samp
y1:y1.ID  1743.758  1.124e+03 2351.973     5000
y2:y1.ID       -12.174 -4.464e+01   20.891     5000
y1:y2.ID       -12.174 -4.464e+01   20.891     5000
y2:y2.ID             3.876  3.188e-06    7.734     3514

R-structure:  ~us(trait):units

                          post.mean l-95% CI u-95% CI eff.samp
y1:y1.units   791.876  696.035    882.1     4568
y2:y1.units         2.582   -1.169      6.1     5000
y1:y2.units         2.582   -1.169      6.1     5000
y2:y2.units             1.000    1.000      1.0        0

Location effects: cbind(y1, y2) ~ trait - 1 + at.level(trait, 1):loca

                               post.mean l-95% CI u-95% CI eff.samp  pMCMC
traity1                  456.6321 419.7208 494.7115     5239 <2e-04 ***
traity2                      -11.5170 -14.4037  -8.8785     4137 <2e-04 ***
at.level(trait, 1):location1   11.6986  -4.0319  28.8620     5000 0.1676
at.level(trait, 1):location2   37.8950  22.5601  53.8510     5000 <2e-04 ***
traity1:age               -3.4116  -3.9302  -2.9405     5209 <2e-04 ***
traity2:age                    0.1779   0.1383   0.2184     4670 <2e-04 ***
at.level(trait, 1):size     4.4378   1.9564   6.9510     5421 0.0004 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cutpoints:
                      post.mean l-95% CI u-95% CI eff.samp
cutpoint.traity2.1     2.890    2.473    3.308     4696
cutpoint.trait y2.2     4.600    4.074    5.186     4474
cutpoint.trait y2.3     7.451    6.565    8.327     4071
cutpoint.trait y2.4    12.053   10.131   14.039     4045


herit_y1<- M3$VCV[,'y1:y1.animal']/( M3$VCV[,'y1:y1.animal']+ M3$VCV[,'y1:y1.ID']+ M3$VCV[,'y1:y1.units'])
mean(herit_y1)
=0.31 (0.13-0.48)

herit_y2<- M3$VCV[,'y2:y2.animal']/( M3$VCV[,'y2:y2.animal']+ M3$VCV[,'y2:y2.ID']+ M3$VCV[,'y2:y2.units'])
mean(herit_y2)
=0.68 (0.42-0.95)

corr.gen<- M3$VCV[,'y1:y2.animal']/sqrt(M3$VCV[,'y1:y1.animal']* M3$VCV[,y2:y2.animal'])
mean(corr.gen)
=-0.36 (-0.70-0.004)


~~~~~~~~~~~~~~~~~~~~~~~~~~~~
summary table for bivariate model with trait:location instead of at.level(trait, 1):location

Location effects: cbind(y1, y2) ~ trait - 1 + trait:location + trait:age + at.level(trait, 1):size

                              post.mean l-95% CI u-95% CI eff.samp  pMCMC
traity1                 468.8456 433.1704 499.3057   761.34 <0.001 ***
traity2                     -11.6934 -14.6120  -9.0283    88.10 <0.001 ***
traity1:location2      23.5282  10.8358  36.8648  1000.00 <0.001 ***
traity2:location2           0.9043  -0.2198   2.1732   341.51  0.128
traity1:location3     -14.6928 -30.8528   1.0206   887.36  0.074 .
traity2:location3           0.8966  -0.4732   2.3384   453.18  0.218
traity1:age              -3.3947  -3.9142  -2.9634   839.89 <0.001 ***
traity2:age                   0.1718   0.1356   0.2125    92.35 <0.001 ***
at.level(trait, 1):size    4.4933   2.2326   7.3976  1000.00  0.002 **



	[[alternative HTML version deleted]]



More information about the R-sig-mixed-models mailing list