[R-sig-ME] MCMCglmm Different Outputs for ginverse Inputs

[Cameron So]

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]]