Hello MCMCglmm users,
I am estimating the Va and heritability of a trait in separate univariate and bivariate models.
Some background beforehand: In my experiment, individuals are exposed to separate treatments.
The pedigree for each treatment is quite similar given I am using plants + multiple replicates.
Therefore, I have two univariate models: a subset of the dataset for each treatment. This has allowed
me to compare the Va and heritability estimates between treatments using the univariate models.
In the multivariate model, I include the entire dataset so I can evaluate cross-environment COVa.
I am noticing that my heritability estimates between the univariate and bivariate models differ
substantially. Could anyone offer any advice on what I am observing?
NOTE: I am not having issues of autocorrelation, nor convergence. Additionally, the univariate
Models include a second random effect term (matID) but this shouldn?t affect the additive (animal)
term that drastically. The priors obviously do differ between the univariate and bivariate models
because of the binary character of ?treatment?.
Below are my models:
## UNIVARITE MODEL EXAMPLE ##
prior8.3b <- list(R = list(V = 1, nu = 0.002),
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)))
HW_model8.3b <- MCMCglmm(leaf ~ plot, random = ~animal + matID,
ginverse = list(animal = Ainv),
family = "gaussian", data = heated.leaf, prior = prior8.3b,
nitt = 2100000, thin = 1000, burnin = 100000, verbose = T, pr = TRUE)
HW_herit8.3b <- HW_model8.3b$VCV[, "animal"]/(HW_model8.3b$VCV[, "animal"] + HW_model8.3b$VCV[, "matID"] + HW_model8.3b$VCV[, "units"])
mean(HW_herit8.3b)
#h2 = ~ 0.176
## BIVARIATE MODEL EXAMPLE ##
prior6.2 <- list(R = list(V = diag(1), fix = 1),
G = list(G1 = list(V = diag(2), nu = 2, alpha.mu = c(0,0), alpha.V = diag(c(1,1)))))
PL_model6.2 <- MCMCglmm(leaf ~ plot:treatment + treatment, random = ~us(treatment):animal,
ginverse = list(animal = Ainv), rcov=~units,
family = "gaussian", data = plastic.leaf, prior = prior6.2,
nitt = 4100000, thin = 2000, burnin = 100000, verbose = T, pr = TRUE, trunc = TRUE)
herit_PL9.2_A <-PL_model9.2$VCV[,'treatmentA:treatmentA.animal']/ (PL_model9.2$VCV[,'treatmentA:treatmentA.animal'] + 1)
mean(herit_PL9.2_A)
#h2 = ~0.73
Thank you for any help in advance!
Cameron
MCMCglmm Contrasting h2 in univariate and bivariate model
2 messages · Cameron So, Walid Mawass
Hello Cameron, Based on what you have sent, I can comment on a few things that might help. 1- the second model you are using is not a bivariate model since you only have a single response variable. What you have constructed is a univariate model with a random interaction, in your case its the interaction of the random term animal with the treatment variable which I presume has only two levels since you specified a 2x2 matrix for the random term. 2- You have different priors between the two models since you fixed the residual term to 1 which is not needed here since you are modeling the same response variable leaf with a Gaussian distribution. This could explain the differences you got in the heritability estimates. 3- is there a reasoning for using trunc=TRUE in your second model? hope this helps to clear things up for you and good luck
Walid Mawass Ph.D. student in Evolutionary Biology Population Genetics Laboratory University of Qu?bec at Trois-Rivi?res 3351, boul. des Forges, C.P. 500 Trois-Rivi?res (Qu?bec) G9A 5H7 Telephone: 819-376-5011 poste 3384 On Wed, Apr 29, 2020 at 3:00 PM Cameron So <cameron.so at mail.utoronto.ca> wrote: > Hello MCMCglmm users, > > I am estimating the Va and heritability of a trait in separate univariate > and bivariate models. > Some background beforehand: In my experiment, individuals are exposed to > separate treatments. > The pedigree for each treatment is quite similar given I am using plants + > multiple replicates. > Therefore, I have two univariate models: a subset of the dataset for each > treatment. This has allowed > me to compare the Va and heritability estimates between treatments using > the univariate models. > In the multivariate model, I include the entire dataset so I can evaluate > cross-environment COVa. > > I am noticing that my heritability estimates between the univariate and > bivariate models differ > substantially. Could anyone offer any advice on what I am observing? > > NOTE: I am not having issues of autocorrelation, nor convergence. > Additionally, the univariate > Models include a second random effect term (matID) but this shouldn?t > affect the additive (animal) > term that drastically. The priors obviously do differ between the > univariate and bivariate models > because of the binary character of ?treatment?. > > Below are my models: > > ## UNIVARITE MODEL EXAMPLE ## > > prior8.3b <- list(R = list(V = 1, nu = 0.002), > 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))) > > HW_model8.3b <- MCMCglmm(leaf ~ plot, random = ~animal + matID, > ginverse = list(animal = Ainv), > family = "gaussian", data = heated.leaf, prior = > prior8.3b, > nitt = 2100000, thin = 1000, burnin = 100000, verbose > = T, pr = TRUE) > > HW_herit8.3b <- HW_model8.3b$VCV[, "animal"]/(HW_model8.3b$VCV[, "animal"] > + HW_model8.3b$VCV[, "matID"] + HW_model8.3b$VCV[, "units"]) > mean(HW_herit8.3b) > #h2 = ~ 0.176 > > ## BIVARIATE MODEL EXAMPLE ## > > > prior6.2 <- list(R = list(V = diag(1), fix = 1), > G = list(G1 = list(V = diag(2), nu = 2, alpha.mu = > c(0,0), alpha.V = diag(c(1,1))))) > > PL_model6.2 <- MCMCglmm(leaf ~ plot:treatment + treatment, random = > ~us(treatment):animal, > ginverse = list(animal = Ainv), rcov=~units, > family = "gaussian", data = plastic.leaf, prior = > prior6.2, > nitt = 4100000, thin = 2000, burnin = 100000, > verbose = T, pr = TRUE, trunc = TRUE) > > herit_PL9.2_A <-PL_model9.2$VCV[,'treatmentA:treatmentA.animal']/ > (PL_model9.2$VCV[,'treatmentA:treatmentA.animal'] + 1) > mean(herit_PL9.2_A) > #h2 = ~0.73 > > > Thank you for any help in advance! > > Cameron > > > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-mixed-models at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models > [[alternative HTML version deleted]]