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

Matthew Wolak mew0099 @end|ng |rom @uburn@edu
Wed Feb 19 03:46:50 CET 2020


Hi Cameron,

Q1:
In your second model (`HW_model9.2b`) `animal` does not represent a term for the additive genetic effects for each level of "animal" (well plant, I take it), correlated with one another as described by the relatedness matrix (when you supply `ginverse = list(animal = Ainv...`). Because you dropped the ginverse list entry for animal, the model assumes the effects associated with each "animal" are independent (i.e., a diagonal matrix instead of the relatedness matrix).

Depending on whether you have repeated measures or not, this might be modelling the exact same effects as the residuals.

Perhaps with this second model you were trying to utilize the old (and I think less preferred) method of indicating which term represented the additive genetic variance by using the special `animal` label, but you *also need to give an object to the `pedigree` argument.


Q2:
MCMCglmm is constructing inverse Wishart priors for the variance components. The inverse Gamma is a special case of the inverse Wishart (when `V = 1`) with the shape and scale parameters equal to `nu / 2` (so in your example, `...nu=0.002` yields inverse Gamma shape and scale parameters of 0.001. I think you can think of the inverse Wishart as the multivariate generalization of the inverse Gamma. So different arguments for `nu` are just different inverse Wishart distributions.

However, what you have specified (e.g., because `alpha.V = 1000` along with `V=1, nu=1`) are priors for a parameter expanded model that essentially causes the model to use priors from a different distribution. What you have actually given to your model are marginal priors for the variances that are F-distributions (see Gelman 2006 <-- full reference in either a Hadfield paper or the course notes).

Best of luck, I would be excited to see whether you see dominance-by-environment interactions!

Sincerely,
Matthew

****************************************************
Matthew E. Wolak, Ph.D.
Assistant Professor
Department of Biological Sciences
Auburn University
306 Funchess Hall
Auburn, AL 36849, USA

http://qgevoeco.com
Email: matthew.wolak using auburn.edu
Tel: 334-844-9242


________________________________
From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> on behalf of Cameron So <cameron.so using mail.utoronto.ca>
Sent: Tuesday, February 18, 2020 4:31 PM
To: r-sig-mixed-models using r-project.org <r-sig-mixed-models using r-project.org>
Subject: [R-sig-ME] MCMCglmm Different Outputs for ginverse Inputs

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

_______________________________________________
R-sig-mixed-models using r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

	[[alternative HTML version deleted]]



More information about the R-sig-mixed-models mailing list