A decision problem

A decision problem is defined in terms of an outcome space, an action space, and an LF (Dawid 2007). Broadly defined, the goal of decision theory is to help choose among actions whose consequences cannot be completely anticipated, typically because they depend on some future or unknown state of the world; in our framework, the actions are the choices of the offspring that will be the parents in the next breeding cycle, which depend on the unknown BVs to be predicted. Expected loss theory handles this choice by assigning a quantitative loss to each decision, a probability to each state of the world, and then selecting an action that minimizes the expected value of the resulting loss. This idea has proven to be a widely applicable description of rational behavior (Parmigiani and Inoue 2009).

Let A be the action space, and be an action. No action can be taken without potential losses. The LF is denoted by and represents the associated penalty when a decision maker takes action , and the real state of nature is . In a Bayesian statistical framework, the expected loss is the expectation of the LF with respect to the posterior measure, i.e., When comparing two actions, and , after data X have been observed, the preferred action is the one for which the posterior expected loss is smaller. An the action that minimizes the posterior expected loss is called a Bayes action.

Loss functions in GS

In our particular case, the action space comprises the set of candidate lines for selection; we expect to assign a loss given our preferences, and then select the best individuals based on the expected loss. To explain the idea of selection under the decision theory approach, we first focused on a single quantitative trait where three LFs are proposed. In addition, we focused on traits whose phenotypic values we wanted to increase in successive selection cycles. Figure 1 illustrates the idea of selection under LFs for a univariate trait; Figure 1a depicts the distribution of the base population that is truncated at , meaning that the selected individuals will have mean . Parameters and are essential for defining the LF. Figure 1b depicts the three LFs we will describe later; these three functions are minimized when the means of selected candidates is close to , which reflects breeders’ preferences for those individuals with higher response to selection. The gray lines show the theoretical distribution of the trait of interest for individuals according to the breeder’s preference, as well as three possible distributions of candidates. The distribution of the selected candidate is in red because its mean is the closest to (i.e., higher response to selection) with almost the same variance of the base population; that is, under normality, the divergence between the truncated distribution and the distribution of the selected candidate is the lowest. This idea can be generalized to a multi-trait setting. Below we describe three LFs that we used when selecting candidates to be parents of the next generation.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

A) Classic idea of selection by truncation at . Loss functions need to be minimized at in order to favor lines with high response to selection. Losses were standardized by subtracting the minimum value for representation. Loss functions Kullback-Leibler (KL) and Continuous Ranked Probability Score (CRPS) are symmetric on both sides of target , while LinLin loss is asymmetric. B) The solid line in black represents the base population, the solid gray line corresponds to the truncated distribution after censoring at representing the breeder’s preferences. Dashed lines are theoretical distributions of three possible candidates. The candidates’ distribution with mean close to the theoretical (the greater ) and variance similar to that of the parent distribution has the minimum loss.

Univariate Kullback-Leibler (KL) loss function

The KL divergence reflects how different two probability distributions are. The expected loss of information will be minimal if the two distributions approach each other, and it will be zero if they are identical. Therefore, we can measure how close the distribution of the candidates’ BVs is to the hypothetical parental distribution.

Let be a random variable that represents the phenotypic value of a trait of interest in the base population and ; if selection is by truncation (i.e., those realizations of Y falling above a threshold value will be the parents of the next breeding cycle), then is the truncated random variable resulting from left truncation at , which represent the phenotypic values of the selected parents. From properties of normal distribution, , where denotes a truncated normal distribution with parameters and b; but for simplicity, let’s denote . The probability density function of is , for , and . Let be a random variable that represents the candidates’ phenotypic values. If we assume normality, then . Note that, in order to keep the heritability of the trait, we are assuming that Y and have different means but the same variance. A measure of the divergence between the distributions of and is given by the Kullback-Leibler divergence,(1)which we want to minimize in the breeding selection framework.

After evaluating the integral in (1) and simplifying (see Appendix A), the KL loss is given by:(2)where is the mean of the truncated normal distribution and equal to and φ and Φ denote the pdf (probability density function) and cdf (cumulative distribution function) of a standard Gaussian random variable, respectively. By reordering some terms, Equation (2) can be alternatively formulated as:(3a)(3b)where is the selection differential, is the selection response, and is the selection intensity. The second term on the right-hand side of Equation (3a) implies that when R approaches S (increased genetic gains) and the selection intensity increases, the divergence between the truncated distribution and the candidate’s distribution decreases. That is, depends on the intensity of selection and is a decreasing function of (3b).

Multivariate Kullback-Leibler (KL) loss function

In GS there are usually many traits of interest for which multivariate normality is assumed. As in the single-trait case, we can compute the divergence between two multivariate distributions and capture all the uncertainty in multivariate BVs. Then, assuming multivariate normality for phenotypic values of the traits under consideration, the KL divergence accounts for the association between traits, which improves and simplifies the decision procedure, given that all traits share some degree of association and are rarely independent.

Let be the random vector of the phenotypic values of t traits of interest in the base population. We assume that , where is the vector of means in the parental population and P is a positive definitive variance-covariance matrix that captures the association between traits. Now, let be the random vector resulting from left truncation at . Then , where TMVN denotes the Truncated Multivariate Normal distribution. The pdf of is where the normalization factor is . Now, let be the random vector of the phenotypic values of a candidate, and assume that . Let’s define as the differential selection vector and as the vector of the selection response, where G is the genotypic covariance matrix. As in the univariate approach, the multivariate LF based on the multivariate KL divergence is given by (see Appendix A):(4a)(4b)Thus, as the phenotypic and genotypic covariance matrices tend to explain the same amount of variation and association between traits, , the divergence between the parents’ distribution and the candidates’ distribution tends to decrease. Note that the expression is the matrix equivalent of multi-trait heritability (the ratio of the genetic variance in the numerator and the phenotypic variance in the denominator); thus, when , the heritability of each trait is 1 and and .

As in the univariate approach, the multivariate KL loss contains the term that indicates that the joint probability of falling above has less penalty.

Note that P is expressed in terms of its inverse; therefore, those traits with high variance induce a lower penalty, and there is a compromise between the gain from selection and keeping the variability. KL divergence is an appropriate way of selecting because it not only captures the association between traits, but also assigns weights to each trait automatically, thus reducing subjectivism in GS. In addition, Equation (4b) implies that traits with higher variance will have priority in a GS program and will directly preserve the genetic variance as much as possible. Univariate and multivariate KL LFs are represented in Figure 1a (red lines) and Fig. B1a (Appendix B), respectively.

The Continuous Ranked Probability Score and its generalization, the Energy Score

The Continuous Ranked Probability Score (CRPS) function is discussed in detail in Gneiting and Raftery (2007) and Hersbach (2000). The CRPS is a metric that reflects the distance between two random variables in terms of their cumulative distributions. The CRPS in GS acts by choosing candidate lines whose (cumulative) BV distributions are as close as possible to the hypothetical parental distribution that reflects the breeder’s preferences. Nevertheless, the main difference between KL divergence and CRPS is that the latter is less restrictive with lines whose distributions do not perfectly match the distribution that reflects our preferences, as can be seen in Figure 1a (green line).

Keeping the normality assumption, let be a normal random variable with mean and variance ; there is a closed form of this LF (Appendix A), which in our context is defined as(5a)(5b)(5c)where φ and Φ denote the pdf and cdf of a standard Gaussian variable, respectively. Again, CRPS increases with the intensity of selection and is a decreasing function of heritability.

When the interest is in multiple traits, we can use the Energy Score (EnergyS), which generalizes the univariate CRPS (Gneiting and Raftery 2007; Székely and Rizzo 2013) and allows a direct comparison between multi-traits. The Energy Score in terms of a genomic selection goal is expressed as(6)where denote the Euclidean norm, and were previously defined, and denotes an independent random vector with the same distribution as , i.e., . In contrast to KL divergence, with the CRPS and the Energy Score we can avoid the normality assumption while holding the assumption that traits are on an interval scale. This is a big advantage given that frequently the traits of interest in GS are not Gaussian. The EnergyS is depicted in Fig. B1b (Appendix B).

An asymmetric loss function in univariate and multivariate settings

At this point, the previous LFs that we adapted to the GS problem are all symmetric. Therefore, the penalties on the right- and left-hand side of the target at the same distance are equal. Nevertheless, it is natural to think that in conventional and GS plant breeding, those values of greater than might be less penalized if the goal is to increase phenotypic values. Here is where an asymmetric LF plays a key role because it assigns very small penalties to those individuals whose BV distributions are to the right of the hypothetical parental distribution. One simple asymmetric LF is the linear-linear loss (LinLin), whose behavior is linear on both sides of the target; it has the term that induces different penalties (Berk 2011) and is the loss of the quantile regression. The LinLin LF is defined as(7)where ; then this LF is also a decreasing function of . LinLin loss is depicted in Figure 1a (blue line). Note that when , the LinLin loss reduces to a symmetric linear LF.

In the multi-trait setting, we define as the vector of deviations of from and the multivariate LinLin function is 0 when phenotypic and genotypic matrices are identical, . The Multivariate Asymmetric Loss Function (MALF) (Komunjer and Owyang 2011) is the generalization of LinLin LF. For GS purposes, it is expressed as(8)where is the Euclidean norm and τ controls the degree of asymmetry. In a simplified version with the , Equation (8) is expressed as(9)In the univariate case, letting , Equation (9) reduces to twice the LinLin LF, i.e., . Both univariate and multivariate LinLin losses favor selecting those lines with distributions similar to the hypothetical parental distribution (which reflects the breeder’s preferences) as much as possible. The MALF is depicted in Fig. B1c (Appendix B).

Posterior expected loss and their approximations

Each of the LFs proposed above can be used in conventional selection and/or GS to select the best lines according to our preferences. Below we briefly explain the formal approach used to evaluate posterior mean of the LFs proposed above. Given a prediction model (for example, Multi Trait Model (MTM) or Bayesian Ridge Regression (BRR) for a single trait), we have the joint posterior distribution,(10)where is the likelihood function, and is the prior distribution for θ, and elements of θ are an overall mean (μ), marker effects or any other effects (β), and the variance-covariance matrix (P) captures the association between traits (or , the phenotypic variance in a univariate framework), y is a matrix (or vector in a single-trait setting) containing phenotypic records, and X represents the incidence matrix according to the multi-trait or single-trait model.

Let be the random vector of the phenotypic values of a candidate with a covariate vector of markers and distribution . The posterior predictive density corresponding to is given by(11)Hence, given the observed phenotypic and genomic data (), whatever we can predict about the candidate is given by (11) and any choice of the candidates for selection must take into account the uncertainty described by .

According to the Bayesian decision theory, the optimal choice is one that minimizes the posterior expected loss. In this case, once the phenotypical values, y, and the matrix of covariables, X, are observed, we have to average the LF over all the unknowns in the model, that is, θ and the observable but unknown phenotypic value of the o-th candidate. Given that is computed as a function of μ, β and P (or ), hereinafter we will denote all LFs described above as . Thus, the posterior expected value of the LF for candidate o is given by(12)Then, for each o line, we have a posterior expected loss and we will select those lines with the lowest posterior expected losses.

The expected value of the LFs given in (2), (5a) and (7), as well as their generalizations in the multivariate context (4b), (6), and (9) are approximated by Markov Chain Monte Carlo integration by taking m samples from the joint posterior distributions for , (in the univariate context) or for vectors , (in the multivariate context).

Application of univariate and multivariate loss functions in a wheat dataset

To illustrate the application of univariate and multivariate LFs in GS, we used 320 spring wheat lines with records on four traits: grain yield (GY), thousand-kernel weight (TKW), Zn and Fe concentrations in the grain (GZnC and GFeC, respectively) from CIMMYT’s biofortification breeding program and genotyped with DaRT markers (Velu et al. 2016). All traits are positively but low correlated; for example, GY had a correlation of 0.21 with TKW but zero with the other traits, whereas GZnC and GFeC had a 0.26 correlation.

For univariate LFs, we used Bayesian Ridge Regression (BRR) and the BGLR R package (de los Campos and Pérez Rodríguez 2015) in order to obtain posterior distributions of parameters and quantities of interest. We selected the top 10% (32 lines) of the wheat lines whose posterior expected losses were the minimum using the three univariate LFs discussed above (KL, CRPS, and LinLin). Results were contrasted with the top 10% based solely on punctual predicted BVs (Std), i.e., by selecting individuals with the highest punctual predicted BVs. For multivariate LFs (KL, Energy Score and MALF), we fitted the MTM model (de los Campos and Grüneberg 2016). The aim here was to select the “best” lines whose performance across all traits is high. For asymmetric LFs (LinLin and MALF), we fixed the value of α (or τ) equal to 0.9 (according to selection pressure) in order to impose big penalties on those individuals far to the left of the target point; on the other hand, lower penalties were given to those that were equal or greater.

Application of univariate and multivariate loss functions in a simulation scheme

Simulations were used to evaluate the three univariate and multivariate LFs discussed above. In the univariate scenario, we compared results obtained with the selection based on the standard method (Std). In the multivariate scenario, we contrasted results of different LFs. A recurrent selection scheme was used in which the best individuals were selected based on their merit regarding a single quantitative trait or three correlated quantitative traits, and measured through the LFs.

Traits and heritabilities

In the univariate simulation, one trait was genetically simulated, where the gene effects were sampled from a gamma distribution with shape and scale parameters equal to two. Phenotypes were simulated by summing up all true genotypic values and adding a residual effect consistent with the expected heritability fixed at 0.5. That is, , where p is the number of genes, and is the genotype for the j-th gene of the i-line (coded as -1 and 1 for the two homozygotes, respectively). The environmental noise was sampled as , where is the genetic variance.

In the multi-trait simulation, three correlated quantitative traits were genetically simulated assuming a full pleiotropic model. This was carried out by randomly sampling gene effects for all segregating sites from a multivariate normal distribution with zero mean and a previously stated variance-covariance, to ensure a genetic correlation at first generation of -0.37 between trait 1 (T1) and trait 2 (T2); a genetic correlation of 0.34 between traits 2 (T2) and 3 (T3); and a genetic correlation of -0.02 between T1 and T3. To mimic complex and simple quantitative traits, narrow-sense heritabilities of 0.3 and 0.6 were assumed for all traits. Hereinafter, we will always refer to the traits in terms of narrow-sense heritability (), given that in the simulation scheme we simulated a purely additive model and did not include dominance effects.

Selection cycles

The forward-in-time component represented ten cycles of a classic recurrent selection scheme in which the breeder has the ability to select based on the genomic predictions of the breeding values.

In the multivariate LFs framework, genomic predictions of the evaluated lines were carried out using the Multi-Trait Model (MTM) (de los Campos and Grüneberg 2016), where 70% and 30% of the lines were used as the training and testing sets (candidates), respectively. For single-trait simulation, genomic predictions were made using Bayesian Ridge Regression (BRR).

For multi-trait simulation, the simulated phenotypic values were modeled assuming an intercept for each trait as a fixed effect and the predicted breeding value for each line in each trait as random. As in the wheat dataset, we fitted the same MTM model for the multi-trait simulation. In both the univariate and multi-trait simulation, posterior distributions of the predicted breeding values of the candidate lines were ranked according to each selection process (the univariate and multivariate LFs described above). For all LFs described above, approximations to the posterior expected loss were obtained considering 10,000 MCMC samples after a burn-in period of 30,000 samples. For univariate simulation, the top 10% and the top 30% of the lines in the candidate set were selected via their minimum posterior expected losses, while in the multi-trait simulation scheme, the top 10% of the lines in the candidate set were selected with minimum posterior expected losses.

In the LinLin LF and its generalization MALF, we fixed the value of α or τ equal to 0.9 or 0.7 (depending on whether the selection pressure was 10% or 30%) in order to impose big penalties on those individuals in the candidate set that were on the far left-hand side of the target point, while lower penalties were given to those that were equal or greater.

Results are presented as summaries of 20 replicates for each heritability, selection pressure and LF scenario. Individuals selected as parents will have the highest possible BV and the lowest possible value of the corresponding LFs. The simulations were implemented in a C++ program that was compiled, linked and executed within the R version 3.3.3 (R Development Core Team 2017) through the facilities provided by the Rcpp package (Eddelbuettel and Francois 2011).

Data Availability

The durum wheat data including clean imputed markers and phenotypic information for the traits can be found in the following link: http://genomics.cimmyt.org/Decision_theory_GS/