Optimum truncation points for independent culling level selection on a multivariate normal distribution, with an application to dairy cattle selection

is often practiced in breeding programs because extreme animals for some particular traits are rejected by breeders or because records on which genetic evaluation is based are collected sequentially. Optimizing these selection procedures for a given overall breeding objective is equivalent to finding the combination of truncation thresholds or culling levels which maximizes the expected value of the overall genetic value for selected animals. A general Newton-type algorithm has been derived to perform this maximization for any number of normally distributed traits and when the overall probability of being selected is fixed. Using a powerful method for the computation of multivariate normal probability integrals, it has been possible to undertake the numerical calculation of the optimal truncation points when up to 6 correlated traits or stages of selection are considered simultaneously. The extension of this algorithm to the more complex situation of maximizing annual genetic response subject to nonlinear constraints is demonstrated using a dairy cattle model involving milk production and a secondary trait such as type. Consideration is given to three of the four pathways of selection: dams of bulls; sires of bulls; and sires of cows.

(received 1 March 1988, accepted 19 September 1988) Summary &mdash; Independent culling level selection is often practiced in breeding programs because extreme animals for some particular traits are rejected by breeders or because records on which genetic evaluation is based are collected sequentially. Optimizing these selection procedures for a given overall breeding objective is equivalent to finding the combination of truncation thresholds or culling levels which maximizes the expected value of the overall genetic value for selected animals. A general Newton-type algorithm has been derived to perform this maximization for any number of normally distributed traits and when the overall probability of being selected is fixed. Using a powerful method for the computation of multivariate normal probability integrals, it has been possible to undertake the numerical calculation of the optimal truncation points when up to 6 correlated traits or stages of selection are considered simultaneously. The extension of this algorithm to the more complex situation of maximizing annual genetic response subject to nonlinear constraints is demonstrated using a dairy cattle model involving milk production and a secondary trait such as type. Consideration is given to three of the four pathways of selection: dams of bulls; sires of bulls; and sires of cows.
However, this approach is occasionally not used for two main reasons: a) Practitioners sometimes emphasize the need to cull extreme animals (deviants for some traits) because they are deemed undesirable. Then, the selection objective is implicitly recognized as being nonlinear. The cost of such a practice with respect to a strict application of the optimal linear index may not be justifiable when this nonlinearity is unimportant or questionable.
b) The records required to compute the selection index are not always available simultaneously and/or their cost does not justify their collection for all the candidates for selection. Therefore selection schemes involve different stages which correspond to truncations on the joint distribution of all possible records. Within these constraints, it is potentially interesting to evaluate and improve the efficiency of practical selection procedures by determining optimal conditions of application of independent culling level selection. As far as we know, the published works on this topic have been greatly limited by the number of variables considered. Generally speaking, studies on the algebraic derivation of optimum truncation thresholds with corresponding numerical computation have dealt with no more than 2 variables (Namkoong, 1970;Evans, 1980;Cotterill and James, 1981;Smith and Quaas, 1982). Recently, Ducrocq and Colleau (1986) treated examples with 3 variables to illustrate potential uses of a numerical method to compute multivariate normal probabilities. In practice, though, the number of traits or selection stages involved may be significantly larger. Moreover, the algorithms previously considered have been specifically developed for a given number n of variables and their extension to any n is not obvious though algebraic conditions which must be verified at the optimum have been reported (Jain and Amble, 1962;Smith and Quaas, 1982).
Indeed, as a consequence of the limitation of the number of variables considered, the optimization problem has been restricted to rather simple types of objective functions, which may not adequately summarize the overall efficiency of selection schemes. In particular, authors have not considered functions such as the annual genetic gain of the selection objective -computed using Rendel and Robertson's (1950) formula -in relatively complex situations (e.g. involving the 4 paths of transmission of genetic progress). This paper presents basic yet general, i.e. for any n algorithms which can be used for the computation of optimal truncation points for a broad class of objective functions with one or more constraints. Theory is developed and applications are presented for the general multivariate problem considered by Smith and Quaas (1982). A practical example of application in the dairy cattle context is described. Corresponding numerical results are given. Solution of Smith and Quaas' problem in the general case Statement of the problem Let u, x l , ....x,, be n+1 random variables with joint multivariate normal distribution: u is the breeding objective and the x i s are the observed variables.
The problem is to find c l , C2 &dquo;&dquo; Cn so as to maximize Ep (u) subject to: where P is given and represents the overall fraction of candidates selected. Notations Let !&dquo; (x; R n ) be the standard multivariate normal density of dimension n with variance covariance matrix R n .

Let
We need the following recursive definitions for distributions conditioned on q::;;n-1 variates.

Solution
Using the general result of Jain and Amble (1962) we have: Since Q(c 1 ,... c n ) is to be equal to a constant P, the maximization of Ep (u) is tantamount to the maximization of N(c l ,..., c n ). The constraint Q(c,,... c n ) = P is incorporated using the method of Lagrange multipliers (Bass, 1961, p. 928; Smith and Quaas, 1982): the optimal truncation points are those for which the partial derivatives of the function f (w) = N(c 1' &dquo; e n ) + X(Q(c,... c n ) -P) with respect to w'= ( C1' &dquo; '' c n , X) are 0. Â. is called a Lagrange multiplier.
The resulting system of nonlinear equations in w' = ( C1' ... c n . X) is solved iteratively using the multidimensional Newton's method (Dennis and Schnabel, 1983). Denote as wct> the approximate solution at iteration t (w(o) is a given starting value). A better estimate w(t+ 1 )is computed from: The final solution wIt) = w * is obtained when is sufficiently small, where I I h I I denotes any norm.
As long as the starting value w(o) is not too far from w * (generally, w(o) = 0 seems to be a robust initial value), convergence is very fast (quadratic convergence: Dennis and Schnabel, 1983). c * is a local maximum for E (u), provided ' is positive definite, but nothing guarantees that w * is a global maximum for f (w). Now, note that: So we can write: where: Another method exists for the derivation of these expressions, first reasoning on derivatives of integrals and then using Jain and Amble's (1962) formula on conditional distributions. This leads to more compact expressions but may be less flexible for considera-tions on several u (general problem with several constraints on several objectives) (see Appendix).
Expressions (3) to (11) include all the elements required for the computation of the vector 8f(w)/Bw and the matrix (8 2 f(w)/Bw Sw * ) in (1). In particular, It can be observed that the equations (6f(w)/6c i ) = 0 in (12) for 1 !i!n are linear in X. The size of the system of equations to be solved can be easily reduced by absorption of the Lagrange multiplier. For example, we have: is equivalent to: Derivatives with respect to the c i s of the equations in (17) are required for the application of Newton's method as in (2). They are readily derived using (6) to (8).
Numerical applications Studies on independent culling level selection have been mainly limited to 2-trait selection probably because general and efficient programs to compute the multivariate normal probability integrals in Jain and Amble's formula were not available for dimensions > 2. However, easily programmable algorithms exist. In particular, Dutt (1973Dutt ( , 1975 and Dutt and Soms (1976) proposed a general method characterized by good precision when correlation coefficients and truncation points are not too extreme. For more details on this method, its precision and computation times, see Ducrocq and Colleau (1986). Dutt's . technique is well suited for numerical applications of the optimization algorithm presented in this paper when up to 6 selection stages or traits are considered.
Algebraic and numerical results are equivalent to those given by Smith and Quaas (1982).
In (9), Q;ik = 1. The algorithm described here leads to the same results as those presented in Ducrocq and Colleau (1986).
3) . n = 4 to 6. Consider for example, n traits with r; i = (-1 (j/20i) for 1 <_i<_j<_n and with economic weight m i = 1+i/20 1!i!n. Table I presents the truncation points qs on these traits which maximize Ep (u I c1,... c&dquo;) when the overall selected fraction is P = 0.25, 0.025, 0.001. At iteration 0, the c i s were taken equal to 0. The stopping criterion for the Newton's iteration was: where E i was the i th left hand side of system (17).
Convergence was fast and depended on how far the initial value of the truncation points was from the solution. Note, however, that in the examples presented in Table I, correlations between variables are not very high Q r;! ! 1 S 0.3) and the weights of the different traits are of the same magnitude. When this is not the case, the optimal selection procedure may involve no selection at all on one or several of these traits. The same observation applies to small overall selection intensity (Young, 1961;Namkoong, 1970;Smith and Quaas, 1982;Tibau I Font and Ollivier, 1984;Ducrocq and Colleau, 1986). In limiting cases (with very low or very high selection intensity on one or several traits or when correlations are extreme), it should be remembered that the precision of Dutt's algorithm for computation of multivariate probability integrals may be unsatisfactory (Ducrocq and Colleau, 1986). Then alternative methods may have to be used (e.g., Russell et al., 1985). An application in the dairy cattle context Assumptions Dairy cattle selection is performed through a sequence of stages which characterizes the transition from one generation (g) to the next (g+1 ). These first include selection criteria corresponding to the transition between generations g and g+1 (reproductive stage), followed in the course of time by those criteria used during generation g+1, before the next reproductive cycle. Our approach for the optimization of these successive selection stages relies on the assumption of multivariate normality for these 2 criteria when candidates for selection are born. Such an assumption is plausible in the additive polygenic context, especially when heritability values are low (Bulmer, 1980, p. 154; see also Smith and Hammond, 1987, for a discussion on this point). A more strident assumption is that the dispersion parameters of the joint multivariate distribution remain constant through the different selection cycles.
Breeding objective and selection stages Assume that the selection objective in a dairy cattle breed is a linear combination of 2 traits: &dquo;milk production&dquo; and a secondary trait such as &dquo;type&dquo; (both of these traits may be themselves linear combinations of more specific characters). A possible sequence of selection stages which approximates what is often done in practice is the following (Figure 1): 1) Dams of bulls (DB) of generation g are selected based on their estimated breeding values X1 (for milk) and X2 (for type), with respective thresholds C1 and c 2 on the standardized variables. These dams of bulls are mated to sires of bulls (SB) of generation g.
2) The sons of these cows are progeny tested. Sires of cows (SC) and sires of bulls of generation g+1 are then selected according to their estimated breeding values x 3 (for milk) and X4 (for type). Truncation thresholds on these 2 variables are different for SC and SB (c 3 , c 4 and c 5 c 6 , respectively).
Selection of DB can be modelled as if it were performed at birth of the male calves. This is essential in order to be able to invoke the restoring of multivariate normality at each generation. Let RP be the registered (with known pedigree) and recorded population of cows and let y denote the proportion of these cows which can be potential dams of bulls (e.g. !y= 0.53 if Al sons are selected from cows with at least 2 known lactations). If it is assumed (as in Ducrocq, 1984) that an average of n d = 6 potential dams must be selected in order to obtain one male calf entering progeny test, it can be considered that DB selection is performed by truncation on the estimated breeding values x, and X2 of the dams of n b = (y RP)/n d male calves.
The expression of the annual genetic gain given by Rendel and Robertson (1950) is: Selection on the dam of cow path is ignored (I! = 0).
where do is the fraction of the whole population bred to young sires (do = Ty RP/T), i.e.:

Constraints
Three constraints are added here: 1) The fraction Ty of the population bred to young sires is considered as constant, since in practice this is usually the limiting factor for the extension of progeny test. In this example, the number of recorded daughters per young sire n v y is also assumed constant: then, the number ny of young sires progeny-tested each year is fixed, as well as their repeatablity.
2) The number of sires of cows selected each year is determined by the number of cows (= (T-RP) + (1-Ty) RP) to be bred to proven sires in the whole population (T) and the total number of doses produced by a given sire during his lifetime (AI). --,--, 3) The number of sires of bulls retained each year is constrained to be equal to n sB , the number below which problems of inbreeding and reduction of genetic variability are feared.

Numerical methods and results
When constraints (22), (23) and (24) are satisfied, equations (18) to (21) lead to the following result: where L, the sum of the generation intervals over the 4 paths, is a constant in our case.
The first and second derivatives of f (w) are readily obtained using the general formulae given in the preceding sections. The 3 Lagrange multipliers are eliminated through trivial absorption. The nonlinear system to solve then involves 6 unknowns: the 6 truncation thresholds. Solutions obtained using Newton's method are presented in Table III, where parameters take the values given in Table II. The stopping criterion for the Newton's iterations was: . where E i * was the ith left hand side of the absorbed system of equations.
Convergence was fastalways less than 8 iterationsas long as the starting value c ( 0 ) of c = {c ; a ' , 6 was not too far from the solution. In contrast with the numerical results presented in Table I, c !o! = 0 does not lead to convergence, because the values of c in the next iterates are found to be extremely large or low; i.e., totally unrealistic and in regions where the precision of Dutt's method is not good. This behavior is, at least indirectly, a consequence of the very large selection intensity of proven sires. To avoid this situation, more realistic starting values must be chosen. For example, c ( 0 ) can be computed as the truncation points for which the conditions (22), (23) and (24) are satisfied with a given distribution of selection efforts between milk and type, i.e., c ( 0 ) is such that Q (c!O) = d,/a, Q (c,co!! c!<0» = d!! Q (c,tot, c!<o>, C3 (o)) = d 2 la, etc... for some a. In most cases presented in Table III, a = 0.75. For strongly unbalanced weights for the 2 traits (5:1), convergence is obtained only for even larger values of a (e.g., a = 0.99).
It should be noted that selection on type was considered posterior to selection on milk production. This assumption has been chosen because this corresponds to what is done in practice for dams of bulls in France: only the best cows for milk production are evaluated for type. This does not influence the values of optimal thresholds but only their use when results are presented in terms of fraction selected at each stage.
The results in Table III underline the sensitivity of this type of computation to the economic weights assigned to the 2 traits and their genetic correlation. According to the relative economic weight for type, the optimal situation may correspond to virtually no selection on type (weights 5:1) or to a culling on type evaluation of between a quarter and a half of the candidates (weights 2:1) on each path. When the genetic correlation between milk and type varies from slightly positive (0.15) to slighly negative (-0.15) (thus covering the range of values most often indicated in the literature) optimum selection intensities and overall genetic gain for the objective function are not dramatically modified. But then, the genetic trend for type varies from a very favorable increase (for a positive correlation) even when no selection on type is performed, to a not negligible decline (for a negative correlation and weights 3:1 or 5:1 ). These features were also pointed out by Ducrocq (1984).

Conclusion
The previous example clearly demonstates that, at least in certain situations, it is possible to algebraically and numerically develop algorithms to compute optimal selection intensities in relatively complex multistage breeding schemes and with several nonlinear constraints. They are not limited by drastic conditions such as 2-stage selection, uncorrelated traits and/or very simple optimization criteria. Indeed, more complex situations than that presented in the example can be envisaged. For example, the overlap of generations of bulls and cows can be accounted for in the determination of genetic superiorities and generation intervals, extending to the multivariate case the techniques described in Hill (1974), Mocquot (1974, 1976) and Ducrocq and Quaas (1988): dividing each population of candidates into homogeneous cohorts of animals of same age, sex and reproductive role, intracohort selection differentials can be computed under the assumption that selection is performed by retaining all the individuals whose estimated genetic values are above a unique set of truncation thresholds for the n traits or stages considered (Ducrocq, 1984). Then, computation of elements in (18) and their derivatives is tedious but perfectly feasible. However, before undertaking optimization studies of complex schemes, it should be remembered that troublesome problems of convergence may occur such as those found in the application described. In the particular case of the dairy cattle context, convergence problems were also found when we tried to relax the assumption of fixed number of young sires sampled (ny) or the number of daughters per young sire (n!), i.e., their repeatability * . Obviously, further research is needed in a widening field.
Finally, the determination of an optimal selection policy, which is &dquo;optimal&dquo; only from a purely genetic standpoint, should be followed by a sensitivity analysis. It appears that the annual genetic gain AG, varies only slightly over a wide range of values for some parameters. Situations close to the optimum may be of greater interest because they are simple, more practical or less expensive to implement. The knowledge of the true optimum is nevertheless required to assess this &dquo;suboptimality&dquo;.