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

Simona Kralj Fiser @imon@kf @ending from zrc-@@zu@@i
Tue Nov 13 11:39:55 CET 2018


Thank you very much. That helps a lot!
Simona

On Wed, 7 Nov 2018 at 18:33, Walid Mawass <walidmawass10 using gmail.com> wrote:

> 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: 819-376-5011 poste 3384
>
>
> On Wed, Nov 7, 2018 at 12:09 PM Simona Kralj Fiser <simonakf using zrc-sazu.si>
> wrote:
>
> > 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]]
> >
> > _______________________________________________
> > R-sig-mixed-models using r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
>         [[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