[R-sig-ME] MCMCglmm bivariate model
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sun Sep 27 11:28:14 CEST 2015
Hi,
1) It is hard to offer advice on this. It seems from the widths of the
credible intervals that the genetic term and the ID term are hard to
separate. What is the structure of the pedigree like? The scaled
F-distribution with large scale has a long tail. My guess is that this
results in a long tail for the animal term which then forces the ID
term to zero, resulting in a bimodal posterior distribution. Although
a scale of 1000 makes sense for trait 1, it is less reasonable for the
threshold trait. Pierre de Villemereuil makes this point in his
Methods in Ecology and Evolution paper (4 3 260–275 2013). He suggests
a prior that approaches the chi-square distribution (V=1, nu=1000,
alpha.mu=0, alpha.V=1) but you can't use this in a bivariate setting
because both traits have to have a common nu. You could try reducing
the scale of the threshold trait to something more reasonable. For
example,
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(c(1000,100))),G2 =
list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(c(1000,100)))))
or perhaps even less (e.g. 10). Why bimodality increases under the
bivariate specification I'm not sure. I have not seen any papers
exploring the bivariate scaled-F prior theoretically. It would be nice
if some one did! The new version of MCMCglmm (now on-line) can fit
antedependence structures. Here, priors can be placed directly on
regression coefficients (e.g. the breeding value for y2 on y1) which
may turn out to have nicer properties - not sure.
2) There is nothing to worry about, M3 is the correct model and it is
just a reparameterisation of M1. From the bivariate model location 3
trait 1 is the intercept (457), location 1 trait 1 is (457+11 = 468)
and location 2 trait 1 is (457+38 = 495). From the univariate model
location 1 is the intercept (469) location 2 is (469+23 = 491) and
location 3 is (469-14 = 455). They differ slightly of course, but
there is Monte Carlo error, and the model is slightly different (e.g.
there are correlations between traits).
Cheers,
Jarrod
2015 15:47:36 +0000:
> 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]]
>
> _______________________________________________
> 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.
More information about the R-sig-mixed-models
mailing list