Assessment of a Poisson animal model for embryo yield in a simulated multiple ovulation-embryo transfer scheme

Estimation and prediction techniques for Poisson and linear animal models were compared in a simulation study where observations were modelled as embryo yields having a Poisson residual distribution. In a one-way model (fixed mean plus random animal effect) with genetic variance (0';) equal to 0.056 or 0.125 on a log linear scale, Poisson marginal maximum likelihood (MML) gave estimates of 0 '; with smaller empirical bias and mean squared error (MSE) than restricted maximum likelihood (REML) analyses of raw and log-transformed data. Likewise, estimates of residual variance (the average Poisson parameter) were poorer when the estimation was by REML. These results were anticipated as there is no appropriate variance decomposition independent of location parameters in C he linear model. Predictions of random effects obtained from the mode of the joint posterior distribution of fixed and random effects under the Poisson mixed-model tended to have smaller empirical bias and MSE than best linear unbiased prediction (BLUP). Although the latter method does not take into account nonlinearity and does not make use of the assumption that the residual distribution was Poisson, predictions were essentially unbiased. After log transformation of the records, however, BLUP led to unsatisfactory predictions. When embryo yields of zero were ignored in the analysis, Poisson animal models accounting for truncation outperformed REML and BLUP. A mixed-model simulation involving one fixed factor (15 levels) and 2 random factors for 4 sets of variance components was also carried out; in this study, REML was not included in view of highly heterogeneous nature of variances generated on the observed scale. Poisson MML estimates of variance components were seemingly unbiased, suggesting that statistical information in the sample about the variances was adequate. Best linear unbiased estimation (BLUE) of fixed effects had greater empirical bias and MSE than the Poisson estimates from the joint posterior distribution, with differences between the 2 analyses increasing with the genetic variance and with the true values of the fixed effects. Although differences in prediction of random effects between BLUP and Poisson joint modes were small, they were often significant and in favor of those obtained with the Poisson mixed model. In conclusion, if the residual distribution is Poisson, and if the relationship between the Poisson parameter and the fixed and random effects is log linear, REML and BLUE may lead to poor inferences, whereas the BLUP of breeding values is remarkably robust to the departure from linearity in …

2 analyses increasing with the genetic variance and with the true values of the fixed effects. Although differences in prediction of random effects between BLUP and Poisson joint modes were small, they were often significant and in favor of those obtained with the Poisson mixed model. In conclusion, if the residual distribution is Poisson, and if the relationship between the Poisson parameter and the fixed and random effects is log linear, REML and BLUE may lead to poor inferences, whereas the BLUP of breeding values is remarkably robust to the departure from linearity in terms of average bias and MSE.
Optimizing embryo yields could be important for other reasons as well. For instance, with greater yields, the gap in genetic gain between closed and open nucleus breeding schemes could be narrowed (Meuwissen, 1991). Furthermore, because of possible antagonisms between production and reproduction, it may be necessary to use some selection intensity to maintain reproductive performance (Freeman, 1986). Also, if yield promotants, such as bovine somatotropin, are adopted, the relative economic importance of production and reproduction, with respect to genetic selection, will probably shift towards reproduction. Finally, if cytoplasmic or nonadditive genetic effects turn out to be important, it would be desirable to increase embryo yields by selection, so as to produce the appropriate family structures (Van Raden et al, 1992) needed to fully exploit these effects. Lohuis et al (1990) found a zero heritability of embryo yield in dairy cattle. Using restricted maximum likelihood (REML), Hahn (1992) estimated heritabilities of 6 and 4% for number of ova/embryos recovered and number of transferable embryos recovered, respectively, in Holsteins; corresponding repeatabilities were 23 and 15%. Natural twinning ability may be closely related to superovulatory response in dairy cattle, as cow families with high twinning rates tend to have a high ovarian sensitivity to gonadotropins, such as PMSG and FSH (Morris and Day, 1986). Heritabilities of twinning rate in Israeli Holsteins were found to be 2%, using REML, and 10% employing a threshold model (Ron et al, 1990).
Best linear unbiased prediction (BLUP) of breeding values, best linear unbiased estimation (BLUE) of fixed effects, and REML estimation of genetic parameters are widely used in animal breeding research. However, these methods are most appropriate when the data are normally distributed. The distribution of embryo yields is not normal, and it is unlikely that it can be rendered normal by a transformation, particularly when mean yields are low and embryo recovery failure rates are high. Analysis of discrete data with linear models, such as those employed in BLUE or REML, often results in spurious interactions which biologically do not exist ((auaas et al, 1988), which, in turn, leads to non-parsimonious models.
It seems sensible to consider nonlinear forms of analysis for embryo yield. These may be computationally more intensive than BLUP and REML, but can offer more flexibility. The study of Ron et al (1990) suggests that nonlinear models for twinning ability may have the potential of capturing genetic variance for reproduction that would not be usable by selection otherwise. For example, threshold models have been suggested for genetic analysis of categorial traits, such as calving ease (Gianola and Foulley, 1983;Harville and Mee, 1984). In these models, gene substitutions are viewed as occurring in a underlying normal scale. However, the relationship between the outward variate (which is scored categorically, eg, 'easy' versus 'difficult' calving) and the underlying variable is nonlinear and mediated by fixed thresholds. Selection for categorical traits using predictions of breeding values obtained with nonlinear threshold models was shown by simulation to give up to 12% greater genetic gain in a single cycle of selection than that obtained with linear predictors (Meijering and Gianola, 1985). Because genetic gain is cumulative, this increase may be substantial. The use of better models could also improve (eg, smaller mean squared error (MSE)) estimates of differences in embryo yield between treatments.
In the context of embryo yield, an alternative to the threshold model is an analysis based on the Poisson distribution. This is considered to be more suitable for the analysis of variates where the outcome is a count that may take values between zero and infinity. A Poisson mixed-effects model has been developed by Foulley et al (1987). From Foulley and Im (1993) and Perez-Enciso et al (1993). In the Poisson model, the link is the logarithmic function. such that Above, 0' = ![3', u'], 13 is a p x 1 vector of fixed effects, u is a q x 1 vector of breeding values, and w' = [x! z!] is an incidence row vector relating 0 to r¡ i . Let X = fx'l and Z = fz'l be incidence matrices of order n x p and n x q, respectively, such that: In a Bayesian context, Foulley et al (1987) employ the prior densities and where A is the matrix of additive relationships between animals and Q u is the additive genetic variance. Given o,', Foulley et al (1987) calculate the mode of the joint posterior distribution of (3 and u with the algorithm where t denotes iterate number, and where y is the vector of observations. Note that the last term in [8] can be regarded as a vector of standardized (with respect to the conditional Poisson variance) residuals.
Marginal maximum likelihood (MML), a generalization of REML, has been suggested for estimating variance components in nonlinear models H6schele et al, 1987). An expectation-maximization (EM) type iterative algorithm is involved whereby where T = trace(A-1 C u &dquo;), such that and u is the u-component solution to [7] upon convergence for a given o,' value. In !9!, k pertains to the iteration number, and iterations continue until the difference between successive iterates of [9], separated by nested iterates of [7], becomes arbitrarily small. The above implementation of MML is not exact, and arises from the approximation (Foulley et al, 1990) SIMULATION EXPERIMENTS A one-way random effects model Embryo yields in two MOET closed nucleus herd breeding schemes were simulated.
Breeding values (u) for embryo yields for n, and n d base population sires and dams, respectively, were drawn from the distribution u N N(0, I!u), where 0' ; had the values specified later. The dams were superovulated, and the number of embryos collected from each dam were independent drawings from Poisson distributions with parameters: where 1 is a n d x 1 vector of ones, p is a location parameter and u d is the vector of breeding values of the n d dams. In nucleus 1, f ..l = ln(2), whereas in nucleus 2 p = ln(8). Note from [12] that for a given donor dam d i , so, in view of the assumptions, which implies that the location parameter can be interpreted as the mean of the natural log of the Poisson parameters in the population of donor dams. It should be noted, as in Foulley and Im (1993) that

Thus
The sex of the embryos collected from the donor dams was assigned at random (50% probability of obtaining a female), and the probability of survival of a female embryo to age at first breeding was !r = 0.70 in nucleus 1, and 7 r = 0.60 in nucleus 2. This is because research has suggested that embryo quality and yield from a single flush tend to be negatively associated (Hahn, 1992). Thus the expected number of female embryos surviving to age at first breeding produced by a given donor dam d i is, for i = 1, 2, ... n d , and, on average, The genetic merit for embryo yield for the ith female offspring, uo! , was generated by randomly selecting and mating sires and dams from the base population, and using the relationship: where USi and u di are the breeding values of the sire and dam, respectively, of offspring i, and the third term is a Mendelian segregation residual; z -N(0,1).
As with the dams, the vector of true Poisson parameters for female offspring was 71 0 = exp[1p + u o] where u o represents the vector of daughters' genetic values. The unit vector 1 in this case would have dimension equal to the number of surviving female offspring. Embryo yields for daughters were sampled from a Poisson distribution with parameter equal to the ith element of X o .
In each of the 30 replicates of each population, variance components for embryo yield were estimated employing the following methods: 1) Poisson MML as in Foulley et al (1987); 2) REML as if the data were normal; 3) truncated Poisson MML excluding counts of zero, and using the formulae of Foulley et al, (1987); 4) REML-0, ie REML applied to the data excluding counts of zero; and (5) REML-LOG, which was REML applied to the data following a log transformation of the non-null responses while discarding the null responses. Empirical bias and mean squared error (MSE) of the estimates, calculated from the 30 replicates, were used for assessing performance of the variance component estimation procedure. Because the probability of observing a zero count in a Poisson distribution with a mean of 8 is very low, the truncated Poisson and REML-0 analyses were not carried out in nucleus 2. Likewise, breeding values were predicted using the following methods: 1) the Poisson model as in [7] with the true o!, and taking as predictors , where u is the vector of breeding values of sires, dams, and daughters; 2) BLUP (1!* + u * ) in a linear model analysis where the variance components were the average of the 30 REML estimates obtained in the replications and the asterisk denotes direct estimation of location parameters on the observed scale; 3) a truncated Poisson analysis with the true a and predictors as in 1); 4) BLUP-0, as in 2) but excluding zero counts, and using the average of the 30 REML-0 estimates as true variances; and 5) BLUP-LOG, as in 2) after excluding zero counts and transforming the remaining records into logs. The average of the 30 REML-LOG estimates of variance components was used in this case. BLUP-LOG predictors of breeding values were expressed as exp [lti + u] where p and u are solutions to the corresponding mixed linear model equations. Hence, all 5 types of predictions were comparable because breeding values are expressed on the observed scale. As given in !12!, the vector of true Poisson parameters or breeding values for all individuals was deemed to be A = exp[1p + u]. Average bias and MSE of prediction of breeding values of dams and daughters were computed within each data set and these statistics were averaged again over 30 further replicates. Rank correlations between different estimates of breeding values were not considered as they are often very large in spite of the fact that one model may fit the data substantially better than the other (Perez-Enciso et al, 1993).
A mixed model with two random effects The base population consisted of 64 unrelated sires and 512 unrelated dams, and the genetic model was as before. The probability of a daughter surviving to age at first breeding was 7 r = 0.70. Embryo yields on dams and daughters were generated by drawing random numbers from Poisson distributions with parameters: where p is a fixed effect common to all observations, H = {H i } is a 15 x 1 vector of fixed effects, s = { Sj } '&dquo; N(0,Iu£) is a 100 x 1 vector of unrelated 'service sire' effects, u = j Uk } -N(0, A U2 ) is a vector of breeding values independent of service sire effects, and 0 '; and 0 '; are appropriate variance components.
The values of + Hi were assigned such that: Thus, in the absence of random effects, the expected embryo yield ranged from 1 to 15. Each of the 15 values of fl + H i had an equal chance of being assigned to any particular record.
Service sire has been deemed to be an important source of variation for embryo yield in superovulated dairy cows (Lohuis et al, 1990;Hasler, 1992). However, no sizable genetic variance has been detected when embryo yield is viewed as a trait of the donor cow (Lohuis et al, 1990;Hahn, 1992). This  For each of the 4 sets of variance parameters, 30 replicates were generated to assess the sampling performance of Poisson MML in terms of empirical bias and square root MSE. Relative bias was empirical bias as a percentage of the true variance component. Coefficients of variation for REML and MML estimates of variance components were used to provide a direct comparison as they are expressed on different scales. REML estimates were also required in order to compare estimates of fixed effects and predictions of random effects obtained under a linear mpdel analysis with those found under the Poisson model. MML and REML estimates were computed by Laplacian integration (Tempelman and Gianola, 1993) using a Fortran program that incorporated a sparse matrix solver, SMPAK (Eisenstat et al, 1982) and ITPACK subroutines (Kincaid et al, 1982) to set up the system of equations !7!. For REML, this corresponds to the derivativefree algorithm described by Graser et al (1987) with a computing strategy similar to that in Boldman and Van Vleck (1991).
As in the one-way model, averages of REML estimates of the variance components obtained in 30 replicates were used in lieu of the 'true' values (which are not well defined) to compute estimates of fixed effects and predictions of random effects in the linear model analysis; for the Poisson model, the true values of the variance components were used. Empirical biases and MSEs of the estimates of fixed effects obtained with the linear and with the Poisson models were assessed from another 30 replicates within each set of variance components. One more replicate was then generated for each variance component set, from which the empirical average bias and MSE of prediction of service sire and animal random effects were evaluated.
In order to make comparisons on the same scale, the Poisson model predictands of the random effects were defined to be b.exp(s) for service sires and b!exp(u) for additive genetic effects, respectively; b is the 'baseline' parameter: In In the linear mixed model, the predictors were defined to be: and for service sire and genetic effects, respectively. Here the unit vectors 1 are of the same dimension as the respective vectors of random effects and the asterisk is used to denote direct estimation of location parameters on the observed scale.
Estimators for fixed effects were also expressed on the observed scale. The true values of the fixed effects were deemed to be i = exp(p + H i ) for i = 1, 2, ... 15 as in !16!. Estimators for fixed effects under the Poisson model were therefore taken to be exp(ti+!) for i = 1, 2, ... 15. As the linear mixed model estimates parameters on an observable scale, estimators for fixed effects were taken to be R * + H!&dquo;.

One-way model
Means and standard errors of estimates of the genetic variance ( 0 &dquo;) for the five procedures are given in table II and MSEs of the estimates are given in table III. Clearly, estimates obtained with REML and REML-0 were extremely biased; this is so because the genetic components obtained are not on the appropriate scale of measurement (ie the canonical log scale). The problem was somewhat corrected by a logarithmic transformation of the records. For E(A i ) * 2 and Q u = 0.056, the REML-LOG, Poisson and Poisson-truncated estimators were nearly unbiased (within the limits of Monte-Carlo variance), but the Monte-Carlo standard errors were much larger for REML-LOG. For Q 2 u = 0.125, the Poisson estimates were biased downwards (P < 0.05) for both values of E(A i ), while those of REML-LOG were biased upwards and significantly so with E(A; ) x5 8. In a one-way sire threshold model, H6schele et al (1987) also found downward biases for the MML procedure. In spite of these small biases, however, the MSEs of the Poisson estimates (table III) were much lower than those of REML-LOG. The very large (relative to Poisson and REML-LOG) MSEs of the REML and REML-0 procedures illustrate the pitfalls incurred in carrying out a linear model analysis when the situation dictates a nonlinear analysis, or a transformation of the data.
A linear one-way random effects model, however, can be contrived in which case it can be shown that REML may actually estimate somewhat meaningful variance components on the observed scale. Presuming that multiple records on an individual is possible, the variance of Yi! (with subscripts denoting the jth record on the ith individual) can be classically represented as: which from [2] can be written as: such that from [13c] and results presented by Foulley and Im (1993): The covariance between different records on the same individual (ie cov(Y!, Y!!!) can be used to represent the variance of the random effects.
Given independent Poisson sampling conditional on u i , the first term of the above equation is null, and Thus a one-way random linear model that has the same first and second moments as Y ij is where Y2! is the jth record observed on the ith animal, is the overall mean, u* is the random effect of the ith animal and e ! is the residual associated with the jth record on the ith animal. Here (i * = exp()i+o-!/2) ui has null mean and variance a 2 * = exp(2p)exp(u £ ) [exp ( g 2 ) 1] and eij has null mean and variance a e 2* = exp (p -f-Q!/2). The empirical mean REML estimates reported in However, REML estimates of residual variance, rather, of E(A i ), appeared to be biased upwards (P < 0.01) for the higher genetic variance and higher Poisson mean population (table IV). The MSEs of Poisson estimates of average residual variances were much smaller than those obtained with REML, especially in the populations with a higher mean. REML-0 was even worse than REML, both in terms of bias ( MSEs when E(Ai ) * 2, but Poisson had smaller MSEs than BLUP (P < 0.05) for both dams and daughters when E(Aj ) x5 8. These small differences between BLUP and Poisson are somewhat surprising in view of the different scales of prediction, and their practical significance is an open question. It was also surprising that the differences between BLUP and Poisson tended to be more significant with E(Ai ) * 8 than with E(Ai ) x 5 2, since it is known that the Poisson distribution approaches a normal distributions as A i increases (Haight, 1967). However, this may be due to a higher power of the test when detecting a larger difference. Another explanation may be that a lower E(A i ) leads to a lower 'pseudoheritability' and, hence, a lower degree of association between phenotypes and breeding values. In this case, the linear and Poisson models may differ less when predicting breeding values because of a higher degree of shrinkage towards zero. When counts of zero were excluded, the Poisson-truncated method had always smaller MSEs of prediction of breeding values than the BLUP-0 method.

Mixed model
Because variance components estimated by MML and REML are on different scales, empirical coefficients of variation (CV) were used to provide a basis of comparison (fig 1). Clearly, REML estimates were more variable than their Poisson counterparts. No clear pattern with respect to increasing values of the variance components emerged, except that CVs for both MML service sire and genetic estimates seemed relatively more stable while CVs for REML service sire estimates steadily increased for values of ufl larger than 0.0250 (ie 0' ; = 0.0500). For count data with low means, variance components that are estimated under linear mixed-effects models, more general than the one-way random effects model, are virtually meaningless. In Poisson-generated data, these components are highly heteroskedastic from one observation to the next, depending on both experimental design and location parameters (Foulley and Im, 1993).
Relative biases of the MML estimates are given in figure 2. Relative biases were less than 4% for all 4 sets of variance components, with no clear trend with respect to the true values of variance components. Using t-tests, these biases did not differ from zero. Unlike results obtained with threshold models (eg, H6schele et al, 1987;Simianer and Schaeffer, 1989), a small subclass (equal to 1 in the Poisson animal model) did not lead to detectable bias of variance component estimates in the Poisson mixed model.
Relative errors (square roots of the empirical MSEs, expressed as a percentage of the true variance component values) are given in figure 3. Trends with respect to the size of the true variances were somewhat opposite for service sire and genetic variance component estimates. Relative error for the genetic component decreased as the true variance increased, whereas the error of the service sire estimates increased somewhat with the true value of the parameter. The relative errors of MML estimates were almost identical to their empirical CVs (fig 1) because of the small bias, as shown in figure 2.
Empirical biases of fixed-effect estimates obtained with linear and Poisson procedures are given in figures 4-7 for the 4 different sets of true variance components. The 2 methods gave estimates that were biased upwards, but the bias was larger for BLUE in all 4 cases. This apparent paradox can be explained by the fact that BLUE is an unbiased estimator of: and not exp(p + H i ). Empirical biases for the 2 methods, and their difference, increased with increasing values of fixed effects and with higher values of variance components. The upward biases of the Poisson estimates were generally stable across sets of variance components, being always less than 0.5. However, the bias of the BLUEs of the fixed effects increased substantially as variance increased. Although Poisson estimates of fixed effects were manifestly biased upwards at higher embryo yields, the magnitude of their bias was several times smaller than that of BLUE estimates. Empirical MSEs of the fixed effects estimates are depicted in figures 8-11 for each of the 4 sets of variance components. As with empirical biases, MSEs of the linear model and Poisson estimates, and the differences between the MSEs of the 2 procedures tended to increase with increasing values of the variance components and with the true values of the fixed effects. The MSEs of the Poisson estimates were again more stable across sets of variance components and, although tending to increase with the value of the fixed effects, were much smaller than those of BLUE estimates. For example (fig 11), when embryo yield was around 15, the MSE of BLUE was about 6 times larger than that of Poisson estimates.
Because the values of fixed effects were constructed such that: Statistics associated with prediction of genetic effects are depicted in figures 16 and 17, with results presented separately for base generation sires and dams, and for their female progeny surviving to age at first breeding. Poisson-predicted breeding values tended to be biased downwards, whereas BLUP predictions were biased upwards; however, average bias was smaller for Poisson predictions. BLUP sire solutions were significantly biased (P < 0.01) at the 2 highest variance components. Dam and daughter BLUP solutions were significantly biased (P < 0.01 in all cases, except P < 0.05 for dams at ufl = 0.0125). Poisson dam solutions were biased (P < 0.01) only at U2 -0.0250, while Poisson daughter solutions were biased at o l 2= u 0.0250 (P < 0.01) and a = 0.0375 (P < 0.05). Certainly, all predictions reflect uncertainty in the baseline estimates of fixed effects given in [17] and [20] for Poisson and BLUP models, respectively, and upward biases in BLUP predictions of random effects reflect upward biases for the baseline estimates (fig 4-7). Differences in empirical average MSEs of prediction of breeding values between the 2 procedures are depicted in figure 17. Although the differences between the 2 methods were small, a paired t-test gave significant differences for daughters (P < 0.01 for U2 = 0.0125,0.0375, and 0.0500; P < 0.05 for U 2 -0.0250) and for base generation dams, except at or = 0.0125 (P < 0.01 for '7! = 0.0250,0,0500; P < 0.05 for 0 '; = 0.0375). Differences in empirical average MSE of predictions of base generation sires were not significant. CONCLUSIONS This study evaluated the sampling performance of estimators of location and dispersion parameters in Poisson mixed models, and of predictors of breeding values in the context of simulated MOET schemes. Records on embryo yield were drawn from Poisson distributions for 4 populations characterized by appropriate parameter values. Analyses were carried out with the Poisson model and with a linear model in each of these populations.
In general, the estimates of parameters and the predictions of random effects obtained with the Poisson model were better than their linear model counterparts. This is. not surprising because, to begin with, in the Poisson analysis, the data were analyzed with the models employed to generate the records in the simulation. Further, if the distribution is indeed Poisson, there is no appropriate variance decomposition independent of location parameters in the linear model, which makes the variance component estimators employed with normal data somewhat inadequate for estimating dispersion parameters, particularly in the mixed effects model. A log-transformation improved (relative to untransformed REML) the mean squared error performance of REML estimates of genetic variance, but worsened the estimates of residual variance. Similarly, fixed effects were estimated by BLUE with a larger bias and mean squared error than when estimated by the Poisson model.
BLUP was found to be robust in predicting random effects, although Poisson joint modes were significantly better with respect to bias and MSE in many cases. This is remarkable, considering that the BLUP predictors were based on location-parameter-dependent estimates of variance. BLUE, however, showed bias with increasing values of fixed effects and variance components. Truncated-Poisson estimators and predictors always outperformed BLUP and REML when zero counts were excluded from the analysis; truncation of null counts may be common in field records on traits such as litter size (Perez-Enciso et al, 1993).
Estimates of variance components obtained by MML in Poisson mixed models did not exhibit the typical bias due to small subclass sizes encountered often in threshold models. Rather, the downward biases of MML found in the one-way models may be due to small amount of statistical information. After all, MML is, like REML, a biased estimator. However, it should be consistent, because all Bayesian estimators are so under certain forms of selection (Fernando and Gianola, 1986). H6schele et al (1987) attributed the 'small subclass' bias in threshold models to inadequacy of a normal approximation invoked in MML. The same approximation is made when estimating variance with the Poisson model, but its consequences may be less critical here.
In conclusion, REML and BLUE do not perform well when the assumption of a Poisson distribution holds, even if all pertinent factors in the model are included in the analysis. BLUP, however, was found to be robust, although it had a slightly inferior sampling performance. This study did not address the situation when the data are counts, but the distribution is not Poisson. For this instance, estimators based on the assumption of normality may perform better than those that rely on the Poisson distribution, due to central limit theory. It should also be noted that if the operational Poisson model fails to consider all pertinent explanatory factors, the conditional variance may be substantially greater than the conditional mean of an observation; these 2 parameters are defined to be equal as in !2!. Possible diagnostics for investigating this source of overdispersion were discussed by Dean (1989).