[R-sig-ME] Binary response animal models in MCMCglmm
Jarrod Hadfield
j.hadfield at ed.ac.uk
Mon Mar 27 21:49:39 CEST 2017
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
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list