[R-sig-ME] Unsubscribe mailinglist
Hartgerink, D.M. (Marjolein)
d.hartgerink at student.ru.nl
Tue Mar 28 10:51:04 CEST 2017
Dear sir/madam,
I would like to be unsubscribed from this mailing list.
Thank you in advance,
Marjolein Hartgerink
________________________________________
Van: R-sig-mixed-models [r-sig-mixed-models-bounces at r-project.org] namens r-sig-mixed-models-request at r-project.org [r-sig-mixed-models-request at r-project.org]
Verzonden: dinsdag 28 maart 2017 10:31
Aan: r-sig-mixed-models at r-project.org
Onderwerp: R-sig-mixed-models Digest, Vol 123, Issue 23
Send R-sig-mixed-models mailing list submissions to
r-sig-mixed-models at r-project.org
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
or, via email, send a message with subject or body 'help' to
r-sig-mixed-models-request at r-project.org
You can reach the person managing the list at
r-sig-mixed-models-owner at r-project.org
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-mixed-models digest..."
Today's Topics:
1. Re: Binary response animal models in MCMCglmm
(Pierre de Villemereuil)
2. Re: Binary response animal models in MCMCglmm (Jacob Berson)
3. Re: Binary response animal models in MCMCglmm (Jarrod Hadfield)
----------------------------------------------------------------------
Message: 1
Date: Tue, 28 Mar 2017 09:54:23 +1300
From: Pierre de Villemereuil <pierre.de.villemereuil at mailoo.org>
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
Message-ID: <2128458.8BLEj2G8xN at vercors>
Content-Type: text/plain; charset="us-ascii"
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
>
>
>
------------------------------
Message: 2
Date: Tue, 28 Mar 2017 16:24:07 +0800
From: "Jacob Berson" <jacob.berson at research.uwa.edu.au>
To: "'Jarrod Hadfield'" <j.hadfield at ed.ac.uk>, "'Pierre de
Villemereuil'" <pierre.de.villemereuil at mailoo.org>
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Binary response animal models in MCMCglmm
Message-ID: <003b01d2a79c$ad6d3d80$0847b880$@research.uwa.edu.au>
Content-Type: text/plain; charset="us-ascii"
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
>
>
>
------------------------------
Message: 3
Date: Tue, 28 Mar 2017 09:30:43 +0100
From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
To: Jacob Berson <jacob.berson at research.uwa.edu.au>, "'Pierre de
Villemereuil'" <pierre.de.villemereuil at mailoo.org>
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Binary response animal models in MCMCglmm
Message-ID: <37836e79-59d9-9f86-379b-7f0b0e2d1cd0 at ed.ac.uk>
Content-Type: text/plain; charset="windows-1252"; format=flowed
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.
------------------------------
Subject: Digest Footer
_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
------------------------------
End of R-sig-mixed-models Digest, Vol 123, Issue 23
***************************************************
More information about the R-sig-mixed-models
mailing list