J. Anim. Sci. 2004. 82:3465-3473
© 2004 American Society of Animal Science
Use of structured antedependence models for the genetic analysis of growth curves1
F. Jaffrézic
2,
E. Venot,
D. Laloë,
A. Vinet and
G. Renand
INRA Quantitative and Applied Genetics, 78352 Jouy-en-Josas Cedex, France
 |
Abstract
|
|---|
Growth curve analysis is an important issue for many agricultural and laboratory species, for both phenotypic and genetic studies. The aim of this paper is to present the use of a novel statistical approach, namely the structured antedependence (SAD) models, to deal with this issue. The basic idea of these models is that an observation at time t can be explained by the previous observations. These models are especially appropriate to deal with cumulative traits such as growth, as BW at age t clearly depends on BW measures at ages (t –1), (t –2), etc. These models were applied on an INRA experimental Charolais herd data set. The data comprised BW records for 560 cows born over an 11-yr period (from 1988 to 1998) from 60 sires and 369 dams. The proposed SAD models were compared with the well-known random regression (RR) models that are already widely used in various areas of longitudinal data analysis. It was found that the SAD models fit the growth process better with far fewer parameters than the RR models (9 instead of 16 covariance parameters for the phenotypic analysis, and 14 instead of 21 for the genetic analysis). Despite this smaller number of covariance parameters, the likelihood value was found to be much higher with the SAD vs. the RR models, with a difference of 262.9 for the phenotypic analysis with a quartic polynomial for the RR and 751.5 for the genetic analysis with a cubic polynomial for both the genetic and environmental parts of the RR model. The SAD models also proved to be better able to interpolate missing values. Heritability, genetic, and environmental correlation coefficients were estimated for weights from birth to adulthood. The structured antedependence models proved, in this study, to be very appropriate to model growth data in a parsimonious and flexible way.
Key Words: Beef Cattle • Genetic Analysis • Growth Curves • Random Regression Models • Structured Antedependence Models
|
Introduction
|
|---|
Growth curve analysis is an important issue for many agricultural species such as beef cattle, pigs, chicken, and rabbit. In the case of genetic studies, a few methodologies have been proposed. Blasco et al. (2003) proposed using parametric growth curves, such as the Gompertz curve, that have interpretable parameters in terms of growth rate and adult BW. Each parameter of the curve is then decomposed into genetic and environmental components, as in classical quantitative genetics theory, and selection can then be made on these parameters. This approach, however, leads to nonlinear mixed models, and therefore requires the use of Markov chain Monte Carlo methods to estimate the parameters, which can be quite time consuming, especially for large data files.
On the other hand, Meyer (2002, 2004) proposed using random regression (RR) models (Diggle et al., 1994), as in the classical longitudinal data analyses. This approach has already been extensively investigated in the context of test-day records in dairy cattle and proved to have a few drawbacks. In particular, polynomial functions tend to show border effects that lead to the prediction of unreasonable values at the extreme points of the data. A possible way to overcome this drawback would be to use spline functions instead of polynomials; they do, however, still require a quite large number of parameters.
Other approaches have been proposed in the context of longitudinal data analyses, such as structured ante-dependence (SAD) models (Nunez-Anton and Zimmerman, 2000). These models are especially appropriate for cumulative traits, as the observation at time t is defined as a function of the previous observations.
The aim of this paper was to investigate the use of these SAD models for growth curve analysis and apply them to the genetic analysis of an experimental beef cattle data set.
|
Materials and Methods
|
|---|
Structured Antedependence Models
Let Yi(tj) be the observed body weight at age tj for animal i. It is assumed that it can be decomposed into genetic and environmental components, as follows:
 | [1] |
where µ(tj) corresponds to the fixed effects and includes the mean curve in the population. Functions gi(tj) and ei(tj) are Gaussian processes, which are independent of one another, have an expected value of 0 at each age, and covariance functions G(tj,tk) and E(tj,tk), respectively, for two ages tj and tk. They represent the age-dependent genetic and environmental deviations, respectively, for animal i. It is assumed that J measurement times are available with t1 < t2 < ... < tJ. These ages are assumed to be the same for all the individuals; however, the model is able to deal with missing values.
The idea of antedependence models, as originally proposed by Gabriel (1962), is that an observation at time tj can be explained by the previous observations. An antedependence structure of order r is defined by the fact that the kth observation (k > r) given the r preceding ones is independent of all further observations (Gabriel, 1962). Generalizing this concept to genetic analysis, a second-order structured antedependence model for the genetic part gi(tj), for animal i, can be written as:
 | [2] |
 | [3] |
 | [4] |
for j 
2.
Parameters 
1j and 2j are termed antedependence parameters. As originally proposed by Gabriel (1962), these parameters were estimated for each time j. Therefore, if J measurement times were considered, (J –1) parameters 1j and (J –2) parameters 2j had to be estimated. If the number of measurement times is large and the order of antedependence high, that will lead to a very large number of parameters to be estimated. Therefore, Nunez-Anton and Zimmerman (2000) proposed in their SAD models to decrease this number of parameters by fitting a parametric function of the lag time for these antedependence terms. For example, when using an exponential function, for j = 1, ..., J: 1j (tj, tj–1) = exp(–
1 * [tj– tj– 1]) and 2j (tj, tj–2) = exp(–2 *[tj–tj–2]
A further advantage of modeling the antedependence parameters with a continuous function of the lag time is that it allows dealing with unequally spaced data, as well as unbalanced data. It is therefore not required that the animals be weighed at regular time intervals or that the times of measurement are the same for all the animals. To be able to find an appropriate parametric curve, the time spacing should not be completely random. The model can, for example, easily deal with smaller time intervals at earlier rather than late ages.
Although antedependence models are based on an idea similar to autoregressive models, due to the initial condition gi(t0) = 
gi (t0), there is no constraint on the antedependence parameters 1j and 2j, whereas they have to be between –1 and 1 for autoregressive models. Therefore, any parametric function can be used to model these parameters, depending on the biological process studied.
Parameters gi (tj) are assumed normally distributed, with mean 0 and variance 
2g(tj) that are termed "innovation variances." Gabriel (1962) originally suggested estimating one innovation variance for each time of measurement. This again leads to a quite large number of parameters when many different times of measurement are analyzed. In their structured antedependence models, Nunez-Anton and Zimmerman (2000) again suggested using a continuous function of time to model these innovation variances, such as polynomials. For a quadratic polynomial, for instance:
 | [5] |
In this case, therefore, only three parameters will have to be estimated, regardless of the number of times of measurement J. This parametric modeling of the error variances can also be related to the structural model proposed by Foulley and Quaas (1995). The use of a polynomial function to model the variances might, as with the RR models, show some border effects. These border effects will, however, be far less important for the SAD models than for the RR case as the estimated covariance function of the process is not a simple function of these polynomials but rather a complex combination of the antedependence parameters and innovation variances. In practice, simple parametric functions are sufficient to adequately model these variances. When possible, it is best to fit an unstructured antedependence model (UAD) as originally proposed by Gabriel (1962), which will give a nonparametric estimation of the antedependence parameters and innovation variances and will allow one to choose the most appropriate parametric functions of time.
To relate this modeling to the genetic covariance function G(tj,tk), for two ages tj and tk, Pourahmadi (1999) showed that a Cholesky decomposition of the inverse of the covariance matrix can be easily calculated using the antedependence parameters and innovation variances. If there are J measurement times, let G be the genetic covariance matrix, of dimension J x J. It can be shown that:
 | [6] |
where L is a lower triangular matrix with 1 on the diagonal and negatives of the antedependence coefficients rj (r = 1, ..., R for SAD(R), j = 1, ..., J) below the diagonal entries, and D is a diagonal matrix with innovation variances 2j (j = 1, ..., J) as components. An interesting computational property is that the inverse G–1is sparse. Indeed, for a second-order antedependence model, for instance, only the first two subdiagonals are nonzero. Antedependence and innovation variance parameters can be estimated by restricted maximum likelihood procedures. In the following example, these parameters were estimated using the OWN function of ASREML (Gilmour et al., 2002), which requires specification of the covariance matrix and derivatives with respect to each of the parameters. As presently implemented, it is necessary to build the J x J covariance matrix, which will be the same for all the animals. The J times of measurement are therefore assumed to be the same for all the individuals, and should not be too large (20 at most). It is not required, however, that these ages should be equally spaced or that the animals have measurements at all times.
Presentation of the Data
Data analyzed in this study came from an INRA experimental Charolais herd (Mialon et al., 2001). The data set comprised BW records for 560 cows born over an 11-yr period (from 1988 to 1998) from 60 sires and 369 dams. Data were collected monthly from 1998 to 2003, and only measurements at 10 different ages were considered here for each animal, at approximately 0, 112, 224, 364, 540, 720, 900, 1,260, 1,620 and 1,980 d. Although the same ages were considered for each animal, they were unequally spaced and some records were missing. The fixed effects were the year of birth of the animal, twinning effect and the dam age at calving. To obtain the most accurate fit, the mean curve in the population was modeled nonparametrically by fitting one mean at each age by including in the fixed effects µ (tj) of Eq. [1] the variable age as a factor. This was possible because only 10 measurements were considered for each animal and the data had very few missing values. In other practical cases, however, it could be better to use parametric or semiparametric functions of age, such as spline functions.
Analysis
First, a phenotypic analysis of the data was performed to compare various different SAD and RR models. Here, models with antedependence up to order 3 were investigated. The innovation variances were assumed to be changing as a quadratic function of age, and the antedependence parameters were either considered as constant, linear, quadratic, or exponential functions of age. Random regression models based on polynomials up to order 4 (quartic) were used.
Models were compared using the likelihood values and the Bayesian information criterion (BIC; Schwarz, 1978): BIC = ln L –0.5 nc ln(N –p), where ln L is the restricted maximum likelihood value, nc is the number of covariance parameters in the model, p is the number of fixed effects (also equal to rank (X)), and N is the total number of observations. The number of fixed effects in the model was equal to 25 and the total number of observations was 5,455.
This analysis was followed by a genetic study. An animal model was used and the pedigree file comprised 807 animals. Several SAD and RR models were compared for both the genetic and permanent environmental parts. As before, models were compared using the likelihood values and the BIC criterion. The estimated genetic parameters obtained at different ages for the chosen models are presented in the next section.
|
Results
|
|---|
Phenotypic Analysis
As shown in Table 1, it was found that even the simplest first-order SAD model (with only four parameters for the covariance structure) provided a better likelihood value than the quartic RR model (with 16 parameters).
View this table:
[in this window]
[in a new window]
|
Table 1. Likelihood values (Log L) and Bayesian information criterion (BIC) for the phenotypic analysis using different structured antedependence (SAD) models, with antedependence up to order 3, and with the antedependence parameters modeled as constant (const), linear (lin), quadratic (quad), or exponential (exp) functions of time, different random regression (RR) models based on polynomials up to order 4 (quartic), and a completely unstructured model (US) with a 10 x 10 covariance matrix
|
|
The SAD model chosen in this phenotypic study, based on the BIC criterion, is a third-order antedependence model that can be written for any animal i at age t as follows:
 | [7] |
where the antedependence parameters and innovation variances were modeled as follows:
The 10 different ages considered here were coded as t = 1, ..., 10.
The choice of these polynomial functions to model the antedependence parameters was based on a preliminary analysis using a first-order unstructured antedependence model. Figure 1 presents the nonparametric estimates for the 1t antedependence parameter using the first-order unstructured antedependence model and estimates obtained with a first-order SAD model using a quadratic polynomial. It can be seen that the quadratic function is appropriate to model this first antedependence parameter.
View larger version (16K):
[in this window]
[in a new window]
|
Figure 1. Estimated first-order antedependence parameter obtained nonparametrically with the unstructured antedependence model (UAD(1)) and with a quadratic polynomial function with the structured antedependence model (SAD(1)).
|
|
As pointed out in the Materials and Methods section, it might also be possible in other practical examples to use nonlinear functions, such as exponential functions that have asymptotes. Estimation of these nonlinear terms can easily be implemented in ASREML (Gilmour et al., 2002) using the OWN function; it simply requires specifying the derivatives of the chosen function with respect to each of the parameters.
The antedependence model presented in Eq. [7] required only nine parameters for the covariance structure and gave a likelihood value of 1,615.7, whereas the most complicated RR model required 16 parameters, with a likelihood value equal to 1,352.8.
Figure 2 gives the correlation functions estimated with the chosen SAD and RR models. To evaluate the goodness-of-fit of these functions, an unstructured covariance matrix (10 x 10) was also fitted to the data. The estimated correlation function obtained with this saturated model is also presented in Figure 2. One of the main drawbacks of the RR models based on polynomial functions can clearly be seen here, namely border effects (Druet et al., 2003). In fact, although the global shape of the correlation function seems to be close to the unstructured one, a strongly negative correlation was estimated at early ages, which seems to have no biological meaning and is not observed with either the SAD or unstructured model.
View larger version (46K):
[in this window]
[in a new window]
|
Figure 2. Estimated phenotypic correlation functions obtained with the unstructured (US) model and the structured antedependent (SAD) and random regression (RR) models chosen according to the Bayesian information criterion (BIC) values given in Table 1 .
|
|
Conversely, Figure 3 gives the estimated variance functions obtained with the different models, and they all seem quite similar. As expected, the variances were found to quite strongly increase with age (and BW).
View larger version (17K):
[in this window]
[in a new window]
|
Figure 3. Estimated phenotypic variance functions obtained with the unstructured (US) model and the structured antedependent (SAD) and random regression (RR) models chosen according to the Bayesian information criterion (BIC) values given in Table 1.
|
|
It was found that both the chosen SAD and RR models were able to fit individual phenotypic curves well when all the observations were available, even for irregular patterns. The main difficulty of the RR model seemed to be the prediction of missing observations at early ages (due to the poor estimation of the correlation function at these ages). Figure 4 shows, for example, the interpolations obtained with the different models for an animal missing values at ages 112 and 224 d. Both the SAD and the unstructured model predicted the same values, whereas the RR model predicted lower BW values (the difference was approximately equal to 20 kg at age 112 d). The accuracy of the covariance structure is therefore particularly important for dealing with missing values.
To check the interpolation abilities of both models at early ages, observations at age 112 d were deleted from the data for all the animals born since 1995, which represented 220 animals and approximately 40% of the population analyzed. To compare the predictions obtained with the SAD and RR models to the actual BW values at age 112 d for these animals, the Vonesh concordance coefficient (rc) was used (Vonesh et al., 1996) with the following formula:
 | [8] |
where M is the number of animals with deleted values (here M = 220), yi represents the actual observation for animal i at age 112 d, 
i is the predicted value obtained either with the SAD or RR model, and &ymacr; and ŷ are the averages of the observed and predicted values, respectively. This rc coefficient has values between –1 and 1, with a perfect fit at 1 and a lack of fit for negative values.
When calculating this coefficient for the interpolation of the 40% missing values at age 112 d, it was found that the prediction ability was much better with the SAD than with the RR model. In fact, the concordance coefficient was found to be equal to 0.83 for the SAD model for the predicted vs. actual observed BW values, and was 0.44 for the RR model.
The differences in the interpolation abilities were less important for other ages; for example, rc = 0.69 at 224 d of age for the SAD model compared with 0.58 for the RR model, or rc = 0.74 at 540 d of age for the SAD model compared with 0.70 for the RR model.
Finally, the rc coefficient was equal to 0.74 with the SAD model compared with 0.69 with the RR model for the latest age (1,980 d), showing a slightly better extrapolation ability of the SAD model. Figure 5 illustrates this prediction ability for an animal with a "classical" growth curve. In this particular example, the phenotypic prediction obtained with the RR model was decreasing, whereas the BW for this animal would be expected to remain stable, as predicted by the SAD model, or to keep increasing slightly, as actually observed.
View larger version (20K):
[in this window]
[in a new window]
|
Figure 5. Observed and extrapolated phenotypic growth curves obtained with the chosen structured ante-dependent (SAD) and random regression (RR) models for an animal with observations for age greater than 1,620 d deleted.
|
|
Genetic Analysis
Model comparison results for the genetic analysis are presented in Table 2. As for the phenotypic analysis, it was found that the simplest first-order SAD model, with a constant antedependence parameter, had a higher likelihood value than a cubic RR model, while requiring far fewer parameters for the covariance structure (eight instead of 21 for the RR model). The best-fitting SAD model regarding the likelihood value and the BIC criterion was a second-order antedependence, with a quadratic first-order antedependence parameter, a constant second-order antedependence parameter, and quadratic innovation variances for both the genetic and environmental parts.
For the genetic part, the chosen model was as follows (for any animal i, at age t):
 | [9] |
with parameter estimates equal to 1t = 1.29–0.11 t + 0.01 t2, 2 = 0.34, and Var(et) = –15.8 + 33.8 t –1.1 t2.
For the environmental part, the chosen model was as follows:
|
| [10] |
where 1t = 1.50–0.24 t + 0.01 t2, 2 = 0.05, and Var(eit) = –157.4 + 174.4 t –3.5 t2.
Figure 6 shows the estimated genetic correlations obtained with this SAD model and the cubic RR model. The patterns were quite similar, although a more important drop was observed for the RR model at early ages. These genetic correlations were equal to more than 0.95 at late ages and approximately 0.4 at early ages. On the other hand, for the environmental correlations the same drawback was observed for the RR model as in the phenotypic analysis. In fact, as shown in Figure 7, the correlations even became negative at very early ages, although there seems to be no biological reason for this. As expected, the environmental correlations obtained with both methods were found to be much lower than the genetic ones. They ranged from about 0.1 for most ages, to around 0.6 for a few adjacent ages. Due to the very large number of parameters to be estimated, the completely unstructured model (with 110 parameters) did not properly reach convergence in this genetic study, and therefore no "reference" model is available.
View larger version (34K):
[in this window]
[in a new window]
|
Figure 6. Estimated genetic correlation functions obtained with the structured antedependent (SAD) and random regression (RR) models chosen according to the Bayesian information criterion (BIC) values given in Table 2.
|
|
View larger version (34K):
[in this window]
[in a new window]
|
Figure 7. Estimated environmental correlation functions obtained with the structured antedependent (SAD) and random regression (RR) models chosen according to the Bayesian information criterion (BIC) values given in Table 2.
|
|
Figure 8 gives the estimation of the heritability over age for the chosen SAD model. The pattern is characteristic of heritability estimates found in the literature. In a review of univariate estimates of heritability coefficients, Koots et al. (1994) computed average values of 0.35, 0.27, 0.35, and 0.50 for birth weight, weaning weight, yearling weight, and adult cow weight, respectively. A similar pattern was obtained by Bullock et al. (1993) with field data records of Polled Hereford cows: 0.49, 0.24, 0.30, and 0.52, respectively. In two studies of growth curves up to 2 yr of age, Meyer (2002, 2004) also obtained a similar pattern for heritability coefficient estimates, with RR models. In the current study, heritability of birth and adult cow weights (0.54 and 0.70) were among the highest estimates found in the literature. Figure 8 also presents the heritability estimates obtained with the cubic RR model; it was also found to be increasing with age, although the heritability estimated for the adult body weight was slightly lower than with the SAD model (equal to 0.6). The main difference in the pattern of these two heritability curves is for the birth weight. The RR model estimated a very low heritability at birth, equal to 0.05, whereas it is expected to be much higher as previously found in the literature (Bullock et al., 1993; Koots et al., 1994).
View larger version (14K):
[in this window]
[in a new window]
|
Figure 8. Estimated heritability functions over age obtained with the chosen structured antedependent (SAD) and random regression (RR) models, as presented in Table 2.
|
|
|
Discussion
|
|---|
Use of longitudinal models over the simple multitrait model for growth curve analysis has many advantages. In particular, they are able to capture the continuous nature of the process, and therefore allow prediction of the genetic values at any time, even when no measurement is available. They also allow the analysis of unbalanced data and irregularly spaced measurements.
Among them, the SAD models in this study proved to be appropriate to model growth data. They, in fact, offer a parsimonious and flexible way of modeling this process, and were found to be able to better fit the covariance structure of the data than polynomial-based random regression models. It was found in the phenotypic study that even the simplest first-order SAD model (with only four parameters) provided a better fit (higher likelihood value) than a fourth-order RR model (with 16 parameters). Similar results were obtained in the genetic study where a simple first-order SAD model, for both the genetic and environmental parts (with eight parameters), provided a better fit than a cubic RR model (with 21 parameters).
The chosen SAD and RR models proved able to fit individual phenotypic curves very well, when all the observations were available. Differences were noted between the two methodologies for interpolating the missing values. Due to the poor estimation of the correlation structure obtained with the RR models, they generally predicted the missing values less accurately than the SAD models.
To further improve the modeling of the covariance structure, it would also be possible to model variance heterogeneity due to other environmental factors, as proposed by Foulley and Quaas (1995), and Robert-Granié et al. (2002) for simple RR models. This is straightforward in the SAD models, as a structural model is already used to model the changes of the innovation variances over time. Extension of this heterogeneous model can also be applied to the antedependence parameters, to correct, for example, for preferential treatments influencing the growth rate of some animals.
A multivariate extension of the SAD models has been presented by Jaffrézic et al. (2003), which would allow one to study, for example, the dependence between BW and feed intake, or any other time-dependent trait of economic interest.
|
Footnotes
|
|---|
1 Thanks to W. G. Hill for very interesting comments and ideas. Data were provided by the INRA experimental center located near Bourges, France. 
2 Correspondence—e-mail: jaffrezic{at}dga2.jouy.inra.fr.
Received for publication June 17, 2004.
Accepted for publication September 7, 2004.
|
Literature Cited
|
|---|
Blasco, A., M. Piles, and L. Varona. 2003. A Bayesian analysis of the effect of selection for growth rate on growth curves in rabbits. Genet. Sel. Evol. 35:21–41.[Medline]
Bullock, K. D., J. K. Bertrand, and L. L. Benyshek. 1993. Genetic and environmental parameters for mature weight and other growth measures in Polled Hereford cattle. J. Anim. Sci. 71:1737–1741.[Abstract]
Diggle, P. J., K. Y. Liang, and S. L. Zeger. 1994. Analysis of Longitudinal Data. Oxford Univ. Press, Oxford, U.K.
Druet, T., F. Jaffrézic, D. Boichard, and V. Ducrocq. 2003. Modeling lactation curves and estimation of genetic parameters for first lactation test-day records of French Holstein cows. J. Dairy Sci. 86:2480–2490.[Abstract/Free Full Text]
Foulley, J. L., and R. L. Quaas. 1995. Heterogeneous variances in Gaussian linear mixed models. Genet. Sel. Evol. 27:211–228.
Gabriel, K. R. 1962. Ante-dependence analysis of an ordered set of variables. Ann. Math. Stat. 33:201–212.
Gilmour, A. R., B. J. Gogel, B. R. Cullis, S. J. Welham, and R. Thompson. 2002. ASREML User Guide Release 1.0. VSN Int., Hemel Hempstead, U.K.
Jaffrézic, F., R. Thompson, and W. G. Hill. 2003. Structured antedependence models for genetic analysis of multivariate repeated measures in quantitative traits. Genet. Res. 82:55–65.[Medline]
Koots, K. R., J. P. Gibson, C. Smith, and J. W. Wilton. 1994. Analyses of published genetic parameter estimates for beef production traits. 1. Heritability. Anim. Breed. Abstr. 62:309–338.
Meyer, K. 2002. Estimates of covariance functions for growth of Australian beef cattle from a large set of field data. CD-ROM communication No. 11-01 in Proc. 7th World Cong. Genet. Appl. Livest. Prod., Montpellier, France.
Meyer, K. 2004. Scope for a random regression model in genetic evaluation of beef cattle for growth. Livest. Prod. Sci. 86:69–83.
Mialon, M. M., G. Renand, D. Krauss, and F. Ménissier. 2001. Variability of the postpartum recovery of sexual activity of Charolais cows. Livest. Prod. Sci. 69:217–226.
Nunez-Anton, V., and D. L. Zimmerman. 2000. Modeling non-stationary longitudinal data. Biometrics 56:699–705.[Medline]
Pourahmadi, M. 1999. Joint mean-covariance models with applications to longitudinal data: Unconstrained parameterisation. Biometrika 86:677–690.[Abstract/Free Full Text]
Robert-Granié, C., B. Heude, and J. L. Foulley. 2002. Modelling the growth curve of Maine-Anjou beef cattle using heteroskedastic random coefficients models. Genet. Sel. Evol. 34:423–445.[Medline]
Schwarz, G. 1978. Estimating the dimension of a model. Ann. Stat. 6:461–464.
Vonesh, E. F., V. M. Chinchilli, and K. Pu. 1996. Goodness-of-fit in generalized nonlinear mixed-effects models. Biometrics 52:572–587.[Medline]
This article has been cited by other articles:
|
|
|
|
 
M. P. L. Calus, L. L. G. Janss, and R. F. Veerkamp
Genotype by Environment Interaction for Somatic Cell Score Across Bulk Milk Somatic Cell Count and Days in Milk
J Dairy Sci,
December 1, 2006;
89(12):
4846 - 4857.
[Abstract]
[Full Text]
[PDF]
|
|
|
Top
Abstract
Introduction
