[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