Restricted maximum likelihood estimation of genetic parameters for the first three lactations in the Montbéliarde dairy cattle breed

Genetic parameters for the first three lactations have been estimated for the main dairy traits (milk, fat, protein and useful yields adjusted for lactation length, fat and protein contents). Two data sets were analysed, including records on 30 751 cows born from 128 young sires and 52 proven sires. Daughters' performances from the most widely used proven sires were incorporated in order to improve the degree of connectedness among herds. The model fitted young sires as random and proven sires, herd-year, season-year of calving, age at first calving and length of the previous lactation as fixed effects. Relationships among bulls were included. Analysis was by restricted maximum likelihood using an EM-related algorithm and a Cholesky transformation. All genetic correlations were larger than 0.89. Correlations between the first and third lactations were slightly lower than the others. Heritabilities of milk, fat, protein and useful yields ranged from 0.17 to 0.27. Phenotypic correlations between successive lactations were higher than 0.6 and those between lactations 1 and 3 lower than 0.55. Heritabilities of fat and protein contents were higher than 0.44 with phenotypic correlations being stable at about 0.70. The "repeatability model" which considers all lactation records as a single trait can be con


INTRODUCTION
The goal of dairy selection is to improve lifetime production of cows, which implies taking into account the different lactations. Until now, genetic evaluation of the animals has in most cases been made under the assumption that these lactations are influenced by the same genes. In some countries only the first lactations are considered; in others the so-called &dquo;repeatability model&dquo; (Henderson, 1987) in which all lactations are treated as repetitions of one trait is fitted. But the lactations are made at various ages and physiological status of the animals and may therefore be determined somewhat by different genes. The accuracy of the genetic evaluation and thus the effi!ciency of dairy selection might be improved by fitting a multitrait model to the lactations. Reliable estimates of the genetic parameters for the different lactations are needed to appreciate this possible gain in accuracy.
Data usually available for such estimations are selected as breeders cull about one quarter of the animals by the end of each lactation. Their decision is mostly based on dairy performance. Useful methods of estimation of these parameters have been available only recently. Henderson's methods (1953) assume animals are measured for all lactations, thus leading to results biased by the selection. However, the maximum likelihood (ML) (Hartley and Rao, 1967) and restricted maximum likelihood (REML) (Patterson and Thompson, 1971) estimators can take into account this selection (Im et al., 1987), a necessary condition being that the selection process is based only on the observed data or on observed data and independant variables. REML was prefered to ML as it accounts for the loss of degrees of freedom in simultaneous estimation of the fixed effects. Moreover, theoretical studies have shown that the optimum statistical procedure maximising the genetic merit of selected animals consists of estimating variance and covariance components by REML and thereafter applying these estimates in the mixed model equations (Gianola et al., 1986).

Data
Records for the first 3 lactations of Montb6liarde cows whose first calving occurred between 1/09/1979 and 30/08/1982 were extracted from the National Milk Recording files. The conditions of editing are presented in Table 1. Records made after cows changed herds were disregarded (Meyer, 1984). They represented 1.5% of the records for second lactation and 1.3% of the records for third lactation. Cows were nested within herds and the absorption matrix of the herd effects was block diagonal for herds.
Two populations of females were considered. The first was made of daughters of test bulls. It was used to estimate sire components of variance and covariance. The second consisted of daughters of the most widely used proven sires. As these bulls had been selected, they were treated as fixed effects and were not considered for the estimation of sire components. The performances of their daughters were introduced in the analysis in order to improve the accuracy of the estimation through additional information, increased herd size and degree of connectedness between herds.
A total of 180 bulls of which 128 were random test bulls was considered. To simplify computation, records were split into 2 data sets, as did Meyer (1984Meyer ( , 1985a and Swalve and Van Vleck (1987). The first data set consisted of the daughters of sampling bulls born in 1975 and the second of the daughters of sampling bulls born in 1976. For each of the 2 data sets, the most widely used proven sires were determined and their daughters added. Table II summarizes their main characteristics.
The main dairy variables were considered: milk, fat, protein and useful yields, and fat and protein contents. Useful yield (UY) is defined as: Yield traits were corrected multiplicatively for lactation length prior to analysis, according to Poutous and Mocquot (1975) as: Corrected yield = (total yield x 385)/(lactation length + 80) Data were scaled to reduce rounding errors.

Model
The following model was used for each of the 6 variables: where y is the vector of the observations; h is the vector of fixed herd effects (the number of levels of which is shown in Table II); b is the vector of fixed year-season of calving (15 levels for each lactation), age at 1st calving (10 levels for each lactation) and length of the preceding lactation effects (8 levels for each lactation); u is the vector of the sire effects (this effect was treated as fixed when the sire was a proven bull, and as random and normally distributed when the sire was a young bull); and e is the vector of residual effects, assumed normally distributed; X, W and Z are known incidence matrices for the herd effects, the other fixed effects and the sire effects.
Expectations and variances are defined as: where u i is the subvector of u corresponding to the effects of the young bulls and: where A is the relationship matrix of the young bulls, T the matrix of the sire components and * the right direct product (Graybill, 1983). Let n denote the number of animals; with data ordered by lactations within animals, R is block diagonal having n blocks R,! (k = 1, ... n). If the kth cow has made the first 3 lactations, R,! = E where E is the matrix of residual components; if it has been culled before, the rows and columns corresponding in E to the missing records are deleted.

Method
Data were Cholesky transformed (Schaeffer, 1986) but, because the incidence matrix W varied for an animal from one lactation to the next, the vector b could not be transformed and the mixed model equations for this parameter remained unchanged. A combined REML/ML procedure (Meyer, 1983a(Meyer, , 1984(Meyer, , 1985a was used. From a Bayesian viewpoint, only herd effects were integrated in the posterior density. From a classical viewpoint, the inference was based on the likelihood of (n &mdash; r(X )) error contrasts K'y (where K is an n x (n &mdash; r(X)) matrix such that K'X = 0 (Harville, 1977). An algorithm &dquo;related to the E.M.&dquo; (Henderson, 1985) was used. First derivatives of the restricted likelihood function are set to zero (eqn. (7) of Meyer (1986)). For the computation of residual components, use is made of eqn.
(8) of Meyer (1986). For computations, the iterations were stopped when the relative difference between the estimates of the components from one round to the following fell below 1%.
The asymptotic standard errors of the estimates were calculated. This required the computation of the information matrix, which is very extensive. Therefore, in this part of the analysis, the fixed effects of the year-season of calving, of the age at 1st calving and of the length of the preceding lactation were ignored, so that the Cholesky transformation was fully efficient. The relationships among bulls were also ignored in order to reduce the computational requirements. The information matrix I e of the transformed data was first calculated and, after back transformation, the information matrix I of the original data was obtained as it can be showed that: where D is the (6 x 6) matrix whose element (i,j) is (where 0 0c ) is the vector of (transformed) sire and residual variance and covariance components). Computations of the 1, matrix was made using Meyer's algorithm (1983a) and taking advantage of the simplifications the Cholesky transformation made possible.
Because of the computational costs, the two data sets were analyzed separately and the mean of the estimates calculated although the two data sets were not totally independent. Similarly, the asymptotic standard errors of the means of the estimates were obtained as if the asymptotic standard errors of the estimates in the 2 data sets were independent.

RESULTS
The estimates from the 2 data sets differed by less than 1 standard error, except for the variance of the protein content in the third lactation which differed by a little less than 2 standard errors. The asymptotic standard errors in the two data sets differed by less than 0.01 (Table III).
The estimates of the genetic parameters for the yields were similar. The phenotypic variances increased with lactation number. The change of the phenotypic standard deviations was proportional to the increase in the corresponding means and may be considered to be, at least partly, due to a scale effect. By contrast, the genetic components remained nearly constant from the first to the second lactation (except for protein yield where it differed by 26%, but this difference was not significant). Except for fat yield, the genetic component of the third lactation was the highest but the difference was not significant.
The heritabilities for the 3 lactations were therefore slightly but not significantly different (Table IV): the heritabilities for the first and third lactations were similar and higher than for the second lactation. The only exception was protein yield, where the heritabilities for the first and second lactations were equal to 0.18 and slightly smaller than for the third lactation (0.22). All genetic correlations were higher than 0.89. The correlation between the first and third lactations was smaller than the others, which except for useful yield were very similar. The same trend was observed for the phenotypic correlations: the correlations between adjacent lactations (first and second lactations or second and third lactations) was higher than between first and third lactations. All phenotypic correlations were between 0.53 and 0.65. Genetic parameters for the contents measured showed different trends. First the means for the different lactations were very similar (Table II) so there was no scale effect. The genetic components decreased with lactation number, whereas the phenotypic components remained constant. Thus the heritabilities decreased with lactation number, but the differences were not significant. The genetic correlations showed the same trend as for yields and were between 0.90 and 0.96. By contrast, the phenotypic correlations were higher than for yields (between 0.67 and 0.71) and did not vary much.

DISCUSSION
Choice of the method Different iterative algorithms may be used for REML estimation. They differ in convergence speed and computational requirements per round of iteration. Primarily, 3 algorithms have been advocated for analysis on selected data: Fisher's method (Meyer (1983a)), algorithms related to the E.M. algorithm of Dempster et al. (1977) and Meyer's algorithm (Meyer, 1986). Although Fisher's method has the highest convergence speed, it appears to be the most expensive (Meyer, 1986) and was therefore disregarded. Algorithms called &dquo;related to the E.M. algorithm&dquo; by Henderson (1985) converge very slowly but have the property of forcing the estimate within the parameter space. Meyer's &dquo;short cut&dquo; algorithm estimates the residual components via an algorithm related to the E.M. and the genetic components via Fisher's method. Thus the convergence speed is quicker than for algorithms related to the E.M., but the estimates may lie outside the parameter space. The first data set was analysed using Meyer's algorithm, and the estimate of the genetic correlation between first and third lactations was equal to 1.05 and between second and third to 1.09. Such estimates cannot be considered as maximum likelihood estimators (Harville, 1977). They cannot be used in the mixed model equations without using a transformation of the results such as the &dquo;bending&dquo; of Hayes and Hill (1981). In contrast, the E.M. type algorithm gave estimates that were within the parameter space. It required on our data set slightly less computations than Meyer's algorithm, although 8 iterations were needed before convergence was achieved (instead of 6 with Meyer's algorithm), as each iteration took 25% more time with Meyer's algorithm than with the E.M. type algorithm.
The Cholesky transformation makes the absorption of the herd effects quicker. In reference to the time necessary for the absorption of untransformed data, the same process took after Cholesky transformation 43% more time on the first round and 48% less time on the following. This transformation also spares a lot of computations for the estimation of the asymptotic standard errors.

Results of different methods
Maximum likelihood estimators can only take into account selection if it is based on observed data or on observed data and other independent variables. In this analysis, selection occurred both between generations (for the choice of the parents) and within a generation (by the end of each lactation). As the performance of the parents of the animals could not be analysed because of the computational requirements, only the later selection was considered. This is the case for most REML estimates of genetic parameters for lactations. The results are in accordance with studies using maximum likelihood related estimators (Tables V and VI). Except for , the authors used restricted maximum likelihood estimates but the algorithms varied: Fisher's method for Meyer (1983a) and Hagger et al. (1982); &dquo;short cut&dquo; for Meyer (1985a); E.M. type for Colaco et al. (1987), , Simianer (1986b), Swalve and Van Vleck (1987) and Tong et al. (1979). All considered a sire model except for Swalve and Van Vleck (1987), who used an animal model which may better take into account selection, since all relationships are included in the model (Sorensen and Kennedy, 1984) and thus did not observe any decrease in heritability for the second lactation. This decrease that most authors observe might be due to a selection bias. However Swalve and Van Vleck (1987) neglected relationships among herds and thus ignored selection across herds.
Henderson's methods lead to different results as they are affected by selection of the data Meyer and Thompson, 1984). Because of selection at the end of first lactation, they underestimate heritabilities for later lactations. For milk yield, the weighted means of the estimates in the literature are 0.26 for first lactation, 0.20 for second lactation and 0.17 for third lactation (Maijala and Hannah, 1974). In the first data set, Henderson's method III estimates for useful yield were respectively 0.21, 0.08 and 0.19. The first lactation estimate is in accordance with REML estimates because selection has not yet occurred. The third lactation estimate does not differ much either. As the criteria of selection depend less on milk production at the end of the second lactation than of first lactation, the selection bias may be less important. This result may also be, at least partly, due to sampling errors. Similarly, the decrease in the heritabilities for the later lactations for content measures may partly be due to the fact that the selection bias is not well removed for these variables, because the selection criteria are usually based more on milk yield than on content.
The parameters of the first 3 lactations are very similar. The heritabilities for the first and third lactations may, at least for yields, be treated as equal. The slight decrease of the heritability for the second lactation is not significant. It may, at least partly, be due to a selection bias which cannot be totally removed when using a sire model. The determinism of the second lactation may also be slightly different: this performance depends both on the dairy value of the animal and on its ability to recover from both growth and first lactation. The genetic correlations are very high but the correlation between first and third lactations is significantly different from 1. The older the cow is, the more disease it has had to resist and the more its ability to resist is important in the determinism of its lactations. However, the differences in the parameters of the lactations are very small.
It does not seem to be necessary to modify the current French genetic evaluation procedure which fits a repeatability model to the different lactations. All available lactations are taken into account because of the small mean herd size (34.5 cows per herd). Accuracy is increased using all records instead of first records only. This gain is due to both extra genetic information and increased degree of connectedness among herds (Meyer, 1983b). U f ford et al. (1979) reported such an increase even for young bulls whose daughters had only first lactations. Fitting a multi-trait model would imply a very large increase in computational requirements, as time needed for an iterative inversion of the coefficient matrix of the mixed model equations is proportional to the square of its size. But only a very small gain in accuracy could be expected. Simianer (1986a) and Schulte-Coerne (1983) estimate this increase to be less than 1% when 3 lactations are considered and all the genetic correlations are 0.80. The difference is expected to be even smaller in our case because the correlations are higher. However, they restricted their analysis to complete data (i.e. all animals were supposed to have made 3 lactations). In reality, some selection occurs. The selection bias can be totally removed only when the true genetic parameters are used, i.e. with the multi-trait model. But the difference between the 2 models is still expected to be small.