Dear colleagues, How to calculate heritability (additive genetic variation scaled by total variation) and evolvability (additive genetic variation scaled by the mean of the trait) from the results of a mixed model where the genetics is modeled as sire effects? We are studying genetics in reef-building corals. Traits are measured in batches of larvae representing full-sib families (i.e., each measurement represents an average over many full-sibs). Traits are either direct measures (for example, color, a continuous variable), or as a response to some treatment or time (for example, size). in the former (simpler) case, we are after scalar sire effects (i.e., sire-specific intercepts). For MCMCglmm, I came up with the following: mc = MCMCglmm ( color ~1, random = ~ sire , family = "gaussian", ... ) heritability = 4 * mc$VCV [ , ' sire ' ] / ( mc$VCV [ , ' sire ' ] + mc$VCV [ , ' units ' ] ) posterior.mode ( heritability ) HPDinterval ( heritability ) evolvability = 4 * mc$VCV [ , ' sire ' ] / mc$Sol [ , ' (Intercept)' ] posterior.mode ( evolvability ) HPDinterval ( evolvability ) For lmer: mr = lmer (color ~ ( 1 | sire ), ... ) var.sire = VarCorr ( mr ) $ sire [ 1 ] var.resid = as.numeric ( attributes ( summary(m2))$REmat [ attributes (summary(m2)) $REmat [,1] == "Residual" , 3 ] ) # this is a really ugly one! I feel there must be a nicer way to extract residual variance... I will be grateful if you suggest something heritability = 4 * var.sire / (var.sire + var.resid) evolvability = 4 * var.sire / fixef (mr) [1] Are these correct? in particular, the factor of 4 in front of the heritability / evolvability formulae (since we are looking only at sire effects in full-sib families)? If we are OK so far, the second question is about non-gaussian models. How these formulae would be modified if the response variable is binomial, in the form of counts of successes and failures ? Seems like for MCMCglmm there will be no change except family = "multinomial2", is that correct? lmer, on the other hand, does not report residual variance with binomial model, apparently because it is fixed. For logit link function it seems to be fixed at 3.29, so the heritability formula will look like this: heritability = 4 * var.sire / (var.sire + 3.29) correct?... I feel this is way more than enough for one post... I guess I will ask about the models with sire-specific slopes later. cheers Misha Matz
animal model: calculating heritability and evolvability from sire effects
5 messages · mikhail matz, Paul Johnson, David Duffy
You can estimate heritability and evolvability more directly via estimating the additive genetic variance from an animal model as laid out here: http://www.wildanimalmodels.org/tiki-download_wiki_attachment.php?attId=4 Best wishes, Paul -----Original Message----- From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of mikhail matz Sent: 21 December 2011 09:29 To: r-sig-mixed-models at r-project.org Cc: Elisha Meyer Subject: [R-sig-ME] animal model: calculating heritability and evolvability from sire effects Dear colleagues, How to calculate heritability (additive genetic variation scaled by total variation) and evolvability (additive genetic variation scaled by the mean of the trait) from the results of a mixed model where the genetics is modeled as sire effects? We are studying genetics in reef-building corals. Traits are measured in batches of larvae representing full-sib families (i.e., each measurement represents an average over many full-sibs). Traits are either direct measures (for example, color, a continuous variable), or as a response to some treatment or time (for example, size). in the former (simpler) case, we are after scalar sire effects (i.e., sire-specific intercepts). For MCMCglmm, I came up with the following: mc = MCMCglmm ( color ~1, random = ~ sire , family = "gaussian", ... ) heritability = 4 * mc$VCV [ , ' sire ' ] / ( mc$VCV [ , ' sire ' ] + mc$VCV [ , ' units ' ] ) posterior.mode ( heritability ) HPDinterval ( heritability ) evolvability = 4 * mc$VCV [ , ' sire ' ] / mc$Sol [ , ' (Intercept)' ] posterior.mode ( evolvability ) HPDinterval ( evolvability ) For lmer: mr = lmer (color ~ ( 1 | sire ), ... ) var.sire = VarCorr ( mr ) $ sire [ 1 ] var.resid = as.numeric ( attributes ( summary(m2))$REmat [ attributes (summary(m2)) $REmat [,1] == "Residual" , 3 ] ) # this is a really ugly one! I feel there must be a nicer way to extract residual variance... I will be grateful if you suggest something heritability = 4 * var.sire / (var.sire + var.resid) evolvability = 4 * var.sire / fixef (mr) [1] Are these correct? in particular, the factor of 4 in front of the heritability / evolvability formulae (since we are looking only at sire effects in full-sib families)? If we are OK so far, the second question is about non-gaussian models. How these formulae would be modified if the response variable is binomial, in the form of counts of successes and failures ? Seems like for MCMCglmm there will be no change except family = "multinomial2", is that correct? lmer, on the other hand, does not report residual variance with binomial model, apparently because it is fixed. For logit link function it seems to be fixed at 3.29, so the heritability formula will look like this: heritability = 4 * var.sire / (var.sire + 3.29) correct?... I feel this is way more than enough for one post... I guess I will ask about the models with sire-specific slopes later. cheers Misha Matz _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models The University of Glasgow, charity number SC004401
6 days later
Hello Paul - the problem is, I am not sure that the built-in "animal" modeling in MCMCglmm fits my case. I am not measuring values of my traits in individual animals. Instead, I am measuring averages for full-sib families (coral larvae are quite small, so I am measuring them in batches). Even though I can build a pedigree for my families, denoting each full-sib family as an "animal", this would not be correct since, for example, "half-sibs" (in my case, crosses sharing one parent) would share 1/2 of the genes, not 1/4 as would be the case for individual half-sib measurements. Does it make sense?... cheers Misha
On Dec 21, 2011, at 2:10 PM, Paul Johnson wrote:
You can estimate heritability and evolvability more directly via estimating the additive genetic variance from an animal model as laid out here: http://www.wildanimalmodels.org/tiki-download_wiki_attachment.php?attId=4 Best wishes, Paul -----Original Message----- From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of mikhail matz Sent: 21 December 2011 09:29 To: r-sig-mixed-models at r-project.org Cc: Elisha Meyer Subject: [R-sig-ME] animal model: calculating heritability and evolvability from sire effects Dear colleagues, How to calculate heritability (additive genetic variation scaled by total variation) and evolvability (additive genetic variation scaled by the mean of the trait) from the results of a mixed model where the genetics is modeled as sire effects? We are studying genetics in reef-building corals. Traits are measured in batches of larvae representing full-sib families (i.e., each measurement represents an average over many full-sibs). Traits are either direct measures (for example, color, a continuous variable), or as a response to some treatment or time (for example, size). in the former (simpler) case, we are after scalar sire effects (i.e., sire-specific intercepts). For MCMCglmm, I came up with the following: mc = MCMCglmm ( color ~1, random = ~ sire , family = "gaussian", ... ) heritability = 4 * mc$VCV [ , ' sire ' ] / ( mc$VCV [ , ' sire ' ] + mc$VCV [ , ' units ' ] ) posterior.mode ( heritability ) HPDinterval ( heritability ) evolvability = 4 * mc$VCV [ , ' sire ' ] / mc$Sol [ , ' (Intercept)' ] posterior.mode ( evolvability ) HPDinterval ( evolvability ) For lmer: mr = lmer (color ~ ( 1 | sire ), ... ) var.sire = VarCorr ( mr ) $ sire [ 1 ] var.resid = as.numeric ( attributes ( summary(m2))$REmat [ attributes (summary(m2)) $REmat [,1] == "Residual" , 3 ] ) # this is a really ugly one! I feel there must be a nicer way to extract residual variance... I will be grateful if you suggest something heritability = 4 * var.sire / (var.sire + var.resid) evolvability = 4 * var.sire / fixef (mr) [1] Are these correct? in particular, the factor of 4 in front of the heritability / evolvability formulae (since we are looking only at sire effects in full-sib families)? If we are OK so far, the second question is about non-gaussian models. How these formulae would be modified if the response variable is binomial, in the form of counts of successes and failures ? Seems like for MCMCglmm there will be no change except family = "multinomial2", is that correct? lmer, on the other hand, does not report residual variance with binomial model, apparently because it is fixed. For logit link function it seems to be fixed at 3.29, so the heritability formula will look like this: heritability = 4 * var.sire / (var.sire + 3.29) correct?... I feel this is way more than enough for one post... I guess I will ask about the models with sire-specific slopes later. cheers Misha Matz
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models The University of Glasgow, charity number SC004401
On Wed, 28 Dec 2011, mikhail matz wrote:
Hello Paul - the problem is, I am not sure that the built-in "animal" modeling in MCMCglmm fits my case. I am not measuring values of my traits in individual animals. Instead, I am measuring averages for full-sib families (coral larvae are quite small, so I am measuring them in batches). Even though I can build a pedigree for my families, denoting each full-sib family as an "animal", this would not be correct since, for example, "half-sibs" (in my case, crosses sharing one parent) would share 1/2 of the genes, not 1/4 as would be the case for individual half-sib measurements. Does it make sense?...
So you know the trait values for both parents? so you can do a midparent regression? And you don't or can't measure the within-sibship variance? Cheers, David Duffy.
| David Duffy (MBBS PhD) ,-_|\ | email: davidD at qimr.edu.au ph: INT+61+7+3362-0217 fax: -0101 / * | Epidemiology Unit, Queensland Institute of Medical Research \_,-._/ | 300 Herston Rd, Brisbane, Queensland 4029, Australia GPG 4D0B994A v
Hello David -
On Dec 28, 2011, at 4:15 PM, David Duffy wrote:
So you know the trait values for both parents? so you can do a midparent regression? And you don't or can't measure the within-sibship variance?
No, I cannot do midparent regression, unfortunately (this would be SO awesome). I am measuring larval traits, and I cannot measure those in wild-collected adult corals that serve as parents. I cannot measure within-family variance, either - only whole-family averages by measuring batches of larvae (the larvae are too small for most measurements). Please keep those ideas coming.... Misha
On Wed, 28 Dec 2011, mikhail matz wrote:
Hello Paul - the problem is, I am not sure that the built-in "animal" modeling in MCMCglmm fits my case. I am not measuring values of my traits in individual animals. Instead, I am measuring averages for full-sib families (coral larvae are quite small, so I am measuring them in batches). Even though I can build a pedigree for my families, denoting each full-sib family as an "animal", this would not be correct since, for example, "half-sibs" (in my case, crosses sharing one parent) would share 1/2 of the genes, not 1/4 as would be the case for individual half-sib measurements. Does it make sense?...
Cheers, David Duffy. -- | David Duffy (MBBS PhD) ,-_|\ | email: davidD at qimr.edu.au ph: INT+61+7+3362-0217 fax: -0101 / * | Epidemiology Unit, Queensland Institute of Medical Research \_,-._/ | 300 Herston Rd, Brisbane, Queensland 4029, Australia GPG 4D0B994A v