|
|
||||||||
,

,||,#
* Department of Animal Sciences, University of Wisconsin, Madison 53706
Department of Animal and Aquacultural Sciences, Norwegian University of Life Sciences, N-1432 Ås, Norway
Geno Breeding and A.I. Association, N-1432 Ås, Norway
Division of Genetic Epidemiology, Cancer Research UK Clinical Centre, St. Jamess University Hospital, Leeds, LS9 7TF, UK
|| Department of Dairy Science, and
# Department of Biostatistics and Medical Informatics, University of Wisconsin, Madison 53706
1 Corresponding author: nickwu{at}ansci.wisc.edu
| ABSTRACT |
|---|
|
|
|---|
Key Words: mastitis milk yield Monte Carlo method structural equation model
| INTRODUCTION |
|---|
|
|
|---|
Mixed models have been used to infer genetic and environmental correlations between milk yield (MY) and mastitis or SCC (e.g., Monardes and Hayes, 1985; Van Dorp et al., 1998; Heringstad et al., 2006) so as to understand their associations and to help optimize genetic selection programs for improving resistance to mastitis and milk yield simultaneously. These models, however, ignore direct relationships between the 2 phenotypes. High milk yield may increase liability to mastitis and, in turn, the disease may affect milk yield adversely. When mastitis occurs, a more or less severe short-term depression in milk yield may ensue. If mastitis remains unresolved, the effect may be long term, perhaps carrying on to the next lactation (Seegers et al., 2003). These types of direct effects between phenotypes have not been considered in classical mixed model approaches, but can be inferred using the simultaneous and recursive (SIR) models described in a quantitative genetics context by Gianola and Sorensen (2004). A SIR model is one of the many members included in the general concept of "structural equation models," where the main objective is to introduce causal pathways. Compared with a standard mixed model, a SIR model accounts for direct simultaneous and recursive effects between phenotypes. In a 2-trait system, for example, simultaneous effects refer to the presence of reciprocal direct effects between these 2 traits, whereas a recursive specification postulates that one trait affects the other directly, but the reverse does not occur. In dairy cattle, de los Campos et al. (2006b) investigated structural equation models for MY and SCC by using a 2-step likelihood-based approach, and because of the limitations of the software, they ignored covariances between levels of random effects attributable to genetic relationships. Using a similar strategy, de los Campos et al. (2006a) also investigated relationships between SCS and MY in dairy goats.
The objectives of the present research were 1) to extend the SIR models of Gianola and Sorensen (2004) to account for possibly heterogeneous relationships between phenotypes in subpopulations, and 2) to apply SIR models to describe relationships between MY and SCS in first-lactation Norwegian Red (NRF) cows. Competing SIR models were compared using the Bayesian information criterion (BIC). The consequences of accommodating simultaneous and recursive relationships between phenotypes on estimates of genetic parameters (i.e., heritability and genetic correlations) were investigated as well.
| MATERIALS AND METHODS |
|---|
|
|
|---|
![]() | [1] |
where yi is an m x 1 vector containing all phenotypic values for the ith individual; ß and u are vectors containing all "fixed" and "random" effects (this distinction does not arise in a Bayesian context), respectively; Xi and Zi are known incidence matrices of appropriate dimensions, and ei is a vector of residuals. What makes [1] different from a mixed model is the m x m matrix of structural coefficients
, which has the general form
![]() | [2] |
where
jj' is a structural coefficient giving the rate of change of trait j with respect to trait j'. All diagonal elements of
are 1, because each trait is assumed to have no direct effect on itself. If data on all n individuals are stacked, equation [1] expands to
![]() | [3] |
It is assumed that u | G0
N(0,A
G0) and e | R0
N(0,I
R0), where G0 is a genetic covariance matrix, R0 is a residual covariance matrix, A is an additive relationship matrix, and
indicates the Kronecker product.
The above model is extended to allow for heterogeneous
in a population, such that there is a different structural coefficient matrix for each of the subpopulations underlying the heterogeneity. Consider s subpopulations, each with nk measured individuals and with structural coefficient matrix
k in subpopulation k. Then [3] becomes
![]() | [4] |
where yk,i is a vector of phenotypes for the ith individual in the kth subpopulation; Xk,i and Zk,j are the corresponding incidence matrices. The need for heterogeneous coefficients could be due to many factors. In dairy cattle, for example, relationships between MY and mastitis or SCC may vary with level of production. Consequently, an increase in MY could have a different effect on SCC in high-production and low-production cows. Concerns have already been raised that selection for greater milk production may deteriorate resistance to mastitis (e.g., Sordillo and Streicher, 2002).
Fully Conditional Posterior Distributions
In the Bayesian modeling, a multivariate normal prior distribution is assumed with mean vector 1
0 and covariance matrix I
2 for each of the structural coefficient vectors, such that
k |
0,
2
N(1
0, I
2), where
0, = 0
2 = 10000 are known hyperparameters. Likewise, each of the "fixed effects" is assumed to be normal, a priori, with mean ß0 = 0 and variance b2 = 10,000, so that ß | ß0, b2 ~ N(0,I x 1002). The prior distributions for genetic and residual variance-covariance matrices are assumed to be G0 |
G, VG
Wishart1(
G, VG) and R0 |
R, VR
Wishart1(
R, VR), respectively, where Wishart1(
, V) represents an inverse Wishart distribution with degrees of freedom parameter
and scale matrix V. We set
G =
R = 6,
![]() |
![]() |
As in Gianola and Sorensen (2004), given ß and u, the vectors
;kyk,i are assumed to be mutually independent, ssuch that the joint density of all
kyk,i is the product of n normal densities, each pertaining to a single individual. By transforming variables from
kyk,i to yk,i (the determinant of the Jacobian of the transformation is |
k|), the density of the data, given
, ß, u, and R0, is
![]() | [5] |
The structural coefficient matrices
k shown in [2] can be readily set up, given the values in
k. Then the joint posterior distribution of unknown parameters can be factorized as
![]() | [6] |
where
includes all structural coefficient vectors
k (k = 1, 2, ..., s), H collectively contains all known hyperparameters, and Hi denotes a set of hyperparameters affecting a group of parameters (i = ß, G0, R0,
k). Note that all parameters in [6] are assumed independent, a priori, except for the dependence between u and G0.
The fully conditional posterior distributions of the parameters can be derived as in Gianola and Sorensen (2004) by retaining the parts of [6] varying with the parameters (or groups of parameters) of interest and treating the remaining parts as known. In the following section, we present these fully conditional distributions directly and describe a Metropolis-Hastings sampler for heterogeneous structural coefficients, which is different from that in Gianola and Sorensen (2004).
Location Parameters.
The conditional posterior distribution of the location parameters ß and u is
![]() | [7] |
where
![]() |
and
![]() |
Sampling location parameters from [7] can be done either directly, or block by block (Sorensen and Gianola, 2002), or element by element (Wang et al., 1994), or using the single-pass method of Garcia-Cortés and Sorensen (1996).
Variance-Covariance Components.
The conditional posterior distribution of the genetic covariance matrix G0 is the inverse Wishart distribution:
![]() | [8] |
where q is the number of levels of random effects, and all' is the element in row l and column l' of the inverse matrix of A. Similarly, the conditional posterior distribution of the residual covariance matrix R0 is the inverse Wishart distribution:
![]() | [9] |
where ek,i =
kyk,i Xk,iß Zk,iu. If G0 and R0 are assumed to be diagonal, both [8] and [9] reduce to products of independently scaled inverse
2 densities.
Structural Coefficients.
Although the "system" location and dispersion parameters pertain to the entire population, structural coefficients are specific to each of the subpopulations. Consider a subpopulation k. Retaining in [6] the parts that vary with structural coefficients in
k yields
![]() |
where
k includes all elements in
except
k. This implies that the
k vectors are conditionally independent. Thus, the fully conditional distribution of the structural coefficient vector can be shown to have the form
![]() | [10] |
where
![]() |
![]() |
and Yk,i is a working m x [m(m 1)] matrix constructed by noting that
kyk,i = yk,i Yk,i
k for each individual (Gianola and Sorensen, 2004).
Because [10] does not have a recognized form, a Metropolis-Hastings algorithm is proposed to sample elements of
k. Briefly, suppose values of
k for the Markov chain are moving from iteration t 1 to iteration t. Let
k(t1) be a vector containing current values (i.e., at iteration t 1), and let
k* be a vector containing proposed values that are generated using the following proposal density at iteration t:
![]() |
where
is a constant used to empirically tune the acceptance probability within a desired range, usually between 0.2 and 0.5 (Robert, 1995; Jannink and Wu, 2003). Then
k* is accepted with probability
![]() | [11] |
where ELSE refers to the values of all parameters other than those being updated. Note that [11] implies a symmetric move, J(
k* |
k(t1) = J(
k(t1) |
k*), so that the proposal ratio cancels out. If |
k| = 1 (e.g., as in a fully recursive specification), [10] reduces to a multivariate normal distribution, from which the vector
k can be sampled using a Gibbs sampler based on the process
![]() | [12] |
Inferring Genetic (Co)variance Components
In the presence of simultaneity and recursiveness, location and dispersion parameters from SIR models differ from those based on a mixed model, and they must be interpreted as "system parameters." Provided that the matrix
is invertible, one may estimate genetic variance and covariance components in the mixed model from posterior samples of "system" genetic variance and covariance components. Consider a sire model, for example, where the random effect vector u in model [4] is replaced by a sire effect vector s. Let G0 and R0 be the "system" between-sire genetic and residual variance-covariance matrices, respectively. Let G0*, R0*, and P0* be matrices containing the between-sire genetic, residual, and phenotypic variance-covariance components, respectively. Then
![]() | [13] |
![]() | [14] |
![]() | [15] |
Note that heterogeneous structural coefficient matrices lead to heterogeneous genetic and residual matrices in subpopulations, even though the "system" covariance matrices are assumed homogeneous. With relationships [13] through [15], the usual genetic parameters, such as heritability, and genetic and residual correlations can be calculated at each Markov chain Monte Carlo (MCMC) step.
In a bivariate system, Gianola and Sorensen (2004) showed that the apportionment of variance for the first trait into genetic and environmental components depends nontrivially on the structural parameter
12 but not on
21, whereas the opposite is true for the second trait. In multivariate systems involving many traits, the effects of structural coefficients on genetic parameters can be very complex. For example, consider M1 in Figure 1
, and denote G0 = {gij | i = 1, ..., m; j = 1, ..., m}. In this model, because
|
![]() |
then, for example,
![]() |
and
![]() |
where g* 11 and g* 12 are elements in G0* =
1G0
'1. Similar expressions can be obtained for R* 0 =
1R0
'1. Thus, in a multiple-trait system, the apportionment of (co)variance into genetic and environmental components for a given trait depends on direct effects that other traits have on it, as well as on their variance and covariance components. The more traits, the more cryptic their relationships are.
Bayesian Model Comparison
In a Bayesian framework, competing SIR models can be compared using Bayes factors (Jeffreys, 1935; Kass and Raftery, 1995). Consider 2 models, Mj and Mk, with parameter vectors
j and
k, respectively. The Bayes factor in support of Mj over Mk is the ratio between the posterior and the prior odds of the 2 models:
![]() | [16] |
It follows that the Bayes factor is equal to the posterior odd, p(Mj | y)/p(Mk | y), when the prior odds are equal to 1 [i.e., p(Mj) = p(Mk)].
The Bayes factor is difficult to compute. Schwarz (1978) derived the BIC as a large sample approximation to twice the logarithm of the Bayes factor. When several models are compared, it is helpful to compare each of them, say Mj, against a baseline model, say M0. The BIC is calculated as
![]() | [17] |
where
is the average of MCMC sampled log-likelihoods, c is the number of saved MCMC samples, d is the dimension of the corresponding parameter vector
, and n is the sample size. When comparing 2 competing SIR models, say Mj and Mk, then
![]() | [18] |
Thus, by calculating BIC values, one first compares SIR models against the model assuming no simultaneous and recursive effects as the baseline model. Then SIR models can be compared by taking the difference between their BIC values, with the model having the smaller BIC values being preferred. Following Raftery (1995), the BIC difference (
BIC) is interpreted as follows: 0
BIC < 2 (weak evidence), 2
BIC < 6 (positive evidence), 6
BIC < 10 (strong evidence), and
BIC
10 (very strong evidence).
Data Description and Model Specifications
The data, which represented 23,626 first-lactation daughters of 245 NRF sires that had their first progeny test in 1991 and 1992, included test-day records for MY and SCC, and veterinary records (presence or absence) on CM. Only test-day records from the first 120 d of lactation were included and cows with missing test-day records were excluded from the analysis. The 120 d of lactation were arbitrarily divided into 2 equal-length periods: from calving to d 60 (period I), and from d 61 to 120 (period II). For each period, cows were assigned one SCS and one MY test-day record, which were closest in time to the midpoint of that period. Presence or absence of CM was scored based on whether the cow had a veterinary-administered CM treatment in a 15-d period prior to the test-day. Test-day SCC was transformed into SCS, SCS = ln[(SCC/mL) x 103], to achieve an approximately normal distribution. Age at first calving (age) was also considered as a factor in the model, and consisted of 15 classes with age <20 mo as the first class, age >32 mo as the last class, and other classes each representing a single month between the first and last class. Somatic cell score and MY had been adjusted previously for herd effects via empirical best linear unbiased predictions of herd effects obtained using univariate mixed linear models (de los Campos et al., 2006b). The adjusted mean (standard deviation) was 4.040 (1.106) units for SCS1, 3.909 (1.150) units for SCS2, 21.60 (3.55) kg for MY1, and 21.12 (3.48) kg for MY2, respectively, where the number after each variable indicates the lactation period to which it belongs. The frequency of CM1 and CM2 (over test-days) was approximately 3.1 and 1.0%, respectively.
The data were analyzed using 4 SIR models (M1 through M4) and a multivariate mixed model assuming no SIR effects (M0), as baseline. These models are illustrated in Figure 1
. For model M0 (Figure 1A
), it was assumed that 1) correlations existed between sire effects and between residual effects, 2) CM1 affected SCS1, 3) CM2 affected SCS2, and 4) age affected both MY1 and MY2. These assumptions were also applied to all SIR models. In the 4 SIR models (Figure 1B
), simultaneous effects between SCS and MY within each lactation period were assumed, but between-period effects differed among models: M1 = no between-period direct effects; M2 = a direct effect from MY1 to SCS2; M3 = a direct effect from SCS1 to MY2, and M4 = direct effects from MY1 to SCS2 and from SCS1 to MY2.
The analysis was carried out using the SirBayes package (version 1.0), which is programmed in Fortran 95 for Linux. The package is freely available from the authors upon request (nickwu{at}ansci.wisc.edu; Gianola{at}ansci.wisc.edu). For each model, 10 chains were run, each starting from different initial values that were randomly chosen within the parameter space. In preliminary runs, convergence of iterative simulations was monitored based on multiple sequences with over-dispersed starting points (Gelman et al., 2004). Briefly, suppose K parallel sequences are simulated, each of length T (after discarding the burn-in simulations). For each estimand
, between- and within-sequence variances, B and W, were computed respectively, as described in Gelman et al. (2004). Then the marginal posterior variance of the estimand was estimated by a weighted average of W and B, namely,
![]() | [19] |
Finally, convergence of the iterative simulation can be monitored by estimating the factor by which the scale of the current distribution for the estimand might be reduced if the simulations were continued in the limit T
. This potential scale reduction is estimated by
![]() | [20] |
which declines to 1 as T
. This method of monitoring convergence has a key advantage of not requiring the user to examine time series graphs of simulated sequences. Based on the convergence diagnostics (i.e.,
1 < 0.02), it was decided that a single chain should consist of 100,000 iterations. Posterior samples were thinned every 10 iterations after burn-in (i.e., discarding samples from the first 1,000 iterations) and saved for posterior inference. Genetic parameters were calculated for each thinned sample and were saved simultaneously with posterior samples of location and dispersion parameters. Effects of SIR were estimated for all cows jointly using model [3], and for the high (i.e., cows with milk yield greater than the average) and the low (i.e., cows with milk yield less than the average) groups of cows, distinctively, using model [4]. Naturally, using the data to form groups leads to sampling bias, so inferences about heterogeneity of relationships in animals with high and low yield potential may have some bias.
| RESULTS AND DISCUSSION |
|---|
|
|
|---|
The BIC values for each of the 4 SIR models are given in Table 1
. The BIC values ranged from 169.87 (M1) to 150.86 (M4). The large negative BIC values indicated that the SIR models were favored very strongly over the mixed model without SIR relationships (M0). Differences between BIC values indicated strong (
BIC
6.0) and very strong (
BIC
10.0) evidence favoring the model without any between-period effect (M1) over those with between-period effect(s). This was probably because between-period effects were small. Because BIC penalizes complex models by a quantity that is equal to the extra number of parameters times the logarithm of the sample size, it tends to favor simpler, more parsimonious models. Also for this reason, BIC differences strongly favored models M2 and M3 over model M4 (Table 1
).
|
|
|
|
Estimated direct effects from SCS to MY were in agreement with Koldeweij et al. (1999), who reported a decrease of 2.04 kg of milk per unit of increase in log10 (SCC x 103) based on a Swedish data set. Our estimates also agreed with de los Campos et al. (2006b), who estimated simultaneous effects between MY and SCS in the same population using the LISREL software package (Jöreskog and Sörbom, 2005). The method used in the current study differs from that of LISREL in a few respects. First, it is within the Bayesian framework, whereas LISREL uses maximum likelihood to estimate structural coefficients and other parameters. Second, LISREL assumes multivariate normality for all observables, including predictor variables, but this assumption does not apply in the current analysis; the method used in this study assumes joint multivariate normality of the response variables and "random" effects only. Third, LISREL does not accommodate pedigree information, whereas such information can be included in the present method. Nevertheless, both studies obtained similar estimates of direct effects from SCS to MY, yet with some difference. The present study suggested a small, positive reciprocal effect (e.g., from MY to SCS) in period I, which did not appear to be significant in de los Campos et al. (2006b). This difference could be because they estimated direct effects as the average for the whole first lactation, whereas our estimates represented the first 2 periods of lactation, d 0 to 60 and d 61 to 120, respectively. Possibly for the same reason, estimates in period II were aligned with estimates from de los Campos et al. (2006b), but the absolute values of estimates in period I were larger than estimates from de los Campos et al. (2006b). Finally, the present results suggested possible heterogeneity of simultaneous and recursive effects in high- and low-yield cows; this information was not available in de los Campos et al. (2006b).
Heritability and Genetic Correlations
Estimates of heritability based on SIR models were close to those obtained using the mixed model, yet with some minor differences (Figure 4
). For example, posterior mean (standard deviation) of heritability was from 0.12 (0.01) to 0.13 (0.01) for SCS1 and SCS2, regardless of the models. Similarly, posterior mean (standard deviation) of heritability of MY1 was approximately 0.18 (0.02) based on either M0 or SIR models (M1 through M4). Estimated heritability of MY2 was 0.22 (0.02) based on M0, and it ranged from 0.22 (0.02) to 0.23 (0.02) based on SIR models. In general, heritability estimates of MY were larger at lactation period II than at period I.
|
Posterior estimates of the genetic correlation between MY1 and MY2 were similar in all models (Figure 5A
), with the posterior mean (standard deviation) being approximately 0.80 (0.02 to 0.03). Posterior distributions of the genetic correlation between SCS1 and SCS2 showed some differences between models (Figure 5B
). The posterior distributions obtained using SIR models had a lower mean and a larger standard deviation than those from the mixed model (Figure 5B
). The genetic correlation (standard deviation) between SCS1 and SCS2 was 0.68 (0.03) for model M0 and ranged from 0.58 (0.06) to 0.63 (0.04) for the SIR models. Estimated genetic correlations between SCS1 and MY2 were similar across models, with the posterior mean (standard deviation) being 0.21 (0.06) for M0, and 0.22 to 0.24 (0.10 to 0.12) for SIR models.
|
| CONCLUSIONS |
|---|
|
|
|---|
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
Received for publication November 15, 2006. Accepted for publication March 7, 2007.
| REFERENCES |
|---|
|
|
|---|
This article has been cited by other articles:
![]() |
B. Heringstad, X.-L. Wu, and D. Gianola Inferring relationships between health and fertility in Norwegian Red cows using recursive models J Dairy Sci, April 1, 2009; 92(4): 1778 - 1784. [Abstract] [Full Text] [PDF] |
||||
![]() |
E. L. de Maturana, X.-L. Wu, D. Gianola, K. A. Weigel, and G. J. M. Rosa Exploring Biological Relationships Between Calving Traits in Primiparous Cattle with a Bayesian Recursive Model Genetics, January 1, 2009; 181(1): 277 - 287. [Abstract] [Full Text] [PDF] |
||||
![]() |
A. B. Samore, A. F. Groen, P. J. Boettcher, J. Jamrozik, F. Canavesi, and A. Bagnato Genetic Correlation Patterns Between Somatic Cell Score and Protein Yield in the Italian Holstein-Friesian Population J Dairy Sci, October 1, 2008; 91(10): 4013 - 4021. [Abstract] [Full Text] [PDF] |
||||
![]() |
S. Konig, X. L. Wu, D. Gianola, B. Heringstad, and H. Simianer Exploration of Relationships Between Claw Disorders and Milk Yield in Holstein Cows via Recursive Linear and Threshold Models J Dairy Sci, January 1, 2008; 91(1): 395 - 406. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |