[R-sig-ME] MCMCglmm Different Outputs for ginverse Inputs
Cameron So
c@meron@@o @end|ng |rom m@||@utoronto@c@
Tue Feb 18 23:31:26 CET 2020
Hi everyone,
I am using MCMCglmm to estimate additive genetic variance, dominance, and maternal effects in a plant population when exposed to a new environment and normal conditions (different individuals but replicates of genotype). I have two questions, the first being long and the latter being short.
Q1:
Since I am trying to incorporate a dominance random effect in my model, I initially followed the tutorial in Wolak's (2012) paper [ see: https://doi.org/10.1111/j.2041-210X.2012.00213.x ], in which the inverse additive matrix is indicated (see below). However, after checking model diagnostics and ensuring convergence and no autocorrelation, I was observing very low variance for additive, dominance, and subsequently heritability estimates (<0.1 for h2) for numerous univariate trait models. Obviously, this is not what I expected as heritability for various fitness-related traits is usually between 0.1 - 0.7 according to the literature.
Recently, I ran a new model without indicating the inverse additive matrix and let MCMCglmm automatically create it. The output for this model (although there is some high correlation issues that I am currently fixing) had an heritability estimate ~0.7, greater than the previous model.
Would anyone know the reasoning behind these model output differences when the inverse additive matrix is indicated or not? It appears that excluding the Ainv in the model yields a more accurate estimation.
Below is the an example of the R code for the trait plant height with parameter expanded priors.
Shared Priors:
prior9.2 <- 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),
G3 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000)))
Model 1 Example: Including Ainv
Ainv <- inverseA(ped[, 1:3])$Ainv
Dinv <- makeD(ped, 1:3])$Dinv
HW_model9.2 <- MCMCglmm(height ~ plot, random = ~animal + animalDom + matID,
ginverse = list(animal = Ainv, animalDom = Dinv),
family = "gaussian", data = heated.height, prior = prior9.2,
nitt = 2100000, thin = 1000, burnin = 100000, verbose = T, pr = TRUE)
Model 2 Example: Excluding Ainv
listD <- makeD(ped)
Dinv <- listD$Dinv
HW_model9.2b <- MCMCglmm(height ~ plot, random = ~animal + animalDom + matID,
ginverse = list(animalDom = Dinv),
family = "gaussian", data = heated.height, prior = prior9.2,
nitt = 2100000, thin = 1000, burnin = 100000, verbose = T, pr = TRUE)
..
Q2: Regarding priors, I have a shorter question on the parameterization of V and nu for the Inverse Wishart prior. My understanding is that the Inverse Gamma (0.001, 0.001) prior follows (V=1, nu = 0.002). However, for the Inverse Wishart prior, different sources allude to different parameterizations of V and nu. For example, for a multivariate model, once source says it is V=diag(n), nu = n - 0.998 where n is the number of traits, while another source states that it is V=diag(n), nu = n. Would someone be able to clarify my confusion? My apologies as I am quite new to using MCMCglmm and setting appropriate priors.
______
Cameron So
Master's Student | Plant Evolutionary Responses to Climate Change | Weis Lab
Department of Ecology & Evolutionary Biology
University of Toronto
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list