A simulation study for the analysis of uncertain binary responses: Application to first insemination success in beef cattle

A simulation was carried out to investigate the methods of analyzing uncertain binary responses for success or failure at first insemination. A linear mixed model that included, herd, year, and month of mating as fixed effects; and unrelated service sire, sire and residual as random effects was used to generate binary data. Binary responses were assigned using the difference between days to calving and average gestation length. Females deviating from average gestation length lead to uncertain binary responses. Thus, the methods investigated were the following: (1) a threshold model fitted to certain (no uncertainty) binary data (M1); (2) a threshold model fitted to uncertain binary data ignoring uncertainty (M2); and (3) analysis of uncertain binary data, accounting for uncertainty from day 16 to 26 (M3) or from day 14 to 28 (M4) after introduction of the bull, using a threshold model with fuzzy logic classification. There was virtually no difference between point estimates obtained from M1, M3, and M4 with true values. When uncertain binary data were analyzed ignoring uncertainty (M2), sire variance and heritability were underestimated by 22 and 24%, respectively. Thus, for noisy binary data, a threshold model contemplating uncertainty is needed to avoid bias when estimating genetic parameters.


INTRODUCTION
Until recently, limited resources had been allocated to the research of fertility measures in beef cattle. This was, in part, due to the limited availability of data and computational requirements for traits of this nature. Two sources of information, natural service (NS) and artificial insemination (AI) mating data are generally associated with fertility traits. Few studies, however, have looked at the possibility of combining these two sources of information for the purpose of genetic evaluation. Furthermore, these two sources of information have not previously been combined due to the absence of a trait that could be used to describe fertility in both NS and AI mating data. Donoghue et al. [5] suggested a binary trait, calving to first insemination (CFI), which evaluated the probability that a calving event occurred as the result of conception at first insemination. For NS matings, success at first insemination was defined as conceiving during the first 21 days of the breeding season, which corresponded to the first estrous cycle of the female.
A binary trait such as CFI, for genetic evaluation of combined NS and AI mating data, as opposed to a continuous trait such as days to calving (DC), is dictated by the nature of the recording scheme currently used in the beef cattle industry; DC is defined as the time elapsed between the introduction of the bull and the subsequent calving date. In an NS mating situation, the date at which the female is exposed is recorded but the conception date is unknown. In contrast, AI data provides an accurate insemination date but the first exposure is uncertain. Yet, Donoghue et al. [5] proposed that NS and AI data could be combined in one analysis for the trait of CFI due to the high genetic (0.82) correlation between NS and AI mating data for CFI. Although this result is encouraging, it does not alleviate the uncertainty of CFI in the NS data.
With AI data, the exact day of insemination is known whereas in NS data, the number of days to insemination (NDI) resulting in conception is unknown. In order to derive NDI (the time interval between the introduction of the bull and insemination resulting in conception), knowledge about DC and the length of the estrous cycle is used. Since the gestation length (GL) in NS mating data is unknown, an average GL, obtained from AI mating data is deviated from DC to give the approximate NDI. This procedure lends itself to the possibility of uncertainty when females deviate in GL. For example, if a female has an NDI of 20 days, but has a longer than average GL, she would have been classified as failing to conceive at first insemination when the opposite (successful conception at first insemination) was true; thus leading to uncertain binary responses of success or failure at first insemination.
The current study was motivated by the belief that uncertainty is present in NS mating data for a trait such as CFI (at an unknown rate) and to date there have been no studies addressing this particular issue. One proposed way to account for uncertainty in success or failure at first insemination is via a threshold model using fuzzy logic classification, where fuzzy logic classification uses imprecise propositions based on fuzzy set theory to assign partial membership of a set [2]. Therefore, the objective of the current study was to evaluate different methods of analyzing uncertain binary responses with application to success or failure at first insemination in beef cattle.

Statistical analysis and computations
Threshold models are becoming a standard tool for the analysis of discrete data in the field of animal breeding and genetics. Extensive literature on its theoretical basis, implementation and application has been generated in the last twenty years [6,7,12]. More recently, Rekaya et al. [10] have proposed a method for analyzing binary data subject to misclassification using a threshold model. In the present study, an extension of such a method, based on fuzzy logic classification, is presented.

Threshold model for analysis of uncertain binary responses
Let m = (m 1 , m 2 , ..., m n ) be a sample of uncertain independent binary responses of non-observed real data y = (y 1 , y 2 , ..., y n ) , where each y i is Bernoulli with success probability p i that is expressed as a function of some systematic and random effects. Miscoding occurs if some y i is switched, e.g., y i = 0 becomes m i = 1 (i.e., a zero coded as one). Following notation by Rekaya et al. [10], let α = (α 1 , α 2 , ..., α n ) , where α i is an indicator variable for observation i such that α i = 1 if y i is switched, α i = 0 otherwise. Suppose each α i is a Bernoulli with success probability π t (probability of uncertainty) at time where π t is specific for each observation and will be explained further in Section 2.1.2.
Assuming independence between α and y, their joint distribution, given θ and π, is where θ is the vector of systematic and random effects and π = (π t 1 , π t 2 , ..., π t n ) . Thus, where p i (θ) = Φ i (θ), and Φ i (·) is the standard normal cumulative distribution function. Furthermore, the following relationship between y i and m i , given α i , could be established as such that y i = m i when there is no miscoding of the binary response. Using the relationships in equations (2) and (3), the joint probability of α and m, given θ and π, is Finally, prior distributions for θ and π would complete the Bayesian formulation. However, in some situations, π is known or could be inferred from external information. In the next section, fuzzy logic classification was used to determine the vector π.

Application to first insemination success in beef cattle
In beef cattle when natural mating is used, there is no direct recording of conception within the first 21 days of the breeding season. The first 21 days of the breeding season correspond to the length of the estrus cycle in the cattle. Therefore, we are interested in discriminating between cows that conceived in the first 21 days after introduction of the bull -conception occurring during the first estrus cycle -and those that did not conceive in the first 21 days. The only information available from NS matings is DC, which is computed as the time elapsed between the introduction of the bull and the subsequent calving date. Success or failure at first insemination (FIS) is based on the difference between DC and an average GL, where the average GL differed by sex of the calf. Thus, if the difference between DC and average GL is less than or equal to 21 days then FIS = 1, otherwise FIS = 0; in other words, a successful FIS event would be defined as conception occurring during the first 21 days of the breeding season.
Given the variation in GL between cows, and potentially of the estrus cycle, it is possible that some cows have uncertain FIS. Furthermore, the variation in GL, estrus cycle, and the difference between DC and average GL could be used to assess the probability of miscoding for every FIS record. In this simulation study, given the limited information on the variation of the estrus cycle and for the sake of simplicity, only variability in GL was used to assess the uncertainty of FIS. However, for real data sets and depending on the available information, sources of variation other than GL can be used within the framework of the presented method in a straightforward manner.
In the current study, fuzzy logic classification, based on the binary response of FIS and the difference between DC and average GL, was used to calculate the probability of miscoding at time t i . The deviation of average GL from DC gives the approximate NDI (the time interval between the introduction of the bull and insemination resulting in conception). It was then assumed that if the absolute value of the deviation of NDI and 21 was smaller (greater) than a number of days (based on practical and biological knowledge), then the observed FIS response was certain; in other words, if only one standard deviation (5 days) of GL is used to assess uncertainty, cows having an absolute difference between DC and average GL of less than 16 days or greater than 26 days, then there was no uncertainty about the observed FIS response. Otherwise, the binary record was considered as uncertain.
To assess the uncertainty about FIS, a fuzzy logic approach was used. There is no analytical method to determine a fuzzy function. However, such functions have to be derived based on practical and biological knowledge about the data generating process and good heuristics. Given the subjective nature of fuzzy classification, three fuzzy functions were used to assess their validity and the robustness of the method to fuzzy function choice. The following three fuzzy logic functions were used to compute the miscoding at time t i (see Fig. 1): and The interval in equations (5a-c) was based on one standard deviation of five days for GL obtained from AI data. Thus, the interval of uncertain FIS responses ranged from day 16 to 26 after introduction of the bull. In order to evaluate the robustness of the procedure, the interval of uncertain binary responses was increased from 10 to 14 days (nearly 1.5 standard deviation units) and the fuzzy functions were adjusted accordantly (see Fig. 2): and Although the interval of uncertain FIS could be increased to more than 1.5 standard deviation units (cows with very short or long GL), it is biologically reasonable to assume that the number of cows with very short (long) GL is small. Furthermore, increasing the interval of uncertain NDI days to more than 1.5 standard deviation units will unnecessarily increase the number of potential uncertain responses and may not be the best approach. A mixed linear model was used for the analysis of the underlying liability of FIS. In matrix notation the model could be written as: where λ is a vector of unobserved liabilities, β is the vector of fixed effects, s is the vector of unrelated random service sire effects, u is the vector of sire effects, and e is the vector of residual effects. Furthermore, X, Z s , and Z u are the corresponding incidence matrices with the appropriate dimensions.
Using the same notation as in the previous section, let m be defined as the vector of observed uncertain FIS responses and y be defined as the vector of unobserved true FIS responses. The joint probability of α and m, given θ = (β , s , u ) and π, is identical to equation (4) such that where p i (θ) = Φ i (x i β + z si s + z ui u) is the probability of success or failure at first insemination for record i and x i , z si , and z ui are the known row vectors relating the fixed, service sire and additive sire effects to the probability of first insemination success, respectively.
To ensure proper posterior distribution, the following priors were assumed for the parameters in the model: and where I is the identity matrix and A is a known matrix of relationships between sires. Uniform bounded priors, U[0, 1], were assumed for σ 2 s and σ 2 u . The joint posterior density was obtained by the product of densities in equations (7) through (10) leading to defined only within the boundaries of the bounded priors for σ 2 s and σ 2 u . Following Rekaya et al. [10], the conditional posterior distribution of each α i was given by and Thus, each α i (i = 1, 2, ..., n) was sampled from a discrete distribution, with probabilities as in equations (11) and (12). Once all α i (i = 1, 2, ..., n) had been sampled, the true data, y, was generated using the relationship in equation (3). The joint posterior density of θ, σ 2 s , and σ 2 u , given the true data, was obtained as follows: The joint posterior density in equation (13) was augmented with the liabilities for all observations in the data set. All conditional posterior distributions of model parameters were in closed form as described by Albert and Chib [1] and Sorensen et al. [12]. These distributions were normal for the location parameters, truncated normal for each of the liabilities and scaled-inverted chi-square distributions for the dispersion parameters. Liabilities were sampled from their truncated normal distribution using the inverse cumulative distribution function technique [4].
For each round of Gibbs sampling, α i was sampled, taking the value of one or zero. At completion of the sampling process, the probability that observation i was uncertain could be estimated as the ratio between the number of times that α i was one and the total number of samples.

General simulation
A simulation, using a sire model, was carried out to investigate two methods of analyzing uncertain binary responses for FIS. Five over-lapping generations were simulated. The base population included forty unrelated sires and subsequent generations consisted of 115 sires each. Thus, a total of 500 sires were generated. The data set consisted of 10 000 daughter records from bulls in generations 2 through 5. Fixed and random effects were randomly assigned to each record. In all cases, no extreme case (ECP) classes were allowed. A more detailed explanation of ECP classes and the effects on variance component estimation can be found in Moreno et al. [8]. A summary description of the simulated data and pedigree is presented in Table I. A linear mixed model at the liability scale which included herd, year, and month of mating as fixed effects; and unrelated service sire, sire and residual as random effects were used to generate the binary responses. Herd and year/month of mating effects were drawn from uniform distributions, U[0.4, 0.9] and U[-0.4, -0.3], respectively. Service sire effects were drawn from a normal distribution, N(0, Iσ 2 s ), where I is the identity matrix and σ 2 s = 0.03. Transmitting abilities of the base population sires were sampled from a normal distribution, N(0, Iσ 2 u ), where I is the identity matrix and σ 2 u = 0.02. The remaining sire effects were sampled from a normal distribution with the mean equal to one half of the grandsire's transmitting ability and variance equal to (3/4)σ 2 u . The residual variance was sampled from a normal distribution, N(0, Iσ 2 e ), where I is the identity matrix and σ 2 e = 1.0. The liability for each daughter record was calculated as the sum of all effects included in the model. Given the liability, the cycle in which the daughter conceived was assigned. A daughter was randomly assigned to one of five cycles, where cycle 1 is the first 21 days after introduction of the bull (i.e., the first estrous cycle occurring after the bull was introduced), cycle 2, 3, and 4 were the second, third and fourth 21 days, respectively, after the introduction of the bull, and cycle 5 was non-calvers (failed to calve). The NDI was randomly sampled for cows conceiving in cycles 1 to 4 using the integer value in the following formula: where c is the cycle and u was drawn from a uniform distribution, U[0, 1]. The NDI obtained in equation (14) was then used to generate two data sets. The first data set consisted of binary responses (D1) for FIS assigned using simulated NDI, where the binary responses were set to one if NDI was less than or equal to 21 days, and zero otherwise (including non-calvers). The second data set (D2) consisted of DC values. The DC values were obtained by generating a GL from a normal distribution with mean (standard deviation) of 280 (5.1) days for bull calves and 279 (5.2) days for heifer calves and adding the simulated NDI from equation (14). Sex of calf was randomly assigned and non-calvers (cycle 5) were assigned a subjective DC value equal to 980 days.
To mimic the field data, DC records in D2, assuming that NDI and GL were unknown, were used to generate the binary trait of FIS. For that purpose, an average GL specific to each sex (280 and 279 days for male and female calves, respectively) was subtracted from the DC in D2. If the difference was less than or equal to 21 days, then FIS = 1, otherwise FIS = 0 (including non-calvers). This third data set was labeled D3.
Three methods of accounting for uncertainty in success or failure at first insemination were investigated: (1) analysis of D1 (certain data -no uncertainty present) using a threshold model (M1); (2) analysis of D3 (uncertain data), without consideration of uncertainty, using a threshold model (M2); and (3) analysis of D3 (uncertain data), using a threshold model with fuzzy logic classification (M3 and M4). Three different fuzzy logic functions were used to account for uncertain FIS responses. For M3, the use of fuzzy function in equations (5a), (5b), and (5c) to account for uncertainty will be denoted as M3-a, M3-b, and M3-c, respectively. Similarly, the use of equations (6a), (6b), and (6c), in M4, will be denoted as M4-a, M4-b, and M4-c, respectively. Ten replicates were simulated for analysis using M1, M2, M3a-c, and M4a-c.

Field data based simulation
In order to evaluate the validity of the proposed method with field data, a real data based simulation, using a sire model, was carried out. The structure of a true available data set after stringent editing was used for the simulation. Editing consisted in removing herds and service sires with less than five records as well as daughters with unknown sires. The resulting data set consisted of 24 961 daughter records from 1 667 sires. A description of the field data structure and pedigree is presented in Table I. A linear mixed model which included herd, year, and month of mating as fixed effects; and unrelated service sire, sire and residual as random effects were used to generate the binary responses. Herd, year, and month of mating effects were drawn from uniform distributions, U[0.10, 0.60], U[-0.35, -0.25], and U[-0.40, -0.30], respectively. Service sire, sire, and residual effects were sampled in the same manner as the original simulation in Section 2.2.
A daughter was randomly assigned to one of five cycles using the sum of all effects included in the model, where the five cycles were defined as in Section 2.2. The NDI was randomly sampled using equation (14) and three data sets (D1, D2, and D3) were generated as in Section 2.2. The simulated data was then analyzed using the three methods (M1, M2, and M3) exactly in the same manner as in Section 2.2. For M3 and given the lack of differences between the three fuzzy functions used in the first simulation, only one fuzzy function (5a) was used to account for uncertainty with the real data based simulation and will be denoted as M3a. Ten replicates were simulated.

Convergence
Convergence diagnostics were based on the method of Raftery and Lewis [9] as implemented in the Bayesian Output Analysis (BOA) program [11]. The required burn-in period was always less than 2000 iterations for all parameters and all analyses. Thus, a total chain length of 75 000 iterations of the Gibbs sampler was run with a conservative burn-in of 25 000 iterations. The remaining 50 000 iterations were retained without thinning for post-Gibbs analysis. Point estimates of the mean and standard deviation were obtained from the BOA software [11] for all parameters and analyses. Furthermore, the high posterior density 95% [HPD (95%)] intervals were calculated using the algorithm described by Chen and Shao [3] as implemented in the BOA program [11].

Variance components
The posterior mean, standard deviation and the HPD (95%) interval for service sire and sire variance, based on 10 replicates, is presented in Table II. The point estimates of service sire (0.033) and sire (0.018) variance using M1 were, as expected, similar to the true values (0.03 and 0.02, respectively) used in the simulation. For M2, point estimates of the service sire (0.026) and sire variance (0.014) were underestimated compared to the true values of these parameters and the point estimates using M1 (0.033 and 0.018, respectively). Furthermore, the true values were barely included in the HPD (95%) interval of M2 indicating a potential bias in the estimation of genetic parameters. The point estimate for sire variance obtained from M3a-c as well as M4a-c (ranging from 0.018 to 0.020) were virtually the same as M1 (0.018), whereas the service sire estimates were slightly smaller (ranging from 0.028 to 0.029) when compared with M1 (0.033). For both parameters, the point estimates were similar to the true values used in the simulation.

Heritability estimates
The posterior mean, standard deviation and the HPD (95%) interval for heritability, based on 10 replicates, is presented in Table III. The point estimate of heritability obtained using M1 (0.070) was slightly lower than the true value (0.076) used in the simulation. The estimate of heritability obtained from M2 (0.053) was clearly underestimated. The true value of heritability (0.076) was barely included in the HPD (95%) interval similarly to the variance components. Using M3a-c and M4a-c, the point estimates of heritability (ranging from 0.070 to 0.0740) were virtually the same as the true value used in the simulation and the point estimate obtained using M1 (0.070).

Pearson correlations
The results of the Pearson correlations for service sire effects and sire breeding values were averaged over the 10 replicates and are presented in Table IV. Pearson correlations between true breeding values and predicted breeding values obtained from M1 and M2 were 0.50 and 0.45, respectively. When the model contemplated uncertainty using fuzzy logic classification (M3a-c and M4a-c), Pearson correlations were 0.46 between true and predicted breeding values. For service sire effects, Pearson correlations between true and estimated service sire effects were 0.61 and 0.58 for M1 and M2, respectively. Pearson correlations were 0.58 between true and estimated service sire effects when the model (M3a-c and M4a-c) contemplated uncertainty using fuzzy logic classification.  ; and the value in parentheses is the true value used in the simulation. c M1 = analysis of binary data using a threshold model; M2 = analysis of uncertain binary data, without consideration for uncertainty, using a threshold model; M3 = analysis of uncertain binary data, accounting for uncertainty from day 16 to 26 after introduction of the bull, using a threshold model with one of three fuzzy logic classification functions (functions a, b or c); and M4 = analysis of uncertain binary data, accounting for uncertainty from day 14 to 28 after introduction of the bull, using a threshold model with one of three fuzzy logic classification functions (functions a, b or c).

Variance components
The posterior mean, standard deviation and the HPD (95%) interval for service sire and sire variance, based on 10 replicates, is presented in Table IV. The point estimates of service sire (0.033) and sire (0.018) variance using M1 were, as expected, similar to the true values (0.030 and 0.020, respectively) used in the simulation. The point estimates for sire and service sire variance obtained from M3a were 0.017 and 0.030, respectively. The estimates obtained from M3a were quite similar to the corresponding estimates obtained from M1 (0.018 and 0.033, respectively). For both parameters, the point estimates from M1 and M3a were similar to the true values used in the field data based simulation. Similar to the general simulation, point estimates of the service sire (0.023) and sire variance (0.013) obtained from M2 were underestimated compared to the true values of these parameters and the point estimates using M1  Table IV. Summary of the posterior mean (PM), standard deviation (PS), and lower (HL) and upper (HU) bounds of the high posterior density 95% interval for variance components obtained using the field data based simulation a . (0.033 and 0.018, respectively). Furthermore, the true value of service sire and sire variance were barely included in their respective HPD (95%) intervals using M2 indicating a potential bias estimating these parameters.

Heritability estimates
The posterior mean, standard deviation and the HPD (95%) interval for heritability, based on 10 replicates, is presented in Table V. Similar to the general simulation, the point estimate of heritability obtained using M1 (0.069) was slightly lower than the true value (0.076) used in the simulation. The estimate of heritability obtained from M2 (0.051) was clearly underestimated. In fact, the true value of heritability (0.076) was similar to the upper bound of the HPD (95%) interval. Using M3a, the point estimate of heritability (0.065) was slightly smaller than the true value used in the simulation, but similar to the point estimate obtained using M1 (0.069).

Pearson correlations
The results of the Pearson correlations for service sire effects and sire breeding values were averaged over the 10 replicates and are presented in Table VI for the simulation using a field data structure. Pearson correlations between ; and the value in parentheses is the true value used in the simulation. c M1 = analysis of binary data using a threshold model; M2 = analysis of uncertain binary data, without consideration for uncertainty, using a threshold model; and M3a = analysis of uncertain binary data, accounting for uncertainty from day 16 to 26 after introduction of the bull, using a threshold model with fuzzy logic classification function a. true and predicted breeding values obtained from M1, M2, and M3 were 0.40, 0.36, and 0.37, respectively. The Pearson correlations between true and estimated service sire effects were slightly higher than the correlations between breeding values. The correlations between true and estimated service sire effects obtained from M1, M2, and M3 were 0.46, 0.40, and 0.41, respectively.

Parameter estimation
The results of this study indicate that a statistical model that contemplates uncertainty of binary responses (M3a-c and M4a-c) resulted in more accurate estimation of variance components and heritability when compared with a standard threshold model (M2). There was virtually no difference between the point estimates obtained from M3a-c and M4a-c with the true values of the parameters (Tabs. II and III). Furthermore, point estimates of genetic parameters using the correct binary data (M1) and the uncertain data analyzed with fuzzy logic classification (M3a-c and M4a-c) were practically the same.
When uncertain binary data were analyzed with a standard threshold model with no consideration for uncertainty (M2), sire variance, service sire variance, and heritability were underestimated by 22, 21, and 24%, respectively, when compared with M1. The service sire variance using M2 was slightly less Table VI. affected by the uncertainty of the binary responses. This result was expected because of the large number of records per service sire (on average, 33 records per service sire). In fact, the true value for service sire variance (0.030) was well centered within the HPD (95%) interval for M2. The true values for sire variance (0.020), and heritability (0.076), however, were closer to the upper bound of the corresponding HPD (95%) interval using M2. Based on the normal like shape of the posterior distributions of these parameters, this result indicates that the estimates of sire variance and heritability using model M2 are less likely to be equal to the true values. Thus, point estimates of these parameters will underestimate the true values and could be biased. Conversely, the true value of sire variance, service sire variance, and heritability were well centered within the HPD (95%) interval using M1, M3a-c, and M4a-c; furthermore indicating that a threshold model contemplating uncertain binary responses resulted in a more accurate estimation of genetic parameters.
Similar results were obtained when the simulation was conducted following the structure of an existing field data set. In fact, the statistical model that contemplated uncertainty of binary responses (M3a) resulted in a more accurate estimation of variance components and heritability when compared with a standard threshold model (M2). There was very little difference between the point estimates obtained from M3a with the true values of the parameters (Tabs. IV and V). Furthermore, point estimates of genetic parameters using the correct binary data (M1) and the uncertain data analyzed with fuzzy logic classification (M3a) were practically the same. However, when uncertain binary data were analyzed with a standard threshold model with no consideration for uncertainty (M2), sire variance, service sire variance, and heritability were underestimated by 28, 30, and 26%, respectively, when compared with M1. Sire and service sire variances using M2 were affected by the uncertainty of the binary responses. The underestimation of the genetic parameters was more pronounced in the field data case in part due to the smaller number of records per service sire and sire (on average, 16.6 records per service sire and 15 records per sire) when compared with the general simulation (on average, 33 records per service sire and 21.7 records per sire). Furthermore, the true value for service sire variance (0.030) was barely contained within the HPD (95%) interval for M2. The true value for sire variance (0.020), and heritability (0.076), however, were similar to the upper bound of the corresponding HPD (95%) interval using M2. These results suggest that the estimates of the genetic parameters using model M2 are likely to be biased.
The results of this study suggest that when analyzing binary data with uncertainty, a standard threshold model could lead to biased inferences. The problem could be avoided, or at least attenuated, by a statistical model that contemplates uncertainty. Moreover, for field data with low heritability, the use of external information is helpful in discriminating between uncertain, which is coded correctly or incorrectly, and miscoded data. In this study, the mean and standard deviation of the GL were used to assess the probability of miscoding of success or failure at first insemination via fuzzy logic classification. There did not appear to be large differences between the intervals (14 to 28 days vs. 16 to 26 days) or fuzzy logic functions (two linear and one non-linear) used to account for uncertainty of success or failure at first insemination. Furthermore, if extra information is available on potential other sources of uncertainty of FIS, such information should be included in the fuzzy function.

Pearson correlations
Although there were no observable differences in the Pearson correlations (Tab. VI) between true and estimated effects using M1, M2, and M3 in the current study, we would expect differences to be more noticeable in field data when the numbers of records per service sire or sire is limited. In such a case, a change in the status of one observation could have a major impact on the estimated effects. In NS mating data for small herds, it would not be unreasonable to expect service sires and sires to have fewer numbers of records per sire than those observed in the current study.
Although there were very small differences in the Pearson correlations (Tab. VI) between true and estimated effects using M2 and M3a using a field data structure in the simulation, the correlations were smaller than those obtained from the general simulation. As with the parameter estimation, this could be due to the fact that the number of records per service sire or sire is limited. Thus, a change in the status of one observation could have a more important impact on the estimated effects.