Title:	Probability Theory for Selecting Candidates in Plant Breeding
Version:	1.0.4.6
Description:	Use probability theory under the Bayesian framework for calculating the risk of selecting candidates in a multi-environment context. Contained are functions used to fit a Bayesian multi-environment model (based on the available presets), extract posterior values and maximum posterior values, compute the variance components, check the model’s convergence, and calculate the probabilities. For both across and within-environments scopes, the package computes the probability of superior performance and the pairwise probability of superior performance. Furthermore, the probability of superior stability and the pairwise probability of superior stability across environments is estimated. A joint probability of superior performance and stability is also provided.
URL:	https://github.com/saulo-chaves/ProbBreed, https://saulo-chaves.github.io/ProbBreed_site/, https://saulo-chaves.github.io/ProbBreed/
BugReports:	https://github.com/saulo-chaves/ProbBreed/issues
License:	AGPL (≥ 3)
Depends:	R (≥ 3.5.0)
Imports:	ggplot2, lifecycle, methods, Rcpp (≥ 0.12.0), RcppParallel (≥ 5.0.1), rlang, rstan (≥ 2.32.0), rstantools (≥ 2.4.0), stats, utils
LinkingTo:	StanHeaders (≥ 2.32.0), rstan (≥ 2.32.0), BH (≥ 1.72.0-2), Rcpp (≥ 0.12.0), RcppEigen (≥ 0.3.3.3.0), RcppParallel (≥ 5.0.1)
Suggests:	knitr, rmarkdown
Encoding:	UTF-8
UseLTO:	true
NeedsCompilation:	yes
RoxygenNote:	7.3.2
LazyData:	true
Biarch:	true
SystemRequirements:	GNU make
Packaged:	2025-03-17 21:18:52 UTC; Agroe
Author:	Saulo Chaves [aut, cre], Kaio Dias [aut, cph], Matheus Krause [aut]
Maintainer:	Saulo Chaves <saulo.chaves@ufv.br>
Repository:	CRAN
Date/Publication:	2025-03-17 22:10:04 UTC

The 'ProbBreed' package.

Description

ProbBreed uses probability theory under the Bayesian framework for calculating the risk of selecting candidates in a multi-environment context. Contained are functions used to fit a Bayesian multi-environment model (based on the available presets), extract posterior values and maximum posterior values, compute the variance components, check the model’s convergence, and calculate the probabilities. For both across and within-environments scopes, the package computes the probability of superior performance and the pairwise probability of superior performance. Furthermore, the probability of superior stability and the pairwise probability of superior stability across environments is estimated.

Author(s)

Maintainer: Saulo Chaves saulo.chaves@ufv.br (ORCID)

Authors:

• Kaio Dias kaio.o.dias@ufv.br (ORCID) [copyright holder]

• Matheus Krause mdkrause@iastate.edu (ORCID)

References

Stan Development Team (NA). RStan: the R interface to Stan. R package version 2.32.6. https://mc-stan.org

Dias, K. O. G, Santos J. P. R., Krause, M. D., Piepho H. -P., Guimarães, L. J. M., Pastina, M. M., and Garcia, A. A. F. (2022). Leveraging probability concepts for cultivar recommendation in multi-environment trials. Theoretical and Applied Genetics, 133(2):443-455. doi:10.1007/s00122-022-04041-y

See Also

Useful links:

• https://github.com/saulo-chaves/ProbBreed

• https://saulo-chaves.github.io/ProbBreed_site/

• https://saulo-chaves.github.io/ProbBreed/

• Report bugs at https://github.com/saulo-chaves/ProbBreed/issues

Bayesian model for multi-environment trials

Description

Fits a Bayesian multi-environment model using rstan, the R interface to Stan.

Usage

bayes_met(
  data,
  gen,
  loc,
  repl,
  trait,
  reg = NULL,
  year = NULL,
  res.het = FALSE,
  iter = 2000,
  cores = 1,
  chains = 4,
  pars = NA,
  warmup = floor(iter/2),
  thin = 1,
  seed = sample.int(.Machine$integer.max, 1),
  init = "random",
  verbose = FALSE,
  algorithm = c("NUTS", "HMC", "Fixed_param"),
  control = NULL,
  include = TRUE,
  show_messages = TRUE,
  ...
)

Arguments

Details

The function has nine available models, which will be fitted according to the options set in the arguments:

1. Entry-mean model : fitted when repl = NULL, reg = NULL and year = NULL:

   y = \mu + g + l + \varepsilon

   Where y is the phenotype, \mu is the intercept, g is the genotypic effect, l is the location (or environment) effect, and \varepsilon is the error (which contains the genotype-by-location interaction, in this case).

2. Randomized complete blocks design : fitted when repl is a single string. It will fit different models depending if reg and year are NULL:

   • reg = NULL and year = NULL :

     y = \mu + g + l + gl + r + \varepsilon

     where gl is the genotype-by-location effect, and r is the replicate effect.

   • reg = "reg" and year = NULL :

     y = \mu + g + m + l + gl + gm + r + \varepsilon

     where m is the region effect, and gm is the genotype-by-region effect.

   • reg = NULL and year = "year" :

     y = \mu + g + t + l + gl + gt + r + \varepsilon

     where t is the year effect, and gt is the genotype-by-year effect.

   • reg = "reg" and year = "year" :

     y = \mu + g + m + t + l + gl + gm + gt + r + \varepsilon

3. Incomplete blocks design : fitted when repl is a string vector of size 2. It will fit different models depending if reg and year are NULL:

   • reg = NULL and year = NULL :

     y = \mu + g + l + gl + r + b + \varepsilon

     where b is the block within replicates effect.

   • reg = "reg" and year = NULL :

     y = \mu + g + m + l + gl + gm + r + b + \varepsilon

   • reg = NULL and year = "year" :

     y = \mu + g + t + l + gl + gt + r + b + \varepsilon

   • reg = "reg" and year = "year" :

     y = \mu + g + m + t + l + gl + gm + gt + r + b + \varepsilon

The models described above have predefined priors:

x \sim \mathcal{N} \left( 0, S^{[x]} \right)

\sigma \sim \mathcal{HalfCauchy}\left( 0, S^{[\sigma]} \right)

where x can be any effect but the error, and \sigma is the standard deviation of the likelihood. If res.het = TRUE, then \sigma_k \sim \mathcal{HalfCauchy}\left( 0, S^{\left[ \sigma_k \right]} \right). The hyperpriors are set as follows:

S^{[x]} \sim \mathcal{HalfCauchy}\left( 0, \phi \right)

where \phi is the known global hyperparameter defined such as \phi = max(y) \times 10.

More details about the usage of bayes_met and other functions of the ProbBreed package can be found at https://saulo-chaves.github.io/ProbBreed_site/. Solutions to convergence or mixing issues can be found at https://mc-stan.org/misc/warnings.html.

Value

An object of S4 class stanfit representing the fitted results. Slot mode for this object indicates if the sampling is done or not.

Methods

sampling

signature(object = "stanmodel") Call a sampler (NUTS, HMC, or Fixed_param depending on parameters) to draw samples from the model defined by S4 class stanmodel given the data, initial values, etc.

See Also

rstan::sampling, rstan::stan, rstan::stanfit

Examples

mod = bayes_met(data = maize,
                gen = "Hybrid",
                loc = "Location",
                repl = c("Rep","Block"),
                trait = "GY",
                reg = "Region",
                year = NULL,
                res.het = TRUE,
                iter = 2000, cores = 2, chain = 4)

Extract outputs from stanfit objects obtained from bayes_met

Description

Extracts outputs of the Bayesian model fitted using bayes_met(), and provides some diagnostics.

Usage

extr_outs(model, probs = c(0.025, 0.975), verbose = FALSE)

Arguments

Details

More details about the usage of extr_outs and other functions of the ProbBreed package can be found at https://saulo-chaves.github.io/ProbBreed_site/.

Value

The function returns an object of class extr, which is a list with:

• variances : a data frame containing the variance components of the model effects, their standard deviation, naive standard error and highest posterior density interval.

• post : a list with the posterior of the effects, and the data generated by the model.

• map : a list with the maximum posterior values of each effect

• ppcheck : a matrix containing the p-values of maximum, minimum, median, mean and standard deviation; effective number of parameters, WAIC2 value, Rhat and effective sample size.

See Also

rstan::stan_diag, ggplot2::ggplot, rstan::check_hmc_diagnostics, plot.extr

Examples

mod = bayes_met(data = maize,
                gen = "Hybrid",
                loc = "Location",
                repl = c("Rep","Block"),
                trait = "GY",
                reg = "Region",
                year = NULL,
                res.het = TRUE,
                iter = 2000, cores = 2, chain = 4)

outs = extr_outs(model = mod,
                 probs = c(0.05, 0.95),
                 verbose = TRUE)

Maize real data set

Description

This dataset belongs to value of cultivation and use maize trials of Embrapa Maize and Sorghum, and was used by Dias et al. (2022). It contains the grain yield of 32 single-cross hybrids and four commercial checks (36 genotypes in total) evaluated in 16 locations across five regions or mega-environments. These trials were laid out in incomplete blocks design, using a block size of 6 and two replications per trial.

Usage

maize

Format

maize

A data frame with 823 rows and 6 columns:

Location

16 locations

Region

5 regions

Rep

2 replicates

Block

6 blocks

Hybrid

36 genotypes

GY

Grain yield (phenotypes)

Source

Dias, K. O. G, Santos J. P. R., Krause, M. D., Piepho H. -P., Guimarães, L. J. M., Pastina, M. M., and Garcia, A. A. F. (2022). Leveraging probability concepts for cultivar recommendation in multi-environment trials. Theoretical and Applied Genetics, 133(2):443-455. doi:10.1007/s00122-022-04041-y

Plots for the extr object

Description

Build plots using the outputs stored in the extr object.

Usage

## S3 method for class 'extr'
plot(x, ..., category = "ppdensity")

Arguments

Details

The available options are:

• ppdensity : Density plots of the empirical and sampled data, useful to assess the model's convergence.

• density : Density plots of the model's effects.

• histogram : Histograms of the model's effects.

• traceplot: Trace plot showing the changes in the effects' values across iterations and chains.

See Also

extr_outs

Examples

mod = bayes_met(data = maize,
                gen = "Hybrid",
                loc = "Location",
                repl = c("Rep","Block"),
                trait = "GY",
                reg = "Region",
                year = NULL,
                res.het = TRUE,
                iter = 2000, cores = 2, chain = 4)

outs = extr_outs(model = mod,
                 probs = c(0.05, 0.95),
                 verbose = TRUE)
plot(outs, category = "ppdensity")
plot(outs, category = "density")
plot(outs, category = "histogram")
plot(outs, category = "traceplot")

Plots for the probsup object

Description

Build plots using the outputs stored in the probsup object.

Usage

## S3 method for class 'probsup'
plot(x, ..., category = "perfo", level = "across")

Arguments

Details

The available options are:

• hpd : a caterpillar plot representing the marginal genotypic value of each genotype, and their respective highest posterior density interval (95% represented by the thick line, and 97.5% represented by the thin line). Available only if level = "across".

• perfo : if level = "across", a lollipop plot illustrating the probabilities of superior performance. If level = "within", a heatmap with the probabilities of superior performance within environments. If a model with reg and/or year is fitted, multiple plots are produced.

• stabi: a lollipop plot with the probabilities of superior stability. If a model with reg and/or year is fitted, multiple plots are produced. Available only if level = "across". Unavailable if an entry-mean model was used in bayes_met.

• pair_perfo : if level = "across", a heatmap representing the pairwise probability of superior performance (the probability of genotypes at the x-axis being superior. to those on the y-axis). If level = "within", a list of heatmaps representing the pairwise probability of superior performance within environments. If a model with reg and/or year is fitted, multiple plots (and multiple lists) are produced. Should this option is set, it is mandatory to store the outputs in an object. (e.g., pl <- plot(obj, category = "pair_perfo", level = "within")) so they can be visualized one at a time. The option level = "within" is unavailable if an entry-mean model was used in bayes_met.

• pair_stabi: a heatmap with the pairwise probabilities of superior stability (the probability of genotypes at the x-axis being more stable than those on the y-axis). If a model with reg and/or year is fitted, multiple plots are produced. Available only if level = "across". Unavailable if an entry-mean model was used in bayes_met.

• joint: a lollipop plot with the joint probabilities of superior performance and stability. Unavailable if an entry-mean model was used in bayes_met.

See Also

prob_sup

Examples

mod = bayes_met(data = maize,
                gen = "Hybrid",
                loc = "Location",
                repl = c("Rep","Block"),
                trait = "GY",
                reg = "Region",
                year = NULL,
                res.het = TRUE,
                iter = 2000, cores = 2, chain = 4)

outs = extr_outs(model = mod,
                 probs = c(0.05, 0.95),
                 verbose = TRUE)

results = prob_sup(extr = outs,
                   int = .2,
                   increase = TRUE,
                   save.df = FALSE,
                   verbose = FALSE)

plot(results, category = "hpd")
plot(results, category = "perfo", level = "across")
plot(results, category = "perfo", level = "within")
plot(results, category = "stabi")
plot(results, category = "pair_perfo", level = "across")
plwithin = plot(results, category = "pair_perfo", level = "within")
plot(results, category = "pair_stabi")
plot(results, category = "joint")

Print an object of class extr

Description

Print a extr object in R console

Usage

## S3 method for class 'extr'
print(x, ...)

Arguments

See Also

extr_outs

Print an object of class probsup

Description

Print a probsup object in R console

Usage

## S3 method for class 'probsup'
print(x, ...)

Arguments

See Also

prob_sup

Probabilities of superior performance and stability

Description

This function estimates the probabilities of superior performance and stability across environments, and probabilities of superior performance within environments.

Usage

prob_sup(extr, int, increase = TRUE, save.df = FALSE, verbose = FALSE)

Arguments

Details

Probabilities provide the risk of recommending a selection candidate for a target population of environments or for a specific environment. prob_sup computes the probabilities of superior performance and the probabilities of superior stability:

• Probability of superior performance

Let \Omega represent the subset of selected genotypes based on their performance across environments. A given genotype j will belong to \Omega if its genotypic marginal value (\hat{g}_j) is high or low enough compared to its peers. prob_sup leverages the Monte Carlo discretized sampling from the posterior distribution to emulate the occurrence of S trials. Then, the probability of the j^{th} genotype belonging to \Omega is the ratio of success (\hat{g}_j \in \Omega) events and the total number of sampled events, as follows:

Pr\left(\hat{g}_j \in \Omega \vert y \right) = \frac{1}{S}\sum_{s=1}^S{I\left(\hat{g}_j^{(s)} \in \Omega \vert y\right)}

where S is the total number of samples \left(s = 1, 2, ..., S \right), and I\left(g_j^{(s)} \in \Omega \vert y\right) is an indicator variable that can assume two values: (1) if \hat{g}_j^{(s)} \in \Omega in the s^{th} sample, and (0) otherwise. S is conditioned to the number of iterations and chains previously set at bayes_met.

Similarly, the within-environment probability of superior performance can be applied to individual environments. Let \Omega_k represent the subset of superior genotypes in the k^{th} environment, so that the probability of the j^{th} \in \Omega_k can calculated as follows:

Pr\left(\hat{g}_{jk} \in \Omega_k \vert y\right) = \frac{1}{S} \sum_{s=1}^S I\left(\hat{g}_{jk}^{(s)} \in \Omega_k \vert y\right)

where I\left(\hat{g}_{jk}^{(s)} \in \Omega_k \vert y\right) is an indicator variable mapping success (1) if \hat{g}_{jk}^{(s)} exists in \Omega_k, and failure (0) otherwise, and \hat{g}_{jk}^{(s)} = \hat{g}_j^{(s)} + \widehat{ge}_{jk}^{(s)}. Note that when computing within-environment probabilities, we are accounting for the interaction of the j^{th} genotype with the k^{th} environment.

The pairwise probabilities of superior performance can also be calculated across or within environments. This metric assesses the probability of the j^{th} genotype being superior to another experimental genotype or a commercial check. The calculations are as follows, across and within environments, respectively:

Pr\left(\hat{g}_{j} > \hat{g}_{j^\prime} \vert y\right) = \frac{1}{S} \sum_{s=1}^S I\left(\hat{g}_{j}^{(s)} > \hat{g}_{j^\prime}^{(s)} \vert y\right)

or

Pr\left(\hat{g}_{jk} > \hat{g}_{j^\prime k} \vert y\right) = \frac{1}{S} \sum_{s=1}^S I\left(\hat{g}_{jk}^{(s)} > \hat{g}_{j^\prime k}^{(s)} \vert y\right)

These equations are set for when the selection direction is positive. If increase = FALSE, > is simply switched by <.

• Probability of superior stability

This probability makes a direct analogy with the method of Shukla (1972): a stable genotype is the one that has a low genotype-by-environment interaction variance [var(\widehat{ge})]. Using the same probability principles previously described, the probability of superior stability is given as follows:

Pr \left[var \left(\widehat{ge}_{jk}\right) \in \Omega \vert y \right] = \frac{1}{S} \sum_{s=1}^S I\left[var \left(\widehat{ge}_{jk}^{(s)} \right) \in \Omega \vert y \right]

where I\left[var \left(\widehat{ge}_{jk}^{(s)} \right) \in \Omega \vert y \right] indicates if var\left(\widehat{ge}_{jk}^{(s)}\right) exists in \Omega (1) or not (0). Pairwise probabilities of superior stability are also possible in this context:

Pr \left[var \left(\widehat{ge}_{jk} \right) < var\left(\widehat{ge}_{j^\prime k} \right) \vert y \right] = \frac{1}{S} \sum_{s=1}^S I \left[var \left(\widehat{ge}_{jk} \right)^{(s)} < var \left(\widehat{ge}_{j^\prime k} \right)^{(s)} \vert y \right]

Note that j will be superior to j^\prime if it has a lower variance of the genotype-by-environment interaction effect. This is true regardless if increase is set to TRUE or FALSE.

The joint probability independent events is the product of the individual probabilities. The estimated genotypic main effects and the variances of GEI effects are independent by design, thus the joint probability of superior performance and stability as follows:

Pr \left[\hat{g}_j \in \Omega \cap var \left(\widehat{ge}_{jk} \right) \in \Omega \right] = Pr \left(\hat{g}_j \in \Omega \right) \times Pr \left[var \left(\widehat{ge}_{jk} \right) \in \Omega \right]

The estimation of these probabilities are strictly related to some key questions that constantly arises in plant breeding:

• What is the risk of recommending a selection candidate for a target population of environments?

• What is the probability of a given selection candidate having good performance if recommended to a target population of environments? And for a specific environment?

• What is the probability of a given selection candidate having better performance than a cultivar check in the target population of environments? And in specific environments?

• How probable is it that a given selection candidate performs similarly across environments?

• What are the chances that a given selection candidate is more stable than a cultivar check in the target population of environments?

• What is the probability that a given selection candidate having a superior and invariable performance across environments?

More details about the usage of prob_sup, as well as the other function of the ProbBreed package can be found at https://saulo-chaves.github.io/ProbBreed_site/.

Value

The function returns an object of class probsup, which contains two lists, one with the across-environments probabilities, and another with the within-environments probabilities. If an entry-mean model was used in ProbBreed::bayes_met, only the across list will be available.

The across list has the following elements:

• g_hpd: Highest posterior density (HPD) of the posterior genotypic main effects.

• perfo: the probabilities of superior performance.

• pair_perfo: the pairwise probabilities of superior performance.

• stabi: a list with the probabilities of superior stability. It contains the data frames gl, gm (when reg is not NULL) and gt (when year is not NULL). Unavailable if an entry-mean model was used in bayes_met.

• pair_stabi: a list with the pairwise probabilities of superior stability. It contains the data frames gl, gm (when reg is not NULL) and gt (when year is not NULL). Unavailable if an entry-mean model was used in bayes_met.

• joint_prob: the joint probabilities of superior performance and stability. Unavailable if an entry-mean model was used in bayes_met.

The within list has the following elements:

• perfo: a list of data frames containing the probabilities of superior performance within locations (gl), regions (gm) and years (gt).

• pair_perfo: lists with the pairwise probabilities of superior performance within locations (gl), regions (gm) and years (gt).

References

Dias, K. O. G, Santos J. P. R., Krause, M. D., Piepho H. -P., Guimarães, L. J. M., Pastina, M. M., and Garcia, A. A. F. (2022). Leveraging probability concepts for cultivar recommendation in multi-environment trials. Theoretical and Applied Genetics, 133(2):443-455. doi:10.1007/s00122-022-04041-y

Shukla, G. K. (1972) Some statistical aspects of partioning genotype environmental componentes of variability. Heredity, 29:237-245. doi:10.1038/hdy.1972.87

See Also

plot.probsup

Examples

mod = bayes_met(data = maize,
                gen = "Hybrid",
                loc = "Location",
                repl = c("Rep","Block"),
                trait = "GY",
                reg = "Region",
                year = NULL,
                res.het = TRUE,
                iter = 2000, cores = 2, chain = 4)

outs = extr_outs(model = mod,
                 probs = c(0.05, 0.95),
                 verbose = TRUE)

results = prob_sup(extr = outs,
                   int = .2,
                   increase = TRUE,
                   save.df = FALSE,
                   verbose = FALSE)

Soybean real data set

Description

This dataset belongs to the USDA Northern Region Uniform Soybean Tests, and it is a subset of the data used by Krause et al. (2023). It contains the empirical best linear unbiased estimates of genotypic means of the seed yield from 39 experimental genotypes evaluated in 14 locations. The original data, available at the package SoyURT, has 4,257 experimental genotypes evaluated at 63 locations and 31 years resulting in 591 location-year combinations (environments) with 39,006 yield values.

Usage

soy

Format

soy

A data frame with 823 rows and 3 columns:

Loc

14 locations

Gen

39 experimental genotypes

Y

435 EBLUEs (phenotypes)

Source

Krause MD, Dias KOG, Singh AK, Beavis WD. 2023. Using soybean historical field trial data to study genotype by environment variation and identify mega-environments with the integration of genetic and non-genetic factors. bioRxiv : the preprint server for biology. doi: https://doi.org/10.1101/2022.04.11.487885