[R-sig-ME] Binary response animal models in MCMCglmm

[Jacob Berson]

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]]