[R-sig-ME] Binary response animal models in MCMCglmm
Jarrod Hadfield
j.hadfield at ed.ac.uk
Tue Mar 28 10:30:43 CEST 2017
Hi Jacob,
Using rcov=~corg(trait):units constrains the residual variances to one
but estimates the residual correlation. However, you shouldn't use fix
in the prior.
Cheers,
Jarrod
On 28/03/2017 09:24, Jacob Berson wrote:
> Hi Jarrod and Pierre
>
> Thank you both very much for your responses.
>
> Pierre, that script is really useful - thanks for making it available.
>
> Jarrod, your answers have been a great help. In regards to 3/, yes you are
> correct that my two responses were male and female traits so each row only
> had one value and one NA. And yes, my motivation for using that prior for
> this model was an attempt at a bivariate form of the prior suggested in de
> Villemereuil et al 2013. Both running the model with nu=3, and a different
> model using us(sex):animal, have given me similar estimates which are much
> closer to what I would expect.
>
> As a follow up...
>
> In the case of a genetic correlation between two binary traits where both
> traits have been measured on the same individual (i.e. each row has a value
> for both responses), is there a way of fixing both of the residual variances
> but estimating the residual correlation? Am I right in my understanding of
> below that using V=diag(2) and fix=1 should always be avoided in this
> circumstance as it will bias the genetic correlation estimate in the
> direction of the residual correlation?
>
> Thanks again for the help, I really appreciate it.
>
> Jacob
>
> -----Original Message-----
> From: Pierre de Villemereuil [mailto:pierre.de.villemereuil at mailoo.org]
> Sent: Tuesday, 28 March 2017 4:54 AM
> To: r-sig-mixed-models at r-project.org; Jacob Berson
> <jacob.berson at research.uwa.edu.au>
> Subject: Re: [R-sig-ME] Binary response animal models in MCMCglmm
>
> 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
>>
>>
>
--
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