[R-sig-ME] Binary response animal models in MCMCglmm
Pierre de Villemereuil
pierre.de.villemereuil at mailoo.org
Mon Mar 27 22:54:23 CEST 2017
Hi Jacob, hi Jarrod,
A bit ago, I wrote a script for a friend to view the priors for multi-response models when one response is binary (hence with fixed residual variance). I decided to share it on GitHub, see here:
https://github.com/devillemereuil/prior-MCMCglmm
It will help you Jacob to see what Jarrod is meaning about priors I think. It might be useful to even more people.
Cheers,
Pierre.
On Monday, 27 March 2017 20:49:39 NZDT Jarrod Hadfield wrote:
> Hi Jacob,
>
> 1/ I usually go for the posterior mode for variances in which the
> posteriors are strongly skewed.
>
> 2/ Compared to family="ordinal", the DIC in threshold models is a level
> closer to what most people are after. But its still a level to high - it
> is asking how well could we predict new observations associated with the
> particular random effects (breeding values in this case) we've happened
> to sample. I don't use DIC.
>
> 3/ Because you have V=diag(2) *and* fix=1 you are fixing the residual
> covariance matrix to an identity matrix. This might be OK if your two
> responses are a male and female trait and so each row only has one value
> and one NA? Even then this makes the algorithm pretty inefficient:
> better to have a single trait indexed by sex. The intersexual genetic
> correlation is then estimated using us(sex):animal.
>
> If this isn't the case then fixing the residual covariance to zero will
> force some of the residual correlation into the genetic term. It is
> possible then to obtain an *estimate* of zero for the genetic
> correlation if the *true* genetic correlation is positive and the
> residual covariance negative.
>
> Much more likely though is that the prior you use is strongly
> constraining the genetic correlation to zero. I guess your motivation
> for using it is that the marginal priors for each variance have nu=1000,
> and this is what Pierre de Villemereuil suggests for probit models in
> his MEE paper (i.e. tending towards a chi-square with 1df)? However, for
> bivariate models the marginal prior for the correlation is flat when
> nu=3 (and pushed outwards towards -1/1 when nu=2). With nu=1000 it is
> pushed inwards towards zero by quite a bit.
>
> 4/ The scale in binary models is not identifiable and so absolute values
> of the genetic variance do not have a meaning outside of the context of
> the residual variance. I've tried to explain this in pictures in the
> supp materials to this paper:
>
> http://onlinelibrary.wiley.com/doi/10.1111/2041-210X.12354/abstract
>
> but not sure if it helps. The other way I try and explain it is that the
> residual variance is not identifiable (Section 2.6 of the course notes).
> If there are better ways of explaining it, it would be good to know!
>
> Cheers,
>
> Jarrod
>
>
> On 27/03/2017 05:18, Jacob Berson wrote:
> > Hi All
> >
> >
> >
> > I have attempted to estimate the heritability of several binary traits, and
> > test for genetic correlations both between two binary traits, as well as
> > between a binary and a Gaussian trait, using the animal model in MCMCglmm.
> >
> >
> >
> > Several issues have come up in my analysis that I'm hoping to get some help
> > with.
> >
> >
> >
> > For some traits the posterior distribution of my heritability estimate is
> > very close to zero, resulting in a non-symmetric density plot (the left tail
> > is essentially cut by the y-axis). With these distributions I get very
> > different point estimates of heritability depending on if I use the
> > posterior mean or posterior mode.
> >
> >
> >
> > 1) I've seen both the mean and mode used in the literature, does anyone
> > have a view on which is the most appropriate in the above circumstance?
> >
> >
> >
> > The abovementioned distributions suggest to me that there is little, if any,
> > additive genetic variance present. However, for some traits the DIC value of
> > a model with "animal" as a random effect is lower than the null model
> > (though after reading
> > https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q1/021875.html and
> > switching from using family="ordinal" to family="threshold" the difference
> > in DIC values was greatly reduced).
> >
> >
> >
> > 2) Are differences in DIC values an appropriate tool for testing model
> > fit when the response is binary? If so, what level of difference is
> > sufficient to reject the null model (I've seen both >5 and >10)?
> >
> >
> >
> > I am getting a posterior distribution centred on zero when testing for
> > intersexual genetic correlations between two binary traits. However, when I
> > look at the data (for example using sire means) it seems to me that there is
> > in fact a strong correlation and that my MCMCglmm results are not correct.
> >
> >
> >
> > 3) Am I trying to do the impossible, i.e., is it still recommended (as
> > was the case a few years ago
> > https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q3/002637.html) not to
> > fit bivariate binary models in MCMCglmm, even when using family="threshold"?
> >
> >
> >
> > Finally, as I understand it, it is not appropriate to report the estimates
> > of the additive genetic variance from binary models because it depends on
> > the residual variance (which I have fixed to 1).
> >
> >
> >
> > 4) Can anyone help a novice like me understand why heritability is more
> > appropriate to report, even though it also incorporates the fixed residual
> > variance?
> >
> >
> >
> > Apologies for the length of this post - after much searching and reading I
> > haven't been able to find solutions to these questions. Any advice on the
> > above and/or feedback on my code below would be very much appreciated.
> >
> >
> >
> > Jacob
> >
> >
> >
> >
> >
> > My code:
> >
> >
> >
> > #Univariate models for binary traits
> >
> >
> >
> > Prior0 <- list(R = list(V = 1, fix = 1))
> >
> >
> >
> > Prior1 <- list(R = list(V = 1, fix = 1),
> >
> > G = list(G1 = list(V = 1, nu = 1000, alpha.mu = 0, alpha.V =
> > 1)))
> >
> >
> >
> > Binary.Null <- MCMCglmm(bin.response ~ 1, family = "threshold", pedigree =
> > Ped,
> >
> > prior = Prior0, data = Data, nitt = 1050000, burnin = 50000,
> > thin = 500, verbose = FALSE)
> >
> >
> >
> > Binary.Va <- MCMCglmm(bin.response ~ 1, random = ~animal, family =
> > "threshold", pedigree = Ped,
> >
> > prior = Prior1, data = Data, nitt = 1050000, burnin = 50000, thin = 500,
> > verbose = FALSE)
> >
> >
> >
> > heritability <- Binary.Va $VCV[,"animal"] / rowSums(Binary.Va [["VCV"]])
> >
> >
> >
> >
> >
> > #Genetic correlation between two binary traits (bin1 and bin2)
> >
> >
> >
> > Prior1rG <- list(R = list(V = diag(2), nu = 0, fix = 1),
> >
> > G = list(G1 = list(V = diag(2), nu = 10001, alpha.mu =
> > c(0,0), alpha.V = diag(2))))
> >
> >
> >
> > Bin.Bin.corr <- MCMCglmm(cbind(bin1, bin2) ~ trait - 1, random =
> > ~us(trait):animal,
> >
> > rcov = ~corg(trait):units, family = c("threshold", "threshold"),
> > pedigree = Ped, prior = Prior1rG,
> >
> > data = Data, nitt = 4050000, burnin = 50000, thin = 2000, verbose =
> > FALSE)
> >
> >
> >
> > Genetic_correlation1 <- Bin.Bin.corr $VCV[, "traitbin1:traitbin2"] /
> > sqrt(Bin.Bin.corr$VCV[,
> >
> > "traitbin1:trait bin1.animal"] * Bin.Bin.corr $VCV[,
> > "traitbin2:traitbin2.animal"])
> >
> >
> >
> >
> >
> > #Genetic correlations between a Gaussian (gau1) and a binary (bin2) trait
> >
> >
> >
> > Prior2rG <- list(R = list(V = diag(2), nu = 0, fix = 2),
> >
> > G = list(G1 = list(V = diag(2), nu = 2, alpha.mu = c(0,0),
> > alpha.V = diag(2)*1000)))
> >
> >
> >
> > Gau.Bin.corr <- MCMCglmm(cbind(gau1, bin2) ~ trait - 1, random =
> > ~us(trait):animal, rcov =
> >
> > ~us(trait):units, family = c("gaussian", "threshold"), pedigree = Ped,
> >
> >
> > prior = Prior2rG, data = Data, nitt = 1050000, burnin = 50000, thin = 500,
> >
> > verbose = FALSE)
> >
> >
> >
> > Genetic_correlation2 <- Gau.Bin.corr$VCV[, "traitgau1:traitbin2"] /
> > sqrt(Gau.Bin.corr $VCV[, "trait
> >
> > gau1:trait gau1.animal"] * Gau.Bin.corr $VCV[,
> > "traitbin2:traitbin2.animal"])
> >
> >
> >
> >
> >
> >
> >
> > Jacob Berson BSc (Hons), PhD Candidate
> >
> > Centre for Evolutionary Biology
> >
> >
> >
> > School of Animal Biology (M092)
> >
> > University of Western Australia
> >
> > 35 Stirling Highway
> >
> > Crawley, WA, 6009
> >
> >
> >
> > Tel: (+61 8) 6488 3425
> >
> >
> >
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
>
More information about the R-sig-mixed-models
mailing list