Skip to content

MCMCglmm: square root of the sampling variance of additive genetic variance

4 messages · Simona Kralj Fiser, Walid Mawass, Pierre de Villemereuil

Original
Simona Kralj Fiser
# 
Hi.
I used MCMCglmm to calculate the heritability of a trait [Va/(Va+Vr)]; e.g.:

prior<-list(G=list(G1=list(V=matrix(p.var*0.5),n=1)),R=list(V=matrix(p.var*0.5),n=1)

model <- MCMCglmm(trait~ 1, random = ~animal, pedigree = pedigree,data =
data, nitt = 5000000, thin = 100, burnin = 150000, prior = prior, verbose =
FALSE)
Iterations = 150001:4999901

 Thinning interval  = 100

 Sample size  = 48500



 DIC: 2032.226



 G-structure:  ~animal



       post.mean l-95% CI u-95% CI eff.samp

animal     78.48    38.18    120.3    48500



 R-structure:  ~units



      post.mean l-95% CI u-95% CI eff.samp

units     84.11     59.5    109.2    48500



 Location effects: trait~ 1



            post.mean l-95% CI u-95% CI eff.samp  pMCMC

(Intercept)     6.918    4.589    9.091    48500 <2e-05 ***
lower    upper

animal 38.18195 120.3350

units  59.50312 109.1574

attr(,"Probability")

[1] 0.95
"units"])
[1] 0.4772017
lower     upper

var1 0.2940021 0.6576105

attr(,"Probability")

[1] 0.95


I have two questions:

1. I am trying to call the standard errors for additive and residual
variance


se(model$VCV[, "animal"]) does not work


I used
[1] 21.36365
[1] 12.8011



I wonder whether the SD (of Va) provides the square root of the sampling
variance of Va. Could you please confirm this? I am interested in
calculating the SE of Va to calculate the SEs of other statistics (e.g.,
CVa).

2. Also, is there a way to plot the posterior distribution of the
heritability (or Va) estimates?

Thank you!

Simona

---
Walid Mawass
# 
Hello,

The MCMCglmm estimates the posterior distribution of the additive and
residual variances, to my knowledge then, there is no standard error
associated with it rather you can output the HPD interval or highest
posterior density interval at a 95 or 98% confidence interval, which you
already have done. I stand to be corrected though.

for your second question, it is quite easy, you just need use the basic
plot function. In your case, plot(model$VCV[, "animal"]), this will return
a plot of the posterior distribution and the iterations of the Markov Chain
to visually check for autocorrelation between iterations.(there are other
packages to plot mcmc outputs if you dont want to go with the basic plot,
bayesplot comes to mind)

Good luck
Pierre de Villemereuil
# 
Hi,

Just to clarify this:
The sd() of the MCMC chain is indeed an estimate of the standard error (+ some error itself coming from the Monte Carlo process), though HPD interval is indeed probably a better (and more Bayesian, in a way) way to measure the uncertainty. 

I might also stand to be corrected! ;)

Cheers,
Pierre.
5 days later
Simona Kralj Fiser
# 
Thank you very much. That helps a lot!
Simona
On Wed, 7 Nov 2018 at 18:33, Walid Mawass <walidmawass10 at gmail.com> wrote: