In the LR method, Legarra and Reverter (2018) used var(x) to denote the variance of a random element, x, sampled from a single realization of the random vector x, and Var(x) to denote the variance-covariance matrix for the vector x. Similarly, they used cov(x, y) to denote the covariance between randomly sampled x _i and y _i and Cov(x, y) to denote the covariance matrix between the vectors x and y. Let u denote the vector of BV of the validation animals, and u ̂ p and u ̂ w denote the vector of EBV of u obtained from partial data and whole data, respectively.

Legarra and Reverter (2018) developed the LR method for best linear unbiased prediction (BLUP) by showing that C o v ( u ̂ w , u ̂ p ) = V a r ( u ̂ p ) , using the results from Reverter et al. (1994) and the result that C o v ( u , u ̂ p ) = V a r ( u ̂ p ) and E ( u ̂ p ) = E ( u ̂ w ) = E ( u ) from Henderson (1982). Using the above results for BLUP, they showed that E [ c o v ( u ̂ w , u ̂ p ) ] = E [ v a r ( u ̂ p ) ] = E [ c o v ( u , u ̂ p ) ] , and thus justified that c o v ( u , u ̂ p ) can be quantified by computing c o v ( u ̂ w , u ̂ p ) .

It is well known that whenever u ̂ p is the conditional mean of the unobservable random variable u _p given observed phenotype, C o v ( u , u ̂ p ) = V a r ( u ̂ p ) (Rao, 1973). Further, based on the double expectation theorem (DeGroot and Schervish, 2002), the condition mean is known to be unbiased in the sense that E ( u ̂ ) = E ( u ) . Thus it follows that, if the conditional mean is used for prediction, E ( u ̂ p ) = E ( u ̂ w ) = E ( u ) . Below, we show that whenever u ̂ p is based on conditional means, C o v ( u ̂ w , u ̂ p ) = V a r ( u ̂ p ) . Then, using the results from Legarra and Reverter (2018), it follows that E [ c o v ( u ̂ w , u ̂ p ) ] = E [ v a r ( u ̂ p ) ] = E [ c o v ( u , u ̂ p ) ] and c o v ( u , u ̂ p ) = c o v ( u ̂ w , u ̂ p ) for predictions using conditional means, regardless of whether the conditional mean is a linear or non-linear function of the data.

Here we show that C o v ( u ̂ w , u ̂ p ) = V a r ( u ̂ p ) when predictions are based on conditional means, which may be non-linear. Let y w = y p y r , where y _w , y _p , and y _r indicate vectors of phenotype records in the whole, partial, and validation (remaining) data set, respectively. In the Supplementary Appendix, we show E y r | y p ( u ̂ w | y p ) = u ̂ p , which represents the expectation over the conditional distribution of y _r given y _p and will be used in the following proof.

Now, we write the C o v ( u ̂ w , u ̂ p ) as: C o v u ̂ w , u ̂ p = E y w u ̂ w − θ w u ̂ p − θ p ′ = E y w u ̂ w u ̂ p − θ p ′ − E y w θ w u ̂ p − θ p ′ = E y w u ̂ w u ̂ p − θ p ′ − θ w E y w u ̂ p − θ p ′ = E y w u ̂ w u ̂ p − θ p ′ − θ w 0 ′ = E y w u ̂ w u ̂ p − θ p ′ = E y w u ̂ w − u ̂ p + u ̂ p u ̂ p − θ p ′ = E y w u ̂ w − u ̂ p u ̂ p − θ p ′ + u ̂ p u ̂ p − θ p ′ , where θ _p and θ _w are the expected values of u ̂ p and u ̂ w , respectively. The first term of this expectation can be shown to be null:

(see the Supplementary Appendix for derivation) E y w u ̂ w − u ̂ p u ̂ p − θ p ′ = E y p E y r | y p u ̂ w − u ̂ p u ̂ p − θ p ′ | y p = E y p E y r | y p u ̂ w − u ̂ p | y p u ̂ p − θ p ′ = E y p u ̂ p − u ̂ p u ̂ p − θ p ′ = 0 , because, as shown in the Appendix that E y r | y p ( u ̂ w | y p ) = u ̂ p . Thus, the C o v ( u ̂ w , u ̂ p ) becomes: C o v u ̂ w , u ̂ p = E y w u ̂ p u ̂ p − θ p ′ = E y p u ̂ p u ̂ p − θ p ′ = V a r u ̂ p . The proof of C o v ( u ̂ w , u ̂ p ) = V a r ( u ̂ p ) shows the LR method also holds for non-linear predictions. This proof is similar in principle to that provided by Belay et al. (2022), but we recognize that it is not limited to BLUP, as invoked in that study, but is applicable to any method of prediction based on the conditional mean (Fernando and Gianola, 1986), including for non-linear models.

To investigate the behavior of the LR method, three partitioning scenarios (Figure 1) were implemented: 1) between animals: phenotype records for days 70–500 of the first 500 individuals comprised the partial set and the phenotypes of the remaining 1,000 individuals were assigned to the validation set, 2) by age within animals: phenotypes for days 70–300 of all 1,500 individuals comprised the partial set and all 1,500 individuals and their phenotypes from days 301–500 were considered as the validation set, and 3) between animals and by age: phenotypes for days 70–300 of the first 500 individuals comprised the partial set and phenotypes for days 301–500 for the remaining 1,000 individuals were assigned to the validation set. The EBV of interest were those for body weight for individuals in the validation set across days, predicted based on the partial set for each of the three scenarios.

As described below, following Legarra and Reverter (2018), two statistics are used to check whether the model is adequate.

1. The first statistic is Δ ̂ = u ̂ p ̄ − u ̂ w ̄ , which is the difference between the mean EBV of individuals in the validation set estimated based on the partial and whole data sets, respectively. This is an estimator of Δ = u ̂ p ̄ − u ̄ . A deviation of Δ ̂ from 0 indicates the model is not adequate.

This estimator was referred to as an estimator of dispersion by Legarra and Reverter (2018) and of inflation by Belay et al. (2022), where 1 indicates no bias. Suppose animals are selected based on EBV to increase the values of a trait. Then if the true regression coefficient is less than 1, the BV of selected candidates is expected to be lower than their EBV, which indicates an upward bias of the EBV of the selected animals. On the other hand, when b ̂ w p is larger than 1, the BV of selected candidates is expected to be higher than their EBV, which indicates a downward bias of the EBV of the selected animals.

Provided the model is adequate, in the LR method, population accuracy is estimated as ρ ̂ p = c o v ( u ̂ w , u ̂ p ) v a r ( u ) ̂ × v a r ( u ̂ p ) , where v a r ( u ) ̂ refers to an estimate of the genetic variance of individuals in the validation set. This estimate was obtained by Gibbs sampling as: v a r ( u ) ̂ = 1 n trn ∑ j = 1 n trn u j 2 − 1 n trn ∑ j = 1 n trn u j 2 , where u _j refers to a sample of the BV of individual j from its posterior distribution in the training set and n _trn is the total number of individuals in the training set (Fernando et al., 2017).

In addition to Δ ̂ , b ̂ w p , and ρ ̂ p , we also calculated the Δ, b _wp , and ρ _p , which are the “true” values of these quantities by using the simulated values of u in place of u ̂ w in the formulas. Note that these “true” values can only be computed in a simulation study, and they are used here to study the performance of the LR method.

The means of Δ and Δ ̂ were calculated for each day of age across all animals in the validation set. These mean values were averaged across days within each replicate to test whether their mean was significantly different from 0 using a t-test. Similarly, true and estimated regression coefficients were averaged across days within each replicate to test whether their mean was significantly different from 1 using a t-test. Additionally, we tested the difference between c o v ( u ̂ w , u ̂ p ) and c o v ( u , u ̂ p ) by first calculating the mean difference across age for each of 20 replications and then used a t-test to determine if these differences were significantly different from zero.

Figure 3 shows the Δ and Δ ̂ for EBV of body weight for each day when the data were partitioned between animals. When the true model (BHGGM) was used, both the Δ and Δ ̂ were symmetrically distributed around 0 for each day, and their mean was not significantly different from 0 (p = 0.84 and = 0.37, respectively). In contrast, when the wrong model (BHQGM) was used, the mean of the Δ was significantly different from 0 (p < 0.001), but the Δ ̂ were symmetrically distributed around 0 for each day and their mean was not significantly different from 0 (p = 0.4).

Figure 4 shows the Δ and Δ ̂ for EBV of body weights for each day when the data were partitioned by age within animals. When the true model was used, both the Δ and Δ ̂ were symmetrically distributed around 0 for each day, and their mean was not significantly different from 0 (p = 0.10 and = 0.09, respectively). When the wrong model was used, the Δ and Δ ̂ were significantly different from 0 (p < 0.001 and = 0.002, respectively). Results for the partitioning between animals and by age were consistent with those in Figure 4 and are shown in Supplementary Figure S1.

Figure 5 shows the true and LR estimates of the regression coefficient of EBV for body weights for each day when the data were partitioned between animals. When the true model was used, both the true and estimated regression coefficients were symmetrically distributed around 1 for each day, and their mean was not significantly different from 1 (p = 0.75 and = 0.53, respectively). When the wrong model was used, the true and LR estimates of the regression coefficient were significantly different from 1 (p < 0.001). Results for the partitioning by age within animals and between animals and by age were consistent with those in Figure 5 and are given in Figure 6; Supplementary Figure S2, respectively.

In Figure 7, the true and LR estimates of prediction accuracy of EBV for body weights for each day when the data were partitioned between animals are presented. The LR estimates of prediction accuracy had a similar pattern as the true estimates of accuracy when using the true model but not when the wrong model was used. When partitioning the data by age within animals, the LR estimates of accuracy showed a similar pattern as the true estimates of accuracy, regardless of the model fitted (Figure 8). We also evaluated the difference between c o v ( u , u ̂ p ) and c o v ( u ̂ w , u ̂ p ) when fitting the true and wrong models for the three data partitioning scenarios (Table 1). There was a non-significant difference (p ≥ 0.74) between c o v ( u , u ̂ p ) and c o v ( u ̂ w , u ̂ p ) when the true model was fitted, but a significant difference (p ≤ 0.004) was observed for each scenario when the wrong model was fitted.

To generalize the LR method for linear or non-linear predictions, we presented a mathematical proof to justify using c o v ( u ̂ w , u ̂ p ) to quantify c o v ( u , u ̂ p ) even when the prediction is non-linear, provided that the fitted model is adequate. In our proof, we assumed that the partial data contains a subset of the phenotypes of the whole data. Belay et al. (2022) showed the LR method is also applicable to BLUP when the partial data contains a subset of the genotypes of the whole data. The proof presented in the current paper is similar in principle to that provided by Belay et al. (2022). Taken together, these two proofs show that the LR method is applicable to predictions based on the conditional mean, regardless of whether the data are partitioned by genotypes or phenotypes and regardless of whether the model is linear or non-linear.

We used simulated longitudinal data to investigate the ability of the LR method to determine whether the fitted model is adequate. In the context of our situation, the non-linear relationship is between the age of the animal and its weight, which makes how the data are partitioned into training and validation sets important. To explore how the strategy for partitioning the data into training and validation sets affects the ability of the LR method to determine whether the fitted model is adequate, three data partitioning strategies were used: between animals, by age within animals, and between animals and by age. We did not simulate selection in our data because we were comparing two models of growth rather than two models of inheritance. When the fitted model for growth is not adequate and the data are partitioned by age within animals, predictions from the training data would not be accurate even without selection. Thus, selection was not required to investigate properties of the LR method in this context. Further, in our analysis of the simulated data, the QTL were used as markers, for simplicity, but this does not influence the conclusions of our study either.

Below, we summarize the implications of the data partitioning strategies on the performance of the LR method, thereby providing guidelines for using the Δ ̂ and b ̂ w p to determine whether the fitted model is adequate. Significant deviations of the LR statistics from their expectations under the correct model were tested by replication of the simulations. In real data analyses, repeated k-fold LR can be used to test significance of model inadequacy using the two LR statistics. Alternatively, if a Bayesian method is employed for the LR analysis, the posterior probability of model inadequacy can be computed from a single partitioning of the data.

When the wrong model (BHQGM) was fitted and the data were partitioned between animals, the Δ was significant, but the Δ ̂ was not able to detect this inadequate model (Figure 3). Macedo et al. (2020) also observed that for a certain misspecification of the model, the LR method was not able to correctly detect and estimate Δ. Figure 5 shows that when the wrong model was used, the true estimate of regression (the regression of u on u ̂ p ) had a significant deviation from 1, and in this case the estimate of the regression coefficient based on the LR method was also significantly different from 1, although differing in magnitude from the true estimate of the regression coefficient. This is also consistent with the results observed by Macedo et al. (2020). The effect of fitting the true or the wrong model on the behavior of EBV against age for the partial and whole data were presented in Figure 2 (left column) for a randomly selected validation individual and data partitioning was by individual rather than age. In this case, when the wrong model was fitted, the EBV from the partial and whole data sets deviated more from the true BV than EBV from the true model did. However, even when the wrong model was used, the EBV from partial and whole data sets were very similar because the data were partitioned by individuals rather than age. This explains why the Δ ̂ based on the LR method was not significant when the wrong model was used. This is because of when the data were partitioned between animals, there is no difference in the distribution of phenotypic values in training and validation sets, because we did not simulate selection and partitioning was not by generation. Thus the mean of u ̂ p was not significantly different from mean of u ̂ w .

Recall that in the LR method, the regression of u on u ̂ p is estimated by the regression of u ̂ w on u ̂ p . In general, we expect that this statistic can be used to determine whether the model that is fitted is inadequate, in any study even when there is no difference in the distribution of phenotypic values in training and validation sets, such as when the data are partitioned by generation and there is no selection. The reason for this is that animals in the validation set will have a range of values for u, u ̂ p , and u ̂ w , and thus there is information to estimate the regression of u on u ̂ p . In fact, the variance of u in the validation set can be even higher without selection compared to with selection, where for example, the validation animals are all from the last generation of selection. In the LR method, these values for u ̂ w and u ̂ p are used to estimate the regression of u on u ̂ p . The variance of this regression estimator is proportional to 1 V a r ( u ̂ p ) , and thus, it is possible that this regression coefficient can be estimated even more accurately when there is no selection.

When the data were partitioned by age within animals, the LR method was able to correctly detect that the model was inadequate using Δ ̂ and b ̂ w p when the wrong model was used (Figures 4, 6). Figure 2 (middle column) shows the EBV of a randomly selected individual when the data were partitioned between animals. When the wrong model was fitted, the EBV estimated from partial and whole data sets both deviated from the true BV but the EBV based on the partial set was quite different from that estimated from the whole set. This illustrates the significant Δ ̂ that was detected by the LR method for this scenario. Results for the partitioning between animal and by age (right column in Figure 2) were similar to those when partitioning by age within animals.

The inconsistency between Δ ̂ and b ̂ w p for different data partitioning strategies suggests that the LR method captures different aspects of the model for different data partitions. When the data were partitioned between animals, both the partial and whole data sets included phenotypes over the range from 70 to 500 days. Thus the fit of a phenotyped individual’s growth curve based on the partial and whole data sets were similar, even for the wrong model, although the fit might deviate from that using the true model. Using the wrong growth model only provides an incorrect fit to the relationship between age and body weight within individual but does not affect the flow of information between relatives. Thus to appropriately test the predictive ability of the fitted growth model using the LR method, we had to predict the body weights of animals that are outside the observed age range for animals in the training set, as was done when the data are partitioned by age. In this case, the partial data set has only body weights measured at ages up to 300 days, whereas the validation data set has body weights measured at ages up to 500 days. Thus when we predict the body weights in the validation set based on the fit of the growth model from the partial data set, we are testing the predictive ability of the growth model, rather than the predictive ability of data from older animals to predict younger animals. In this case, the LR method was able to correctly detect that the model was inadequate when the wrong model was used (bottom right plots in Figures 4, 6).

In general, to properly test the model adequacy using the LR method, we need to use the model to predict the performance of individuals that have values for the relevant predictor variables or combination of predictor variables that were not present in the training data. In our simulated data, the predictor variables included the marker genotypes, as well as age. Let’s define the predicted performance of individual i as y i ̂ = f ( x i ; θ ̂ ) , where f(.) is the linear or non-linear function used for prediction, x _i is a vector of predictor variables for individual i, and θ ̂ is the vector of estimates of model parameters. Below we will use genomic prediction by ridge regression BLUP (RR-BLUP) as an example for illustration. To evaluate the predictive ability of RR-BLUP, the data are partitioned into training and validation sets. The training set is used to fit the predictive model f(.) and to estimate the model parameters θ (i.e., marker effects). By plugging the marker effect estimates θ ̂ and observed marker genotypes x into f(.), the performance of individuals in the validation set can be predicted. The predictive ability of the model is then quantified by comparing the predicted and observed performances of individuals in the validation set. In RR-BLUP, the relevant predictor variables are the marker genotypes. Thus the training and validation sets cannot include the same individual or individuals with the same genotypes across all loci. In our study, the LR method was used to determine whether it could detect an inadequate model when a wrong model was used for the growth curve, where the predictor is age and the dependent variable is BW (longitudinal BW). When the data were partitioned between animals, the training (partial) and validation sets included phenotypes for animals with age ranging from days 70–500, i.e., the same age range as used in the training data was used for the validation data and, therefore, Δ ̂ failed to detect the use of an inadequate model. However, when the data were partitioned by age within animals, the model was trained using phenotypes with ages ranging from days 70–300 and it was tested by predicting body weights for animals with age ranging from days 301–500. In this case, the LR method was able to detect model inadequacy using Δ ̂ . This was even true when the same genotypes were used in both the training and validation sets, because to check if the model used for predicting longitudinal BW is correct, the relevant predictor variable is age.