Exploration of lagged relationships between mastitis and milk yield in dairycows using a Bayesian structural equation Gaussian-threshold model

A Gaussian-threshold model is described under the general framework of structural equation models for inferring simultaneous and recursive relationships between binary and Gaussian characters, and estimating genetic parameters. Relationships between clinical mastitis (CM) and test-day milk yield (MY) in first-lactation Norwegian Red cows were examined using a recursive Gaussian-threshold model. For comparison, the data were also analyzed using a standard Gaussian-threshold, a multivariate linear model, and a recursive multivariate linear model. The first 180 days of lactation were arbitrarily divided into three periods of equal length, in order to investigate how these relationships evolve in the course of lactation. The recursive model showed negative within-period effects from (liability to) CM to test-day MY in all three lactation periods, and positive between-period effects from test-day MY to (liability to) CM in the following period. Estimates of recursive effects and of genetic parameters were time-dependent. The results suggested unfavorable effects of production on liability to mastitis, and dynamic relationships between mastitis and test-dayMYin the course of lactation. Fitting recursive effects had little influence on the estimation of genetic parameters. However, some differences were found in the estimates of heritability, genetic, and residual correlations, using different types of models (Gaussian-threshold vs. multivariate linear).


INTRODUCTION
Multivariate linear models have long been used for multiple-trait genetic evaluation and analysis e.g. [2,18,24]. However, these standard models do not Let y c i ¼ y i;1 ::: y i;t 1 À Á be a vector containing observations for the t 1 continuous characters of the ith individual. Let g i ¼ g i;t 1 þ1 ::: g i;t 1 þt 2 À Á be a vector containing the t 2 binary variables (observable scale) of the ith individual, and y b i ¼ y i;t 1 þ1 ::: y i;t 1 þt 2 À Á be a vector containing the corresponding liability variables (underlying scale), which are assumed to be continuous and normally distributed. The theory of threshold models states that for a binary character, the phenotype of an individual is 1 (e.g., sick) if the underlying liability exceeds a threshold j b and 0 (e.g., healthy) otherwise, so where b ¼ t 1 þ 1; :::; t 1 þ t 2 . The threshold is fixed arbitrarily to center the distribution, so it is not an unknown parameter in a binary threshold model. Let g and y b be vectors containing all binary observations and underlying liabilities, respectively, of all individuals, and let j ¼ j t 1 þ1 ::: j t 1 þt 2 ð Þbe a vector that contains the thresholds for all binary traits. Then, the conditional probability of observing a realization of g, given y b and j, is given by where I(A) is an indicator function, which takes the value 1 if condition A is true and 0 otherwise. Next, consider the joint distribution of the continuous phenotypes and of the liabilities of the binary characters. The unknown liabilities are treated as nuisance parameters, after data augmentation, in the second step of the multilevel modeling. Note that y i ¼ y c i 0 y b i 0 À Á 0 . Assume, further, that variables in y i are affected mutually, so that a phenotype or liability is a linear function of other phenotypes or liabilities, as well as of ''fixed'' and random effects that are relevant. Then, the model is Here, b is a vector of fixed effects; u is a vector of genetic effects; c is a vector of environmental effects (e.g., herds); e i,j is a random residual; x 0 i;j , z 0 i;j , w 0 i;j are incidence row vectors pertaining to the jth trait of the ith individual, and k j,j 0 is Bayesian structural equation Gaussian-threshold model an unknown structural coefficient (i.e., regression coefficient of phenotype j on phenotype or liability j 0 ). If all k's are equal to 0, then (3) is a standard linear model. In matrix form, (3) can be expressed as K y i;1 Á Á Á y i;t 1 y i;t 1 þ1 Á Á Á y i;t 1 þt 2 where X i ¼ x i;1 ::: x i;t 1 x i;t 1 þ1 ::: x i;t 1 þt 2 ð Þ 0 ; Z i ¼ z i;1 ::: z i;t 1 z i;t 1 þ1 ::: z i;t 1 þt 2 ð Þ 0 ; W i ¼ w i;1 ::: w i;t 1 w i;t 1 þ1 ::: w i;t 1 þt 2 ð Þ 0 ; e i ¼ e i;1 ::: e i;t 1 e i;t 1 þ1 ::: e i;t 1 þt 2 ð Þ 0 : The K matrix is a structural coefficient matrix, in which a diagonal element is 1 and an off-diagonal element is Àk jj 0 (j 6 ¼ j 0 ). The conditional distribution of Ky i is assumed multivariate normal, such that or, by changing variables where R 0 is a residual variance-covariance matrix, and k is a vertical concatenation of all off-diagonal elements of K. Conditionally on b, u, and c, the Ky i 's are mutually independent. The same is true of the y i 's, given b, u, c, and k. Thus, For this hierarchical model, the joint distribution of all observed data (including binary scores) and liabilities is Note that, given the liabilities and the thresholds, the vector of discrete outcomes g is independent of y c , the Gaussian phenotypes.

Prior distributions
Following Gianola and Sorensen [10], we assigned multivariate normal prior distributions to structural coefficients and ''fixed'' effects. By assuming an infinitesimal model, the prior distribution of genetic effects is multivariate normal with an unknown genetic covariance matrix G 0 , ujA; G 0 $ N ð0; A G 0 Þ, where A is the additive relationship matrix and represents the Kronecker product. Similarly, the prior distribution of the environmental effects vector is c $ N 0; I D 0 ð Þ , where D 0 is a variance-covariance matrix among environmental effects. The prior distributions of the genetic, environmental, and residual covariance matrices are assumed to be inverted Wishart, Wishart À1 t k ; V k ð Þ, with scaling matrix V k and degrees of freedom parameter t k , where k ¼ G 0 ; D 0 ; R 0 .

Joint posterior distributions
f g be the parameters of the model. The posterior distribution is augmented with the unobserved liabilities such that the joint posterior distribution of all unobservables is where H represents the collection of all known hyper-parameters, and, for example, pðbjH b Þ is the density of the prior distribution of b and H b is a set of known hyper-parameters (i.e., mean b 0 and variance r 2 b 0 ) that the distribution of b depends on.
Bayesian structural equation Gaussian-threshold model

Fully conditional posterior distributions
The fully conditional posterior distributions can be ascertained from (9) by retaining the parts varying with the parameter or group of parameters of interest and treating the remaining parts as known [21].

Liabilities
To obtain the fully conditional posterior distribution of the liability variable (y i,b ) for the bth binary trait of the ith individual, terms in (9) that involve y i,b only are extracted, such that where Here, ELSE refers to data and to the values of all parameters that the conditional distribution of the parameter of interest (y i,b ) depends on. Because the vector y i includes both liabilities and observations on continuous traits for the ith individual, it can be partitioned as where y i,-b represents y i but excluding the liability y i,b . Similarly, X i , Z i , W i , and K are partitioned conformably as Likewise, the residual covariance matrix R 0 is partitioned into a component pertaining to the bth binary trait (r b,b ), vectors containing the covariance components between the bth trait and all other traits (r -b,b and r b,-b ), and the residual covariance matrix of remaining traits (R -b,-b ), as follows: By properties of multivariate Gaussian distributions, the fully conditional posterior distribution of liability y i,b is where Because I i,b indicates whether the liability falls below or above the threshold, (11) represents the density of a normal distribution truncated at j b .

Location parameters
The joint conditional posterior distribution of location parameters is This expression can be recognized as the posterior density of the location parameters in a Gaussian-linear model with proper priors and known dispersion components [21], such that the corresponding distribution is Bayesian structural equation Gaussian-threshold model whereb û c 2 6 6 4 and yÃ ¼ Ky 1 ð Þ 0 Ky 2 ð Þ 0 ::: Ky n ð Þ 0 ð Þ 0 is a pseudo-data vector.

Structural coefficients and dispersion parameters
The fully conditional distribution of k can be derived following Gianola and Sorensen [10] and Wu et al. [26]. Because it does not have a recognizable form, a Metropolis-Hastings algorithm is used to sample k, centering the proposal at their current values [26]. In recursive models (i.e., K is an upper-or lower-diagonal matrix), K j j ¼ 1. Thus, the fully conditional distribution of k reduces to a multivariate normal distribution, and a Gibbs sampler can be used to sample k.
The conditional posterior distribution of the genetic covariance matrix G 0 is inverse Wishart [10]. The fully conditional posterior distribution of the matrix D 0 takes a form similar to that of the genetic covariance matrix.
When there are binary characters, because the variance of the liabilities of each binary character is fixed at 1, the residual covariance matrix R 0 is sampled from a conditional inverse Wishart distribution [14].

Ordered categorical traits
For an ordered categorical character there are two or more thresholds. If the first threshold is fixed, the other(s) have to be estimated. Note that h ¼ j; k; b; u; c; G 0 ; D 0 ; R 0 f g , where j is a vector containing all unknown thresholds. The joint posterior distribution p h; y b jg; y c ; H ð Þremains proportional to (9) if a uniform prior distribution is assigned to j. Thus, all unknown parameters are treated the same as for the case of binary characters, but an extra step is required to sample unknown thresholds during the MCMC steps. The fully conditional posterior distributions of the thresholds are independent, each of which is the collection of all relevant terms in (9). For example, consider the kth threshold for the jth categorical trait. It appears in connection with liabilities corresponding to responses in either the kth category (where the threshold is an upper bound) or the (k + 1)th category (where the threshold is a lower bound). This leads to the use of a uniform process to sample unknown thresholds [21].

Markov chain Monte Carlo sampling
Bayesian analysis via an MCMC implementation is used to infer marginal posterior distributions for parameters of interest. The MCMC sampling procedure consists of iterating through the following loop, after initializing parameters: 1a. Sample liabilities in y b ; 1b. Sample thresholds in j; 2. Sample structural parameters in k, using either the Metropolis-Hastings algorithm or a Gibbs sampler, and then update the ''data'' y Ã i ¼ Ky i ; 3. Sample location parameters in b, u, and c; 4. Sample the genetic covariance matrix G 0 ; 5. Sample the permanent environmental covariance D 0 ; 6. Sample the residual covariance matrix R 0 .
Step 1b is required only when ordered categorical characters are involved.

Transformation from liability to observable scale
In the recursive Gaussian-threshold model, the recursive effects from the categorical character (e.g., disease) to the Gaussian trait (e.g., production) are inferred on the underlying scale (i.e., liability to mastitis). To make interpretation easier these effects should be converted to the observable scale. A straightforward approach for conversion is the one of ''inverse probability'' [7,25]. Here, we present an intuitive approach that measures the difference in means of continuous traits (e.g., MY) between the two categories of a binary trait (e.g., mastitic and healthy), given the realization of underlying liabilities. Denote Here, y Ã i represents adjusted production for individual i (adjusted for all ''fixed'' and random effects, except liability to the disease, l i ), and e i is the Bayesian structural equation Gaussian-threshold model residual term. Then, the difference between means of production between sick (1) and healthy (0) cows can be calculated as where " l 1 and " l 0 are averages of augmented liabilities for sick and healthy cows, respectively, during the MCMC sampling.
2.8. Application to data from Norwegian Red cows The data represented 20 264 first-lactation daughters of 245 Norwegian Red sires that had their first progeny test in 1991 and 1992, and included test-day records for MY and veterinary records on clinical mastitis (CM) cases. Only test-day records from 5 to 180 days after calving were included. Cows with missing test-day records were excluded from the analysis for simplicity.
The 180 days of lactation were divided arbitrarily into three approximately equal-length periods: from day 5 to 60 (period 1), from day 61 to 120 (period 2), and from day 121 to 180 (period 3). For each period, cows were assigned the single MY test-day record that was closest in time to the mid-point of that period. For each test-day, a dummy variable indicating the presence or absence of CM in the 15-day period prior to the test-day was created. According to this definition of CM, a preexisting CM status would affect the following test-day MY, but the reverse would not occur.
Test-day MY decreased monotonically over the three lactation periods. The mean (standard deviation) of test-day MY was 21.40 (4.12) kg, 20.95 (4.02) kg, and 19.99 (4.00) kg at periods 1, 2, and 3, respectively. The presence or absence of CM was scored based on whether or not the cow had a CM treatment in a 15-day period prior to the test-day: 1 if a cow was treated for mastitis in the period and 0 otherwise. The incidence of CM decreased, approximately, from 3.0% at the first period to 0.9% at the second and third periods.

Model specifications
The data were analyzed using a standard multivariate linear sire model (LM), a recursive multivariate linear sire model (R-LM), a standard Gaussian-threshold (GT) sire model, and a recursive Gaussian-threshold (R-GT) sire model. For all models, it was assumed that correlations existed between sire effects as well as between residual effects, and that age at first calving (AGE) and herd affected all traits. AGE (''fixed'' effect) consisted of 15 classes with AGE < 20 months as the first class, AGE > 32 months as the last class, and each month in-between representing a single class. Herds, with 4903 classes, were treated as a random effect in the models, with herd effects affecting MY assumed to be uncorrelated with those affecting CM/liability to CM (LCM). In models R-LM and R-GT, the recursive effects were defined in a lagged manner, such that: CM1/LCM1 ! MY1 ! CM2/LCM2 ! MY2 ! CM3/LCM3 ! MY3, where the number following CM, LCM, and MY indicates the lactation period, and the arrow ! represents a causal relationship.

Analysis of posterior samples
The analyses were carried out using the SirBayes package (version 1.0), which is freely available upon request to the senior author (nickwu@ansci.wisc.edu). A detailed description of the convergence analysis can be found in Wu et al. [26]. Based on the convergence diagnostics results, it was decided that a single chain of 100 000 iterations would be used. Posterior samples from each chain were thinned every 10 iterations after 1000 iterations of burn-in. Genetic parameters were calculated for each thinned sample and saved simultaneously with posterior samples of location and dispersion parameters.

Recursive effects
In the three lactation periods, all recursive effects from LCM/CM to MY had negative posterior means, and those from MY to LCM/CM had positive means (Tab. I). These results suggest that an increased incidence of (or liability to) CM decreased MY at the following test-day, and that the effect from test-day MY to CM/LCM in the next lactation period would be weak. In model R-LM, all recursive effects from CM to MY were considered significant, because their 95% credible intervals did not overlap with zero. In model R-GT, however, only the recursive effect from LCM1 to MY1 could be considered significant, because the 95% credible intervals for the other two recursive effects included zero (Fig. 1a). Estimated recursive effects from CM/LCM to MY showed time-dependent patterns based on both models: the effect was the strongest in the first lactation period, and was reduced substantially in lactation periods 2 and 3. Based on model R-LM, for example, the recursive effect from CM to MY decreased from À0.33 kg per day in lactation period 1 to À0.11 and À0.13 kg per day in lactation periods 2 and 3. A similar trend was observed for the recursive effect from LCM to MY using model R-GT, as illustrated in Figure 1a.
An increase of one unit of LCM in model R-GT, which is equal to 1 residual standard deviation of liability, decreased test-day MY by 0.023 kg per day in lactation period 1, and by À0.002 kg to À0.004 kg per day in lactation periods 2 and 3. An increase in MY resulted in a non-significant increase in liability to CM in the following lactation period (Fig. 1b). The posterior mean of the effects from MY to LCM was between 0.001 and 0.002 liability units (Tab. I). These recursive effects obtained from model R-GT were converted to the observable scale, and the difference in mean test-day MY between the mastitic and healthy cows was À0.20 kg, À0.06 kg, and À0.09 kg per day, respectively, in lactation periods 1, 2, and 3. Converted recursive effects from LCM to MY based on model R-GT were smaller in absolute value than their counterparts based on model R-LM (i.e., from À0.11 kg to À0.33 kg per day), but both results pointed to the same direction, and they indicated the same pattern of influence.

Heritability
The presence of recursive effects in the models did not influence point or interval estimates of heritability for MY, and these estimates were also similar when using the linear or the Gaussian-threshold models (Tab. II). The posterior   heritability estimates for LCM obtained from the two Gaussian-threshold models were considerably higher than those for CM from the two linear models. The posterior means of heritability of LCM were from 0.07 to 0.09, whereas their counterparts of CM varied from 0.03 to 0.05 (Tab. II). The posterior distributions of heritability of MY and LCM were unimodal and approximately symmetric, as illustrated for the R-GT model in Figure 2. The posterior means of herd variances were from 0.25 to 0.36 for MY, and from 0.07 to 0.09 for LCM.

Genetic and residual correlations
Genetic correlations between test-day MY were in good agreement between the linear models and the Gaussian-threshold models (Tab. III). In general, test-day The results were obtained from the Gaussian-threshold models with (solid-lines) or without (dotted lines) the presence of recursive effects.
Bayesian structural equation Gaussian-threshold model MY in the three lactation periods were highly correlated; the posterior means of genetic correlations between test-day MYs ranged from 0.79 to 0.98. As expected, the closer the test-days were, the higher the correlation between MYs was. Posterior distributions of genetic correlations between test-day MYs were highly overlapping between models with or without recursive effects (Fig. 3).  Genetic correlations between LCM in the three lactation periods were also high with the posterior means from the R-GT model ranging from 0.83 to 0.93 (Tab. III). Genetic correlations between CM obtained from the recursive linear models were smaller, ranging from 0.4 to 0.5 (Tab. III). Nevertheless, posterior distributions of genetic correlations between LCM or CM were highly overlapping for the models with or without the presence of recursive effects. Posterior distributions of genetic correlations between LCM from model R-GT are shown in Figure 4.
Genetic correlations between LCM and MY ranged from 0.14 (LCM2-MY3) to 0.62 (LCM1-MY1), based on model R-GT. Genetic correlations between CM and MY were slightly smaller, ranging from 0.08 (CM2-MY3) to 0.37  (CM1-MY1) based on model R-LM. The presence of recursive effects changed only slightly the posterior distributions of genetic correlations between LCM and MY (Fig. 5). The same was true for the two linear models. Residual correlations are given in Table IV. Differences in residual correlations between models with or without recursive effects were minor. The posterior mean of residual correlations between test-day MYs was comparable between the two types of models, ranging from 0.  Many diseases are measured as categorical, rather than quantitative, traits and often as binary response variables. The vast majority of these disease traits have a polygenic basis. Thus, Wright [25] proposed a ''physiological threshold'' theory to explain the link between a continuous latent variable, also referred to as ''liability'' [8], and an observable binary phenotype. The basic assumptions of a threshold model are that: (1) there is an underlying variable whose value is the sum of a normally distributed environmental component and an independent normally distributed genetic component; (2) the ''affected'' character is present in only those, in which the underlying variable exceeds a certain threshold value; and (3) gene substitutions have individually small and strictly additive effects on the underlying variable. A model in which additive action is at the level of some underlying variable (liability scale) may be more sensible than one based on additive gene action on the outward variate (probability scale). Concerning heritability of a binary character, for example, the use of a liability scale circumvents problems arising in the probability scale [7]. Biologically, the linear model assumes a dichotomous (0 or 1) influence of CM whereas, by assuming a threshold model, the effect of LCM on MY is continuous. Mastitis is a complex trait that can be caused by many different pathogens and shows different infection patterns, from mild to very severe clinical cases. It is thus more reasonable to believe that MY is more affected by a severe mastitis than by a mild clinical case. For cows observed as healthy (i.e., CM = 0) there may also be variation in effects on MY, as their health status may vary from completely healthy to almost mastitic (i.e., subclinical mastitis). Therefore, a Gaussian-threshold model is more preferable than a multivariate linear model to describe the relationship between CM and MY.
In this paper, the model of Gianola and Sorensen [10] was extended to describe relationships between Gaussian traits and liabilities of binary traits. Recursive effects from LCM to MY were estimated on the underlying scale of disease (i.e., liability to mastitis) and converted to the probability scale. The conversion method measures the difference in mean MY between mastitic and healthy cows, given the realized liabilities. Converted recursive effects from the Gaussian-threshold model were smaller than those obtained from the linear model. This probably reflects intrinsic differences between these two types of models, i.e., the linear models produce frequency-dependent inferences [9]. Further, in the MCMC sampling, the chains for recursive effects in the linear model did not mix as well as they did in the Gaussian-threshold model, when the incidence of disease was low. Because CM and LCM have different heritabilities, it is also possible that recursive effects between CM and MY based on model R-LM are different from those between LCM and MY based on model R-GT. Nevertheless, both models led to the same conclusion, since recursive effects from both models showed the same pattern and were in the same direction.
At the phenotypic level, the influence of CM on production (e.g., MY) has been documented e.g. [19,22]. There is also evidence that high MY may increase incidence of CM e.g. [16,17]. Using a SIR model, Wu et al. [26] found positive recursive effects from MY to somatic cell score (SCS), and decreasing effects from SCS to MY as lactation proceeded. These results may reflect a relationship between CM incidence and the magnitude of the recursive effect to MY. Incidence of CM decreases as lactation progresses [1,3], and so did the recursive effect from CM to the following test-day MY, possibly because of reduced variation. Both the recursive linear model and the recursive Gaussian-threshold model showed negative effects from (liability to) CM to MY, and these effects decreased from lactation period 1 to later periods.

Estimation of genetic parameters under recursive relationships
The presence of recursive effects in the models did not affect point or interval estimates of heritability of LCM or MY. This conclusion was in agreement with previous studies, in which the estimates of heritability obtained assuming SIR relationships [6,26] were similar to those from standard models. Estimated heritabilities for test-day MY were in agreement with previous reports in the same population [6]; estimated heritabilities for liability to CM were in agreement with those of Chang et al. [5] for the same lactation periods, and slightly higher 354 X.-L. Wu et al. than those of Heringstad et al. [11], who estimated heritability of CM in the course of lactation in the same population using a longitudinal threshold model.
Wu et al. [26] found that estimates of some genetic and residual correlations from SIR models could differ considerably from those obtained using standard mixed models. In the present analysis, however, similar estimates of genetic and residual correlations were found regardless of the presence of recursive effects in the models. The observed difference in genetic correlations could be data driven. However, there was a discrepancy between the linear models and the Gaussianthreshold models in estimates of genetic and residual correlations involving (liability to) CM.
Genetic correlations between (liabilities to) CM in the three lactation periods were similar to some previous reports [5,11]. The positive, moderate to high, genetic correlations between LCM/CM and MY were in agreement with previous studies [4,12], and indicate the involvement of common genetic factors or pathways in genetic expression and regulations of these two traits [13,20]. From the viewpoint of genetic selection, the positive genetic correlations between liability to mastitis and MY are unfavorable, because selection for higher MY would be associated with an increased liability to CM. It is known that there is an antagonistic genetic correlation between mastitis and milk production e.g. [4,12] but knowledge is limited regarding how this association evolves in the course of lactation. Thus, the present application represents an effort toward obtaining a dynamic picture of these relationships.