JDS
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


J. Dairy Sci. 2007. 90:3508-3521. doi:10.3168/jds.2006-762
© 2007 American Dairy Science Association ®

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Interpretive Summary
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Wu, X.-L.
Right arrow Articles by Gianola, D.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Wu, X.-L.
Right arrow Articles by Gianola, D.

Inferring Relationships Between Somatic Cell Score and Milk Yield Using Simultaneous and Recursive Models

X.-L. Wu*,1, B. Heringstad{dagger},{ddagger}, Y.-M. Chang§, G. de los Campos* and D. Gianola*,{dagger},||,#

* Department of Animal Sciences, University of Wisconsin, Madison 53706
{dagger} Department of Animal and Aquacultural Sciences, Norwegian University of Life Sciences, N-1432 Ås, Norway
{ddagger} Geno Breeding and A.I. Association, N-1432 Ås, Norway
§ Division of Genetic Epidemiology, Cancer Research UK Clinical Centre, St. James’s 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
 TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
A Bayesian analysis via Markov chain Monte Carlo methods extending the simultaneous and recursive model of Gianola and Sorensen (2004) was proposed to account for possible population heterogeneity. The method was used to infer relationships between milk yield and somatic cell scores of Norwegian Red cows. Data consisted of test-day records of milk yield and somatic cell score of first-lactation cows during the first 120 d of lactation. Results suggested large negative direct effects from somatic cell score to milk yield and small reciprocal effects from milk yield to somatic cell score. The direct effects were larger in the first 60 d of lactation than in the subsequent period. Bayesian model selection strongly favored the simultaneous and recursive models for milk yield and somatic cell score over the corresponding mixed model without considering simultaneity or recursiveness. Estimated effects between milk yield and somatic cell score seemed to be yield-dependent, larger in higher producing cows than in lower producing cows. Heritability estimates from the simultaneous and recursive models were similar to those from the mixed model, but some genetic correlations differed considerably among models.

Key Words: mastitis • milk yield • Monte Carlo method • structural equation model


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
Mastitis, a bacteria-related inflammation of the mammary gland, is a serious and costly health problem in dairy cattle worldwide (Hillerton and Berry, 2005). Recording of clinical mastitis (CM) cases is often difficult in practice; further, most mastitis occurs as a low-grade infection, a subclinical state. In dairy cattle, SCC has been used as an indirect measurement of mastitis, because it is easy to record and it has a reasonably large genetic correlation with mastitis (Shook and Schutz, 1994). Somatic cells are primarily polymorphonuclear leukocytes, which migrate to the mammary gland in large numbers during infection. Thus, an elevated SCC in milk is indicative of clinical or subclinical mastitis.

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
 TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
Simultaneous and Recursive (SIR) Models
Consider m traits that are measured in n individuals. The SIR mixed model for data on individual i (Gianola and Sorensen, 2004) can be represented as


Formula 1[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 {Lambda}, which has the general form


Formula 2[2]

where {lambda}jj' is a structural coefficient giving the rate of change of trait j with respect to trait j'. All diagonal elements of {Lambda} 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


Formula 3[3]

It is assumed that u | G0 ~ N(0,A {otimes} G0) and e | R0 ~ N(0,I {otimes} R0), where G0 is a genetic covariance matrix, R0 is a residual covariance matrix, A is an additive relationship matrix, and {otimes} indicates the Kronecker product.

The above model is extended to allow for heterogeneous {lambda} 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 {Lambda}k in subpopulation k. Then [3] becomes


Formula 4[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{lambda}0 and covariance matrix I{tau}2 for each of the structural coefficient vectors, such that {lambda}k | {lambda}0, {tau}2 ~ N(1{lambda}0, I{tau}2), where {lambda}0, = 0 {tau}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 | {nu}G, VG ~ Wishart–1({nu}G, VG) and R0 | {nu}R, VR ~ Wishart–1({nu}R, VR), respectively, where Wishart–1({nu}, V) represents an inverse Wishart distribution with degrees of freedom parameter {nu} and scale matrix V. We set {nu}G = {nu}R = 6,


Formula 4


Formula 4

As in Gianola and Sorensen (2004), given ß and u, the vectors {Lambda};kyk,i are assumed to be mutually independent, ssuch that the joint density of all {Lambda}kyk,i is the product of n normal densities, each pertaining to a single individual. By transforming variables from {Lambda}kyk,i to yk,i (the determinant of the Jacobian of the transformation is |{Lambda}k|), the density of the data, given {lambda}, ß, u, and R0, is


Formula 5[5]

The structural coefficient matrices {Lambda}k shown in [2] can be readily set up, given the values in {lambda}k. Then the joint posterior distribution of unknown parameters can be factorized as


Formula 6[6]

where {lambda} includes all structural coefficient vectors {lambda}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, {lambda}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


Formula 7[7]

where


Formula 7

and


Formula 7

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:


Formula 8[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:


Formula 9[9]

where ek,i = {Lambda}kyk,iXk,iß – Zk,iu. If G0 and R0 are assumed to be diagonal, both [8] and [9] reduce to products of independently scaled inverse {chi}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 {lambda}k yields


Formula 9

where {lambda}k includes all elements in {lambda} except {lambda}k. This implies that the {lambda}k vectors are conditionally independent. Thus, the fully conditional distribution of the structural coefficient vector can be shown to have the form


Formula 10[10]

where


Formula 10


Formula 10

and Yk,i is a working m x [m(m – 1)] matrix constructed by noting that {Lambda}kyk,i = yk,iYk,i{lambda}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 {lambda}k. Briefly, suppose values of {lambda}k for the Markov chain are moving from iteration t – 1 to iteration t. Let {lambda}k(t–1) be a vector containing current values (i.e., at iteration t – 1), and let {lambda}k* be a vector containing proposed values that are generated using the following proposal density at iteration t:


Formula 10

where {delta} 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 {lambda}k* is accepted with probability


Formula 11[11]

where ELSE refers to the values of all parameters other than those being updated. Note that [11] implies a symmetric move, J({lambda}k* | {lambda}k(t–1) = J({lambda}k(t–1) | {lambda}k*), so that the proposal ratio cancels out. If |{Lambda}k| = 1 (e.g., as in a fully recursive specification), [10] reduces to a multivariate normal distribution, from which the vector {lambda}k can be sampled using a Gibbs sampler based on the process


Formula 12[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 {Lambda} 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


Formula 13[13]


Formula 14[14]


Formula 15[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 {lambda}12 but not on {lambda}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 1Go, and denote G0 = {gij | i = 1, ..., m; j = 1, ..., m}. In this model, because


Figure 1
View larger version (12K):
[in this window]
[in a new window]

 
Figure 1. Five models used in the analysis: (A) M0 = a mixed model assuming no simultaneous and recursive effects, and (B) M1–M4 = 4 possible relationships between milk yield (MY) and SCS. AGE = age at first calving; CM1 and CM2 = clinical mastitis at the first and second lactation periods; SCS1 and SCS2 = SCS in the first and second lactation periods; MY1 and MY2 = milk yield in the first and second lactation periods; S1 to S4 = sire effects; E1 to E4 = residual effects. A single headed arrow (->) indicates a causal relationship, 2 reciprocal arrows ({rightleftarrows}) indicate simultaneous effects, and a double headed connector ({leftrightarrow}) represents a correlation.

 

Formula 15

then, for example,


Formula 15

and


Formula 15

where g* 11 and g* 12 are elements in G0* = {Lambda}1G0{Lambda}'–1. Similar expressions can be obtained for R* 0 = {Lambda}1R0{Lambda}'–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 {theta}j and {theta}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:


Formula 16[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


Formula 17[17]

whereFormula 17 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 {theta}, and n is the sample size. When comparing 2 competing SIR models, say Mj and Mk, then


Formula 18[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 ({Delta}BIC) is interpreted as follows: 0 ≤ {Delta}BIC < 2 (weak evidence), 2 ≤ {Delta}BIC < 6 (positive evidence), 6 ≤ {Delta}BIC < 10 (strong evidence), and {Delta}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 10–3], 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 1Go. For model M0 (Figure 1AGo), 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 1BGo), 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 {phi}, 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,


Formula 19[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 -> {infty}. This potential scale reduction is estimated by


Formula 20[20]

which declines to 1 as T -> {infty}. 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., Formula 20 – 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
 TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
Bayesian Assessment of SIR Models
Comparing models with different levels of simultaneity and recursiveness is challenging because the number of models potentially entertained can be very large. In SIR models, as illustrated in Figure 1BGo, model selection can be expressed in terms of which of the arrow(s) are activated (added) or deactivated (removed). The number of arrows, however, increases dramatically with phenotypic variables (i.e., traits). Considering 4 traits, for example, permutating the 4 dependent variables leads to 24(4–1) = 4,096 competing models to explain the relationships between these 4 phenotypic variables. Thus, comparing a complete set of competing models is not feasible when there are many traits in the analysis. Nevertheless, one can make use of theory or prior knowledge to substantially reduce the number of models under consideration. In our application, we assumed a priori that simultaneous relationships between SCS and MY would exist within periods without being sure about influences between lactation periods. This allowed us to reduce the number of competing models to only 4 in the current study.

The BIC values for each of the 4 SIR models are given in Table 1Go. 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 ({Delta}BIC ≥ 6.0) and very strong ({Delta}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 1Go).


View this table:
[in this window]
[in a new window]

 
Table 1. Model comparison based on Bayesian information criterion (BIC) differences1
 
Direct Effects Between SCS and MY
The actual acceptance rate for structural coefficients ranged from 0.1670 to 0.1697, which were slightly lower than the recommended range between 0.2 and 0.5 (Robert, 1995; Jannink and Wu, 2003). The lower acceptance rate may be because all the structural coefficients were accepted or rejected as a block. Estimated SIR effects between MY and SCS are given in Table 2Go. Direct effects from SCS to MY were strong and negative (approximately –1.8 kg per unit of SCS in period I, and from –0.9 to –1.0 kg per unit of SCS in period II), whereas the reciprocal effects from MY to SCS were small and positive in period I and not different from zero in period II. Results in period I indicate that an increase in MY augments the risk of mastitis, and mastitis increases SCC. Because CM was treated as a covariate in the model equations for SCS, the elevation of SCC as MY increased may be indicating a higher risk of subclinical infection as yield increases. Posterior distributions of direct effects from SCS to MY were similar for the 4 SIR models in both periods, and posterior means were larger in absolute value in period I than in period II (Figure 2AGo). Posterior estimates of direct effects from MY1 to SCS1 in period I were all centered at approximately 0.04 units of SCS/kg of milk yield. Their counterparts in period II showed minor differences between models (Figure 2BGo), which were probably due to chance, because 0 entered with appreciable density in the posterior distributions. On average, the absolute values of direct effect from SCS to MY decreased by 47% from the first to the second lactation period, whereas the effect from MY to SCS decreased from being mildly positive in period I to near zero in period II. This is reasonable because more than 60% of mastitis cases occur during the first 60 d after calving (Heringstad et al., 1999).


View this table:
[in this window]
[in a new window]

 
Table 2. Posterior mean (SD) of structural coefficients obtained using simultaneous and recursive (SIR) mixed models1
 

Figure 2
View larger version (20K):
[in this window]
[in a new window]

 
Figure 2. Posterior distributions of simultaneous and recursive effects at 2 lactation periods for (A) a direct effect from SCS to milk yield (MY) and (B) a direct effect from MY to SCS, for the following models: Model 1 (M1) = assuming simultaneous effects between SCS and MY within each lactation period; M2 = M1 + a recursive effect from MY1 to SCS2; M3 = M1 + a recursive effect from SCS1 to MY2; M4 = M1 + a recursive effect from MY1 to SCS2 + a recursive effect from SCS1 to MY2. Effect from MY to SCS in units of SCS per kilogram of milk yield; effect from SCS to MY in kilograms of milk per unit of SCS.

 
Estimated simultaneous effects between MY and SCS were different in the high (above average) and low (below average) milk yield groups. As shown in Figure 3A and 3BGo, the posterior distributions of direct effects between MY and SCS for the high group were sharper (i.e., less variation), and posterior means had larger absolute values than those in the low group. The direct effects of SCS on MY in the high group were more strongly negative than those in the low group. This may indicate that the effect of subclinical mastitis (CM was accounted for in the model) on MY was stronger in high-producing cows, although it may reflect scale effects as well. A negative correlation between milk production and resistance to mastitis has been reported (Detilleux et al., 1995). It is likely that resistance to mastitis will deteriorate in populations that are selected for increased production (Sordillo and Streicher, 2002) unless breeding objectives place some weight on resistance to the disease. However, results from this study would suggest that the strength of this association may be yield-dependent.


Figure 3
View larger version (21K):
[in this window]
[in a new window]

 
Figure 3. Posterior distributions of (A) a direct effect from SCS to milk yield (MY) and (B) a direct effect from MY to SCS for the high (above average) and low (below average) milk yield groups, obtained from the simultaneous model, M1. Effect from MY to SCS in units of SCS per kilogram of milk yield; effect from SCS to MY in kilograms of milk per unit of SCS.

 
Between-period effects are given in Table 2Go. The posterior means for the effect from MY1 to SCS2 were positive, and the effect from SCS1 to MY2 was negative. However, these effects cannot be considered statistically significant, because their 95% credible intervals included zero.

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 10–3) 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 4Go). 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.


Figure 4
View larger version (26K):
[in this window]
[in a new window]

 
Figure 4. Posterior distributions of heritability for (A) SCS (SCS1) and milk yield (MY1) in the first lactation period, and (B) SCS (SCS2) and milk yield (MY2) in the second lactation period for the following models: M0 = mixed model assuming no simultaneous and recursive effects; M1 = assuming simultaneous effects between SCS and MY within each lactation period; M2 = M1 + a recursive effect from MY1 to SCS2; M3 = M1 + a recursive effect from SCS1 to MY2; M4 = M1 + a recursive effect from MY1 to SCS2 + a recursive effect from SCS1 to MY2.

 
The estimates of heritability of SCS were in agreement with previous reports in the same population. In NRF, Ødegård et al. (2004) found that heritability of SCS was 0.11, and de los Campos et al. (2006b) reported estimates of heritability of SCS ranging from 0.09 to 0.11. However, heritability estimates from this study for test-day milk yield in the first 60-d period of lactation (MY1) were higher than those for the entire first-lactation period obtained by de los Campos et al. (2006b) in the same population, possibly suggesting larger genetic variation (relative to environmental variance) at early lactation stages.

Posterior estimates of the genetic correlation between MY1 and MY2 were similar in all models (Figure 5AGo), 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 5BGo). The posterior distributions obtained using SIR models had a lower mean and a larger standard deviation than those from the mixed model (Figure 5BGo). 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.


Figure 5
View larger version (16K):
[in this window]
[in a new window]

 
Figure 5. Posterior distributions of genetic correlations (A) between milk yield in the first (MY1) and second (MY2) periods, (B) between SCS in the first (SCS1) and second (SCS2) periods, and (C) between MY1 and SCS2, for the following models: M0 = a standard mixed model; M1 = assuming simultaneous effects between SCS and MY within each lactation period; M2 = M1 + a recursive effect from MY1 to SCS2; M3 = M1 + a recursive effect from SCS1 to MY2; M4 = M1 + a recursive effect from MY1 to SCS2 + a recursive effect from SCS1 to MY2.

 
A larger difference was observed in the posterior distributions of the genetic correlation between MY1 and SCS2 (Figure 5CGo); the posterior means (standard deviations) were 0.19 (0.06) for model M0 and 0.52 to 0.53 (0.06 to 0.08) for the SIR models. Moderate to large differences in genetic correlations were also observed between SCS1 and MY1 and between SCS2 and MY2. The genetic correlation based on M0 was 0.14 (0.06) between SCS1 and MY1, and 0.19 (0.06) between SCS2 and MY2, but their counterparts based on SIR models were approximately twice as large. All genetic correlations between MY and SCS were positive. This corroborates the idea that increasing milk yield by selection would deteriorate resistance to subclinical mastitis. However, there is some evidence that the genetic correlation between milk yield and SCC is close to zero in the second lactation (Koivula et al., 2005).


    CONCLUSIONS
 TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
The simultaneous and recursive models of Gianola and Sorensen (2004) were extended to account for population heterogeneity. The extended method assumed different simultaneous and recursive effects for different subpopulations, leading to stratum-specific estimates of genetic and residual variance-covariance components. Proposed within the Bayesian framework via an MCMC implementation, this method has provided greater modeling flexibility than maximum likelihood approaches. This method was applied to estimate direct effects between SCS and MY in first-lactation NRF cows. The results suggested large negative effects from SCS to MY, and small or nonsignificant reciprocal effects from MY to SCS. These direct effects decreased from lactation period I to period II, and differed between high and low milk yield groups of cows. Heritability estimates from SIR models were similar to those from the mixed models, but some genetic correlations differed considerably among models.


    ACKNOWLEDGEMENTS
 TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
This research was financially supported by the Babcock Institute for International Dairy Research and Development, University of Wisconsin-Madison, and by grants NRICGP/USDA 2003-35205-12833, NSF DEB-0089742, and NSF DMS-044371. Access to the data was given by the Norwegian Dairy Herd Recording System and the Norwegian Cattle Health Service in agreement number 004.2005. We thank 2 anonymous reviewers and the editors for their helpful comments.

Received for publication November 15, 2006. Accepted for publication March 7, 2007.


    REFERENCES
 TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 


de los Campos, G., D. Gianola, P. Boettcher, and P. Moroni. 2006a. A structural equation model for describing relationships between somatic cell score and milk yield in dairy goats. J. Anim. Sci. 84:2934–2941.[Abstract/Free Full Text]

de los Campos, G., D. Gianola, and B. Heringstad. 2006b. A structural equation model for describing relationships between somatic cell count and milk yield in first-lactation dairy cows. J. Dairy Sci. 89:4445–4455.[Abstract/Free Full Text]

Detilleux, J. C., M. E. Kehrli, J. R. Stabel, A. E. Freeman, and D. H. Kelley. 1995. Study of immunological dysfunction in periparturient Holstein cattle selected for high and average milk production. Vet. Immunol. Immunopathol. 44:251–267.[CrossRef][Medline]

Garcia-Cortés, L. A., and D. Sorensen. 1996. On a multivariate implementation of the Gibbs sampler. Genet. Sel. Evol. 28:121–126.[CrossRef]

Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin. 2004. Bayesian Data Analysis. 2nd ed. Chapman & Hall/CRC, London, UK.

Gianola, D., and D. Sorensen. 2004. Quantitative genetic models for describing simultaneous and recursive relationships between phenotypes. Genetics 167:1407–1424.[Abstract/Free Full Text]

Heringstad, B., D. Gianola, Y. M. Chang, J. Ødegård, and G. Klemetsdal. 2006. Genetic associations between clinical mastitis and somatic cell score in early first-lactation cows. J. Dairy Sci. 89:2236–2244.[Abstract/Free Full Text]

Heringstad, B., G. Klemetsdal, and J. Ruane. 1999. Clinical mastitis in Norwegian cattle: Frequency, variance components, and genetic correlation with protein yield. J. Dairy Sci. 82:1325–1330.[Abstract]

Hillerton, J. E., and E. A. Berry. 2005. Treating mastitis in the cow—A tradition or an archaism. J. Appl. Microbiol. 98:1250–1255.[CrossRef][Medline]

Jannink, J.-L., and X.-L. Wu. 2003. Estimating allelic number and identity in state of QTLs in interconnected families. Genet. Res. 81:133–144.[CrossRef][Medline]

Jeffreys, H. 1935. Some tests of significance, treated by the theory of probability. Proc. Camb. Philol. Soc. 31:203–222.

Jöreskog, K., and D. Sörbom. 2005. LISREL for Windows [computer software]. Scientific Software International, Lincolnwood, IL.

Kass, R. E., and A. E. Raftery. 1995. Bayes factors. J. Am. Stat. Assoc. 90:773–795.[CrossRef]

Koivula, M., E. A. Mantysaari, E. Negussie, and T. Serenius. 2005. Genetic and phenotypic relationships among milk yield and somatic cell count before and after clinical mastitis. J. Dairy Sci. 88:827–833.[Abstract/Free Full Text]

Koldeweij, E., U. Emanuelson, and L. Janson. 1999. Relation of milk production loss to milk somatic cell count. Acta Vet. Scand. 40:47–56.[Medline]

Monardes, H. G., and J. F. Hayes. 1985. Genetic and phenotypic relationships between lactation cell counts and milk yield and composition of Holstein cows. J. Dairy Sci. 68:1250–1256.[Abstract/Free Full Text]

Ødegård, J., B. Heringstad, and G. Klemetsdal. 2004. Short communication: Bivariate genetic analysis of clinical mastitis and somatic cell count in Norwegian dairy cattle. J. Dairy Sci. 87:3515–3517.[Abstract/Free Full Text]

Raftery, A. E. 1995. Bayesian model selection in social research. Sociol. Methodol. 25:111–163.[Medline]

Robert, G. O. 1995. Markov chain concepts related to sampling algorithms. Pages 45–58 in Markov Chain Monte Carlo in Practice. W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, ed. Chapman & Hall, London, UK.

Schwarz, G. 1978. Estimating the dimension of a model. Ann. Stat. 6:461–464.[CrossRef]

Seegers, H., C. Fourichon, and F. Beaudeau. 2003. Production effects related to mastitis and mastitis economics in dairy cattle herds. Vet. Res. 34:475–491.[CrossRef][Medline]

Shook, G. E., and M. M. Schutz. 1994. Selection on somatic cell score to improve resistance to mastitis in the United States. J. Dairy Sci. 77:648–658.[Abstract]

Sordillo, L. M., and K. L. Streicher. 2002. Mammary gland immunity and mastitis susceptibility. J. Mammary Gland Biol. Neoplasia 7:135–146.[CrossRef][Medline]

Sorensen, D., and D. Gianola. 2002. Likelihood, Bayesian, and MCMC Methods in Quantitative Genetics. Springer-Verlag, New York, NY.

Van Dorp, T. E., J. C. Dekkers, S. W. Martin, and J. P. Noordhuizen. 1998. Genetic parameters of health disorders, and relationships with 305-day milk yield and conformation traits of registered Holstein cows. J. Dairy Sci. 81:2264–2270.[Abstract]

Wang, C. S., J. J. Rutledge, and D. Gianola. 1994. Bayesian analysis of mixed linear models via Gibbs sampling with an application to litter size in Iberian pigs. Genet. Sel. Evol. 26:91–115.[CrossRef]


This article has been cited by other articles:


Home page
J DAIRY SCIHome page
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]


This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Interpretive Summary
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Wu, X.-L.
Right arrow Articles by Gianola, D.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Wu, X.-L.
Right arrow Articles by Gianola, D.


HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS