[R-sig-ME] MCMCglmm: square root of the sampling variance of additive genetic variance

[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)

summary(model)

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

HPDinterval(model$VCV)

lower upper

animal 38.18195 120.3350

units 59.50312 109.1574

attr(,"Probability")

[1] 0.95

herit <- model$VCV[, "animal"]/(model$VCV[, "animal"] + model$VCV[,

"units"])

mean(herit)

[1] 0.4772017

HPDinterval(herit, 0.95)

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

sd(model$VCV[, "animal"])

[1] 21.36365

sd(model$VCV[, "units"])

[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

---

[[alternative HTML version deleted]]

[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
--
Walid Mawass
Ph.D. candidate in Cellular and Molecular 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: [phone] poste 3384

Post by Simona Kralj Fiser
Hi.
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)

summary(model)

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

HPDinterval(model$VCV)

lower upper
animal 38.18195 120.3350
units 59.50312 109.1574
attr(,"Probability")
[1] 0.95

herit <- model$VCV[, "animal"]/(model$VCV[, "animal"] + model$VCV[,

"units"])

mean(herit)

[1] 0.4772017

HPDinterval(herit, 0.95)

lower upper
var1 0.2940021 0.6576105
attr(,"Probability")
[1] 0.95
1. I am trying to call the standard errors for additive and residual
variance
se(model$VCV[, "animal"]) does not work
I used

sd(model$VCV[, "animal"])

[1] 21.36365

sd(model$VCV[, "units"])

[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
---
[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]

[Pierre de Villemereuil]

Hi,

Post by Walid Mawass
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.

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.

[Simona Kralj Fiser]

Thank you very much. That helps a lot!
Simona

Post by 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
--
Walid Mawass
Ph.D. candidate in Cellular and Molecular 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: [phone] poste 3384

Post by Simona Kralj Fiser
Hi.

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

Post by Simona Kralj Fiser
model <- MCMCglmm(trait~ 1, random = ~animal, pedigree = pedigree,data =
data, nitt = 5000000, thin = 100, burnin = 150000, prior = prior,

verbose =

Post by Simona Kralj Fiser
FALSE)

summary(model)

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

HPDinterval(model$VCV)

lower upper
animal 38.18195 120.3350
units 59.50312 109.1574
attr(,"Probability")
[1] 0.95

herit <- model$VCV[, "animal"]/(model$VCV[, "animal"] + model$VCV[,

"units"])

mean(herit)

[1] 0.4772017

HPDinterval(herit, 0.95)

lower upper
var1 0.2940021 0.6576105
attr(,"Probability")
[1] 0.95
1. I am trying to call the standard errors for additive and residual
variance
se(model$VCV[, "animal"]) does not work
I used

sd(model$VCV[, "animal"])

[1] 21.36365

sd(model$VCV[, "units"])

[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
---
[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]