Maeen Network Mirror :: Help for package bayesammi

Type:

Package

Title:

Bayesian Estimation of the Additive Main Effects and Multiplicative Interaction Model

Version:

1.0.0

Date:

2026-11-04

Description:

Bayesian estimation of the additive main effects and multiplicative interaction (AMMI) model of multi-environment trials by the Gibbs sampler of Crossa, Perez-Elizalde, Jarquin, Cotes, Viele, Liu and Cornelius (2011) <doi:10.2135/cropsci2010.06.0343>: posterior distributions of the overall mean, genotype and environment effects, singular values, genotype and environment scores (von Mises-Fisher conditionals on the constrained spheres) and error variance; several chains with burn-in and thinning; posterior summaries with highest posterior density intervals; the convergence diagnostics of Raftery-Lewis, Heidelberger-Welch and Gelman-Rubin; the AMMI biplot with highest posterior density regions for inference on genotypes and environments; least-squares AMMI; and a plain-language interpretation of every result, computed from its numbers. Beyond the article: unbalanced trials with missing plots and cells (data augmentation), environment-specific error variances, the Bayesian GGE (sites regression) model of Perez-Elizalde, Jarquin and Crossa (2012) <doi:10.1007/s13253-011-0063-9>, probabilities of superior performance and stability (Dias et al. 2022 <doi:10.1007/s00122-022-04041-y>), posterior ranks, Bayesian stability measures, predictions and the winner of every environment, model choice by DIC and WAIC with posterior predictive checks, draws for the 'posterior', 'bayesplot' and 'coda' packages, a compiled sampler, and a built-in 'shiny' application. Bayesian stability analysis: the Finlay-Wilkinson, Eberhart-Russell and Shukla models by Gibbs sampling, and the posterior distributions of classical stability indices (Wricke, Shukla, Lin-Binns, Annicchiarico, Eskridge's safety-first) and of AMMI-based indices (ASV, SIPC, EV, Za, WAAS, WAASY). Every result prints the references of its method.

License:

GPL-2

URL:

https://myaseen208.com/bayesammi/, https://github.com/myaseen208/bayesammi

BugReports:

https://github.com/myaseen208/bayesammi/issues

Encoding:

UTF-8

Language:

en-US

LazyData:

true

Depends:

R (≥ 4.2)

Imports:

coda, collapse, data.table, ggplot2, ggrepel, grDevices, MASS, Rcpp, scales, stats, utils

LinkingTo:

Rcpp, RcppArmadillo

Suggests:

bayesplot, bslib, DT, flextable, knitr, lme4, loo, magick, officer, pkgload, posterior, quarto, ragg, R.rsp, shiny, svglite, testthat (≥ 3.0.0), withr

VignetteBuilder:

quarto, R.rsp

Config/testthat/edition:

Config/skill/prefix:

bam

RoxygenNote:

7.3.1

NeedsCompilation:

yes

Packaged:

2026-11-10 02:44:51 UTC; myaseen208

Author:

Muhammad Yaseen

[aut, cre], Jose Crossa [aut, ctb], Sergio Perez-Elizalde [aut, ctb], Diego Jarquin [aut, ctb], Jose Miguel Cotes [aut, ctb], Kert Viele [aut, ctb], Genzhou Liu [aut, ctb], Paul L. Cornelius [aut, ctb], Julian Garcia Abadillo Velasco [aut, ctb]

Maintainer:

Muhammad Yaseen <myaseen208@gmail.com>

Repository:

CRAN

Date/Publication:

2026-11-10 09:30:23 UTC

bayesammi: Bayesian Estimation of the Additive Main Effects and Multiplicative Interaction Model

Description

bayesammi estimates the additive main effects and multiplicative interaction (AMMI) model of multi-environment trials in plant breeding by Bayesian methods, with the Gibbs sampler of Crossa et al. (2011). Every parameter has a posterior distribution, so the AMMI biplot gets statistical inference: highest posterior density (HPD) regions show which genotypes and environments interact significantly and which form homogeneous groups. Beyond the article, the package fits unbalanced trials, environment-specific error variances and the GGE (sites regression) model, turns the posterior into decisions (winners, probabilities of superior performance and stability, ranks, stability measures, predictions), compares and checks models, and runs every analysis in a built-in application. It also fits the Finlay-Wilkinson, Eberhart-Russell and Shukla stability models by Gibbs sampling and gives posterior distributions of classical and AMMI-based stability indices. The package is for breeders and agronomists who analyse genotype-by-environment interaction and want uncertainty statements about their conclusions.

Details

Functions, in the order of an analysis:

ge_mean, ge_var, g_eff, e_eff, ge_eff: cell means and variances, main effects, interaction effects and their singular value decomposition.
ge_ammi, ge_model: least-squares AMMI and the mixed model of the error variance.
bam_priors, bayes_ammi: priors and the Gibbs sampler of the AMMI or GGE model, with a common or environment-specific error variances and missing plots or cells (chains, burn-in, thinning, posterior summaries with HPD intervals, convergence diagnostics, DIC and WAIC).
biplots: the biplot with HPD regions.
bam_predict, bam_winners, bam_superiority, bam_ranks, bam_stability: decisions for breeders.
bam_compare, bam_ppc, bam-draws: model choice, posterior predictive checks, and the draws for posterior, bayesplot and coda.
bam_stab_priors, bayes_stability: Bayesian Finlay-Wilkinson, Eberhart-Russell and Shukla stability models; bam_indices, bam_ammi_indices: posterior distributions of classical and AMMI-based stability indices.
bam_cite: the references of the methods behind a result (also printed under every result).
bam_app: the built-in application.
bam_interpret: a plain-language interpretation of a result.
matrix_k, orthnorm: orthonormal contrasts and Gram-Schmidt orthonormalisation.

Every analysis returns the named list documented on its help page, with a class and methods (print, summary, plot, as.data.table); see bam-methods. Utilities: bam_theme, bam_save, bam_notation, bam_reference.

Data. Maiz: grain yield of 9 maize genotypes in 20 sites with 4 replicates (the example of Crossa et al. 2011).

Vignettes. Start with vignette("bayesammi_01_getting_started", package = "bayesammi"); browseVignettes("bayesammi") lists all of them, each as HTML and as PDF.

Notation

Scalars are upper-case letters in normal face, indices lower-case italic, vectors lower-case bold and matrices upper-case bold; Greek letters follow the same rule. The full table is bam_notation().

`G`, `E`, `B`, `K`	numbers of genotypes, environments, replicates, bilinear terms
`\mathbf{Y}`	`G \times E` table of cell means
`\mu`, `\boldsymbol{\tau}`, `\boldsymbol{\delta}`	overall mean, genotype and environment effects
`\lambda_k`, `\boldsymbol{\alpha}_k`, `\boldsymbol{\gamma}_k`	singular value and genotype and environment score vectors of term `k`
`\sigma^2`, `\sigma^2_j`	common and environment-specific error variance
`\eta_{ij}`	expected value of genotype `i` in environment `j`
`g_i`, `h_j`, `1 + b_i`	genotype effect, environmental index and slope of the stability models
`d_{ij}`, `\sigma^2_{d,i}`	deviation from the regression and its (genotype-specific) variance

Author(s)

Muhammad Yaseen (myaseen208@gmail.com), ORCID 0000-0002-5923-1714.

References

Crossa, J., Perez-Elizalde, S., Jarquin, D., Cotes, J. M., Viele, K., Liu, G. and Cornelius, P. L. (2011). Bayesian estimation of the additive main effects and multiplicative interaction model. Crop Science, 51(4), 1458–1469. doi:10.2135/cropsci2010.06.0343

Perez-Elizalde, S., Jarquin, D. and Crossa, J. (2012). A general Bayesian estimation method of linear-bilinear models applied to plant breeding trials with genotype x environment interaction. Journal of Agricultural, Biological, and Environmental Statistics, 17(1), 15–37. doi:10.1007/s13253-011-0063-9

Viele, K. and Srinivasan, C. (2000). Parsimonious estimation of multiplicative interaction in analysis of variance using Kullback-Leibler information. Journal of Statistical Planning and Inference, 84(1–2), 201–219.

Maize multi-environment trial

Description

Grain yield of 9 maize (Zea mays L.) genotypes in 20 sites, a randomized complete block design with 4 replicates in every site: the example of Crossa et al. (2011), also used by Burgueno et al. (2008).

Usage

Maiz

Format

A data.frame with 720 rows and 7 columns:

site: site (environment), 1 to 20.
entry: genotype, 1 to 9.
fold10, fold5, fold3: fold numbers of 10-, 5- and 3-fold cross-validation (by genotype).
rep: replicate within the site, 1 to 4.
y: grain yield (kg/ha).

Details

The grand mean is 4858 kg/ha; the least-squares singular values of the interaction are 5923 and 3070 for the first two terms (Crossa et al. 2011).

Change from the earlier versions. The help page described another data set (1320 rows; columns Gen, Institute, Rep, Block, Env, Yield).

Author(s)

Muhammad Yaseen (myaseen208@gmail.com), Jose Crossa (j.crossa@cgiar.org)

Source

International Maize and Wheat Improvement Center (CIMMYT).

References

Burgueno, J., Crossa, J., Cornelius, P. L. and Yang, R.-C. (2008). Using factor analytic models for joining environments and genotypes without crossover genotype x environment interaction. Crop Science, 48(4), 1291–1305.

Examples

data(Maiz)
str(Maiz)
ge_mean(
    .data = Maiz
  , .y    = y
  , .gen  = entry
  , .env  = site
  )$grand_mean

The draws of a fit for posterior, bayesplot and coda

Description

The draws of bayes_ammi() are a coda::mcmc.list (component draws). These methods hand them to the packages posterior and coda; every plot of bayesplot (mcmc_trace(), mcmc_areas(), mcmc_rank_overlay(), ...) then works on them, as does posterior::summarise_draws() (rank-normalised \hat{R}, bulk and tail effective sample sizes).

Usage

## S3 method for class 'bam_bayes_ammi'
as_draws_df(x, ...)

## S3 method for class 'bam_bayes_ammi'
as_draws(x, ...)

## S3 method for class 'bam_bayes_ammi'
as.mcmc.list(x, ...)

Arguments

x

a result of bayes_ammi().

...

unused.

Value

as_draws_df(): a posterior::draws_df; as_draws(): a posterior::draws_list; as.mcmc.list(): the coda::mcmc.list.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Gabry, J., Simpson, D., Vehtari, A., Betancourt, M. and Gelman, A. (2019). Visualization in Bayesian workflow. Journal of the Royal Statistical Society A, 182(2), 389–402. doi:10.1111/rssa.12378

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. and Burkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved R-hat for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718. doi:10.1214/20-BA1221

Examples

data(Maiz)
fit <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 500
    , .k      = 2
    , .burnin = 100
    , .thin   = 1
    , .chains = 2
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )
if (requireNamespace("posterior", quietly = TRUE)) {
  dr <- posterior::as_draws_df(x = fit)
  posterior::summarise_draws(
      .x = posterior::subset_draws(
          x        = dr
        , variable = c("mu", "lambda[1]", "lambda[2]")
        )
    )
}
if (requireNamespace("bayesplot", quietly = TRUE)) {
  bayesplot::mcmc_trace(x = fit$draws, pars = c("lambda[1]", "lambda[2]"))
}

Methods for bayesammi results

Description

The analysis functions of bayesammi return the named list they have always returned, now with the class c("bam_<analysis>", "bam_result", "list"). Components are read as before (x$ANOVA, x[["Means"]]); the methods below print, summarise and plot the result and turn any component into a data.table.

Usage

## S3 method for class 'bam_result'
print(x, digits = 4L, ...)

## S3 method for class 'bam_result'
summary(object, ...)

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

## S3 method for class 'bam_result'
format(x, ...)

## S3 method for class 'bam_result'
as.data.table(x, keep.rownames = FALSE, ..., component = 1L)

## S3 method for class 'bam_result'
as.data.frame(x, row.names = NULL, optional = FALSE, ..., component = 1L)

## S3 method for class 'bam_result'
plot(x, ...)

Arguments

x, object

a bayesammi result.

digits

significant digits of printed numbers (default 4).

...

passed on (to print() of the components, or to the plotting function for plot()).

keep.rownames, row.names, optional

as in data.table::as.data.table() and as.data.frame(); row names of a matrix or data frame component are always kept, in a first column term.

component

name or number of the component to convert (default: the first one).

Value

print() returns x invisibly; summary() an object of class summary.bam_result; format() a character string; as.data.table() and as.data.frame() the chosen component as a table; plot() a ggplot object.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

Posterior distributions of AMMI-based stability indices

Description

bam_ammi_indices() computes the AMMI stability value and the other stability indices built from the singular values and genotype scores of an AMMI or GGE fit on every posterior draw: ASV (Purchase et al. 2000), SIPC (Sneller et al. 1997), EV (Zobel 1994), Za (Zali et al. 2012), WAAS and the mean-stability index WAASY (Olivoto et al. 2019).

Usage

bam_ammi_indices(fit, higher_better = TRUE, .weight_y = 50)

Arguments

fit

a result of bayes_ammi().

higher_better

TRUE (default) when larger values of the trait are better (for WAASY).

.weight_y

weight \theta_Y (0-100) of the mean in WAASY (default 50, equal weights).

Details

With K multiplicative terms, singular values \lambda_k, genotype scores \alpha_{ik}, coordinates \mathrm{PC}_{ik} = \lambda_k^{1/2}\alpha_{ik} and the share of term k in the fitted interaction \theta_k = \lambda_k^2/\sum_m\lambda_m^2, on every draw:

ASV: \sqrt{(\lambda_1^2/\lambda_2^2\,\mathrm{PC}_{i1})^2 + \mathrm{PC}_{i2}^2}, with the ratio of the posterior means of the singular values (as in bam_stability(); a ratio of draws explodes when \lambda_2 comes near zero); needs K \ge 2.
SIPC: \sum_k|\mathrm{PC}_{ik}|.
EV: \sum_k\alpha_{ik}^2/K.
Za: \sum_k\theta_k|\alpha_{ik}|.
WAAS: \sum_k\theta_k|\mathrm{PC}_{ik}| (the weights sum to one).
WAASY: (\theta_Y R^Y_i + (100 - \theta_Y)R^W_i)/100, where R^Y_i rescales the genotype means to 0-100 (100 for the best) and R^W_i rescales WAAS to 0-100 (100 for the most stable); larger is better.

All but WAASY are smaller for more stable genotypes. The shares \theta_k are those of the K fitted terms; Zali et al. (2012) and Olivoto et al. (2019) use the shares of the significant terms in the whole interaction, which differ by a constant factor for Za and not at all for WAAS.

Value

An object of class c("bam_ammi_indices", "bam_result", "list") with indices and means, as in bam_indices(). Methods: print(), plot(), bam_interpret(), bam_cite().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Purchase, J. L., Hatting, H. and van Deventer, C. S. (2000). Genotype x environment interaction of winter wheat (Triticum aestivum L.) in South Africa: II. Stability analysis of yield performance. South African Journal of Plant and Soil, 17(3), 101–107.

Sneller, C. H., Kilgore-Norquest, L. and Dombek, D. (1997). Repeatability of yield stability statistics in soybean. Crop Science, 37(2), 383–390.

Zobel, R. W. (1994). Stress resistance and root systems. In Proceedings of the Workshop on Adaptation of Plants to Soil Stress, INTSORMIL Publication 94-2, 80–99. University of Nebraska-Lincoln.

Zali, H., Farshadfar, E., Sabaghpour, S. H. and Karimizadeh, R. (2012). Evaluation of genotype x environment interaction in chickpea using measures of stability from AMMI model. Annals of Biological Research, 3(7), 3126–3136.

Olivoto, T., Lucio, A. D. C., da Silva, J. A. G., Marchioro, V. S., de Souza, V. Q. and Jost, E. (2019). Mean performance and stability in multi-environment trials I: combining features of AMMI and BLUP techniques. Agronomy Journal, 111(6), 2949–2960. doi:10.2134/agronj2019.03.0220

Examples

data(Maiz)
fit <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 600
    , .k      = 2
    , .burnin = 100
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )
ai <-
  bam_ammi_indices(
      fit           = fit
    , higher_better = TRUE
    , .weight_y     = 50
    )
ai$means
plot(x = ai)

Bayesian AMMI and GGE analysis in a web browser

Description

bam_app() starts the application of bayesammi: load a trial (the example Maiz or a CSV file with one row per plot), fit the Bayesian AMMI or GGE model, check convergence, draw the biplot with HPD regions, and turn the posterior into decisions (which genotype wins where, probabilities of superior performance and stability, ranks, Bayesian stability, predictions of missing cells), and compare models. Every page shows an interpretation computed from the results (bam_interpret()), the output as R prints it, formatted tables, plots and the R code that repeats the analysis, each with a Copy button.

Usage

bam_app(launch.browser = interactive(), ...)

bam_app_object()

Arguments

launch.browser

passed to shiny::runApp(): TRUE (default in an interactive session) opens the application in the web browser.

...

further arguments to shiny::runApp(), such as port.

Details

The application needs the suggested packages shiny, bslib and DT. Plots are sent to the browser as WebP images when magick is installed (smaller than PNG; PNG otherwise) and can be downloaded as PNG, PDF or WebP. Text is set in IBM Plex Sans and output in IBM Plex Mono when the browser can reach Google Fonts, in the system fonts otherwise.

Value

bam_app() is called for its side effect and returns the value of shiny::runApp() when the application is closed; bam_app_object() returns the application object (for shiny::runApp(), shiny::testServer() or deployment on a server).

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

Examples

if (interactive()) {
  bam_app(
      launch.browser = TRUE
    )
}

References of a result or of the package

Description

bam_cite() gives the references of the methods behind a result of the package, or the references of every method in the package, as text or BibTeX. citation("bayesammi") cites the package and its main methods.

Usage

bam_cite(x = NULL, style = c("text", "bibtex"))

Arguments

x

a result of a function of the package, or NULL (default) for every reference of the package.

style

"text" (default) or "bibtex".

Value

An object of class bibentry (printed in the chosen style).

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

Examples

bam_cite(
    x     = NULL
  , style = "text"
  )
data(Maiz)
g <-
  ge_ammi(
      .data = Maiz
    , .y    = y
    , .gen  = entry
    , .env  = site
    , .rep  = rep
    )
bam_cite(
    x     = g
  , style = "bibtex"
  )

Compare Bayesian AMMI and GGE models

Description

bam_compare() compares fits of bayes_ammi() that differ in the number of multiplicative terms, the model (AMMI or GGE), the error variances or the priors, by the deviance information criterion (DIC) and the widely applicable information criterion (WAIC). Smaller values predict new data better.

Usage

bam_compare(...)

Arguments

...

fits of bayes_ammi() or bayes_stability() to the same data, named or not (names label the rows).

Details

With D(\boldsymbol{\theta}) = -2\log p(\mathbf{y} \mid \boldsymbol{\theta}) over the observed plots,

\mathrm{DIC} = D(\bar{\boldsymbol{\theta}}) + 2p_D, \qquad p_D = \overline{D(\boldsymbol{\theta})} - D(\bar{\boldsymbol{\theta}})

(Spiegelhalter et al. 2002), the deviance evaluated at the posterior means of the cell means and variances; and

\mathrm{WAIC} = -2\left(\sum_o \log\frac{1}{S}\sum_s p(y_o \mid \boldsymbol{\theta}^{(s)}) - \sum_o \mathrm{Var}_s\log p(y_o \mid \boldsymbol{\theta}^{(s)})\right)

(Watanabe 2010; Gelman et al. 2014), accumulated over all draws during sampling. All fits must use the same data.

The standard error of a difference in WAIC comes from the plot-by-plot differences of the contributions w_o to WAIC: \mathrm{se} = \sqrt{n\,\mathrm{Var}_o(w_o^{(m)} - w_o^{(\mathrm{best})})} (Vehtari et al. 2017). Because the models are fitted to the same plots, it is smaller than the standard error of either WAIC. A difference of less than about two of its standard errors does not separate two models.

Value

An object of class c("bam_compare", "bam_result", "list") with table (data.table sorted by WAIC: model, terms, variances, DIC, pD, WAIC, p_waic, se (standard error of WAIC), lambda_last_lower (lower HPD limit of the last singular value), lambda_last_mean, dWAIC (difference from the best) and dWAIC_se (its standard error)). Methods: print(), plot(), bam_interpret().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Spiegelhalter, D. J., Best, N. G., Carlin, B. P. and van der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society B, 64(4), 583–639. doi:10.1111/1467-9868.00353

Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11, 3571–3594.

Gelman, A., Hwang, J. and Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24(6), 997–1016. doi:10.1007/s11222-013-9416-2

Vehtari, A., Gelman, A. and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413–1432. doi:10.1007/s11222-016-9696-4

Examples

data(Maiz)
fits <-
  lapply(
      X   = 1:2
    , FUN = function(K) {
        bayes_ammi(
            .data   = Maiz
          , .y      = y
          , .gen    = entry
          , .env    = site
          , .rep    = rep
          , .nIter  = 600
          , .k      = K
          , .burnin = 100
          , .thin   = 1
          , .chains = 1
          , .priors = bam_priors(
                mu_mean     = 0
              , mu_var      = 1e+8
              , tau_mean    = 0
              , tau_var     = 1e+8
              , delta_mean  = 0
              , delta_var   = 1e+8
              , lambda_mean = 0
              , lambda_var  = 1e+8
              , nu          = 1
              , s2          = NULL
              )
          , .level  = 0.95
          , .seed   = 1
          , .align  = TRUE
          , .model  = "AMMI"
          , .het    = FALSE
          , .waic   = TRUE
          )
      }
    )
bam_compare(
    AMMI1 = fits[[1]]
  , AMMI2 = fits[[2]]
  )

Posterior distributions of classical stability indices

Description

bam_indices() computes five classical stability and superiority indices on every posterior draw of the cell means of a fit, so that each index has a posterior mean, an HPD interval and, for every genotype, a probability of being the best. The indices are those of Wricke (1962), Shukla (1972), Lin and Binns (1988), Annicchiarico (1992) and Eskridge (1990).

Usage

bam_indices(fit, higher_better = TRUE, .prob = 0.25, .alpha = 0.05)

Arguments

fit

a result of bayes_ammi() or bayes_stability().

higher_better

TRUE (default) when larger values of the trait are better.

.prob

error probability p of Annicchiarico's index (default 0.25).

.alpha

probability \alpha of the safety-first index (default 0.05).

Details

With \eta_{ij} the cell mean of genotype i in environment j (one posterior draw), \bar{\eta}_{i\cdot}, \bar{\eta}_{\cdot j} and \bar{\eta}_{\cdot\cdot} the genotype, environment and overall means, and the interaction \epsilon_{ij} = \eta_{ij} - \bar{\eta}_{i\cdot} - \bar{\eta}_{\cdot j} + \bar{\eta}_{\cdot\cdot}:

wricke: ecovalence W_i = \sum_j \epsilon_{ij}^2; smaller is more stable.
shukla: stability variance \sigma^2_i = \frac{G}{(G-2)(E-1)}W_i - \frac{\sum_k W_k}{(G-1)(G-2)(E-1)}; smaller is more stable.
lin_binns: superiority P_i = \sum_j(\eta_{ij} - M_j)^2/(2E), M_j the best cell mean of environment j; smaller is better.
annicchiarico: reliability I_i = \bar{Z}_i - z_{1-p}s_{Z_i} with Z_{ij} = 100\,\eta_{ij}/\bar{\eta}_{\cdot j}; larger is better (with higher_better = FALSE: \bar{Z}_i + z_{1-p}s_{Z_i}, smaller better).
safety_first: Eskridge's index with the variance across environments, \bar{\eta}_{i\cdot} - z_{1-\alpha}s_i; larger is better (with higher_better = FALSE: \bar{\eta}_{i\cdot} + z_{1-\alpha}s_i, smaller better). bayes_stability() gives the model-based versions of the regression and Shukla models.

Because the indices are computed from the posterior of the cell means, they measure the stability of the expected performance, free of the plot error; the HPD intervals show how well the trial determines them.

Value

An object of class c("bam_indices", "bam_result", "list") with

indices: data.table in long form: index, Gen, mean, lower, upper (HPD interval), prob_best (posterior probability that the genotype is the best on this index) and rank_mean (expected rank, 1 = best).
means: data.table: one row per genotype with the posterior means of the indices and the mean of their expected ranks (rank_overall), sorted by it.

Methods: print(), plot(), bam_interpret(), bam_cite().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Wricke, G. (1962). Uber eine Methode zur Erfassung der okologischen Streubreite in Feldversuchen. Zeitschrift fur Pflanzenzuchtung, 47, 92–96.

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

Lin, C. S. and Binns, M. R. (1988). A superiority measure of cultivar performance for cultivar x location data. Canadian Journal of Plant Science, 68(1), 193–198. doi:10.4141/cjps88-018

Annicchiarico, P. (1992). Cultivar adaptation and recommendation from alfalfa trials in Northern Italy. Journal of Genetics and Breeding, 46, 269–278.

Eskridge, K. M. (1990). Selection of stable cultivars using a safety-first rule. Crop Science, 30, 369–374.

Examples

data(Maiz)
fit <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 600
    , .k      = 2
    , .burnin = 100
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )
ind <-
  bam_indices(
      fit           = fit
    , higher_better = TRUE
    , .prob         = 0.25
    , .alpha        = 0.05
    )
ind$means
plot(x = ind)

Interpretation of a Bayesian AMMI analysis

Description

bam_interpret() turns a result of bayes_ammi(), biplots() or ge_ammi() into a plain-language interpretation: whether the chains have converged, which genotypes and environments differ from the average (main effects whose HPD interval excludes zero), how strong the bilinear terms are, which genotypes and environments interact significantly and in which direction, which are stable, and what to check before acting on the results. Every statement is computed from the numbers of the result.

Usage

bam_interpret(x, higher_better = TRUE, top = 5L, ...)

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

## S3 method for class 'bam_interpretation'
format(x, style = c("text", "markdown"), width = getOption("width", 80L), ...)

Arguments

x

a result of bayes_ammi(), biplots() or ge_ammi().

higher_better

TRUE (default) when larger values of the trait are better.

top

number of genotypes or environments named in lists (default 5).

...

unused.

style

"text" (default) or "markdown".

width

line width of the text (default getOption("width")).

Details

The rules:

convergence: potential scale reduction factor below 1.1 (two or more chains), effective sample size of at least 400, Raftery-Lewis dependence factor below 5 and the Heidelberger-Welch tests passed (Crossa et al. 2011);
main effects: a genotype or environment differs from the average when the HPD interval of its effect excludes zero;
bilinear terms: a term whose HPD interval of \lambda_k reaches down to near zero (lower limit below a tenth of the posterior mean) is weakly supported;
interaction: a genotype (environment) interacts significantly on the first term when the HPD interval of its first score excludes zero; genotypes and environments with significant scores of the same sign interact positively (the genotype does better there than its main effect predicts), of opposite sign negatively;
stability (biplots()): a genotype is stable when its HPD region contains the origin.

Value

An object of class "bam_interpretation": a list with analysis, title, headline, findings, recommendations and cautions (sentences), facts (the numbers and lists behind the sentences), guide (how the method is read) and sources (the full references of the method; bam_cite() gives them as BibTeX). print() shows the text; format(x, style = "markdown") returns it as Markdown.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Examples

data(Maiz)
fm1 <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 500
    , .k      = 2
    , .burnin = 100
    , .thin   = 1
    , .chains = 2
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = c("AMMI", "GGE")
    , .het    = FALSE
    , .waic   = TRUE
    )
txt <-
  bam_interpret(
      x             = fm1
    , higher_better = TRUE
    , top           = 5
    )
txt
cat(format(x = txt, style = "markdown"))

Notation of bayesammi

Description

The symbols used in the help pages and vignettes of bayesammi, with their meaning and the argument or result component that holds them. The notation follows one rule throughout: scalars are upper-case letters in normal face (P, R), indices lower-case italic (i, j), vectors lower-case bold (\mathbf{y}, \boldsymbol{\beta}) and matrices upper-case bold (\mathbf{X}, \boldsymbol{\Sigma}).

Usage

bam_notation(format = c("table", "path", "tex"))

Arguments

format

"table" (default: a data.table with columns symbol (LaTeX), text (plain text), kind (scalar, index, vector, matrix, set, function, operator), meaning and argument), "path" (the path of notation.csv) or "tex" (the path of Notations.tex, the notation appendix used in the articles).

Value

See format.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

Examples

bam_notation(
    format = c("table", "path", "tex")[1]
  )

Posterior predictive checks

Description

bam_ppc() compares statistics of the observed trial with their distribution in trials replicated from the posterior predictive distribution (Gelman et al. 1996). A posterior predictive p-value near 0 or 1 shows a feature of the data that the model does not reproduce, for instance unequal error variances or a missing interaction term.

Usage

bam_ppc(fit, .draws = 1000L, .seed = NULL)

Arguments

fit

a result of bayes_ammi() or bayes_stability().

.draws

number of posterior draws used (default 1000).

.seed

seed for the replicated trials (the random number state of the session is restored).

Details

For a subset of .draws posterior draws, a trial \mathbf{y}^{\mathrm{rep}} with the observed plots is simulated from N(\eta_{ij}, \sigma^2_j), and for each statistic T p = \Pr(T(\mathbf{y}^{\mathrm{rep}}) \ge T(\mathbf{y}) \mid \mathbf{y}). The statistics: variance of the cell means, largest and smallest cell mean, variance of the genotype means, first singular value of the interaction effects of the cell means, pooled within-cell variance, and the ratio of the largest to the smallest within-environment variance.

Value

An object of class c("bam_ppc", "bam_result", "list") with checks (data.table: statistic, observed, mean_rep, p_value) and replicated (data.table of the replicated statistics). Methods: print(), plot(), bam_interpret().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Gelman, A., Meng, X.-L. and Stern, H. (1996). Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6(4), 733–807.

Examples

data(Maiz)
fit <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 600
    , .k      = 2
    , .burnin = 100
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )
pc <-
  bam_ppc(
      fit    = fit
    , .draws = 300
    , .seed  = 1
    )
pc$checks
plot(x = pc)

Predictions of every genotype-environment cell

Description

bam_predict() gives the posterior distribution of the expected response of every genotype in every environment, including the cells that were not observed, and of the genotype means over the environments.

Usage

bam_predict(fit, .level = NULL)

Arguments

fit

a result of bayes_ammi() or bayes_stability() (bam_stability() needs bayes_ammi()).

.level

probability of the HPD intervals (default: that of the fit).

Details

For draw s,

\eta^{(s)}_{ij} = \mu^{(s)} + \tau^{(s)}_i + \delta^{(s)}_j + \sum_{k=1}^{K}\lambda^{(s)}_k \alpha^{(s)}_{ik}\gamma^{(s)}_{jk}

(without \tau_i in the GGE model); the genotype mean over the environments is \bar{\eta}_{i\cdot} = E^{-1}\sum_j\eta_{ij}. A missing cell is predicted from the main effects and from the interaction pattern its genotype and environment show in the other cells; the data augmentation of bayes_ammi() has already propagated its uncertainty into all parameters.

Value

An object of class c("bam_predict", "bam_result", "list"):

cells: data.table: Gen, Env, n (plots observed), observed (their mean, NA when missing), mean, sd, lower, upper.
genotypes: data.table: posterior mean, sd and HPD interval of every genotype mean.

Methods: print(), plot() (which = c("cells", "genotypes")), bam_interpret().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Examples

data(Maiz)
d <- Maiz[!(Maiz$entry == 8 & Maiz$site == 8), ]
fit <-
  bayes_ammi(
      .data   = d
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 1000
    , .k      = 2
    , .burnin = 200
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )
pr <-
  bam_predict(
      fit    = fit
    , .level = 0.95
    )
pr$cells[n == 0]

Prior distributions of the Bayesian AMMI model

Description

bam_priors() sets the hyperparameters of the prior distributions of bayes_ammi(). The defaults are the vague but proper priors of Crossa et al. (2011): normal priors with mean zero and variance 10^8 for the overall mean, the main effects and the singular values, and a scaled inverse chi-square prior with one degree of freedom for the error variance.

Usage

bam_priors(
  mu_mean = 0,
  mu_var = 1e+08,
  tau_mean = 0,
  tau_var = 1e+08,
  delta_mean = 0,
  delta_var = 1e+08,
  lambda_mean = 0,
  lambda_var = 1e+08,
  nu = 1,
  s2 = NULL
)

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

Arguments

mu_mean, mu_var

prior mean and variance of the overall mean \mu.

tau_mean, tau_var

prior mean (a scalar or a vector of length G) and variance of the genotype effects \tau_i.

delta_mean, delta_var

prior mean (scalar or length E) and variance of the environment effects \delta_j.

lambda_mean, lambda_var

prior mean and variance of the singular values \lambda_k.

nu, s2

degrees of freedom \nu and scale S^2_0 of the prior of the error variance (s2 = NULL: the least-squares residual mean square).

x

a "bam_priors" object.

...

unused.

Details

With G genotypes, E environments and K bilinear terms,

\mu \sim N(\mu_0, \sigma^2_\mu), \quad \boldsymbol{\tau} \sim N(\mathbf{m}_\tau, \sigma^2_\tau \mathbf{I}_G), \quad \boldsymbol{\delta} \sim N(\mathbf{m}_\delta, \sigma^2_\delta \mathbf{I}_E),

each main effect restricted to sum zero,

\lambda_k \sim N^{+}(\lambda_0, \sigma^2_\lambda), \quad \lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_K > 0, \qquad \sigma^2 \sim \mathrm{Inv}\textrm{-}\chi^2(\nu, S^2_0),

and spherical uniform priors for the singular vectors on the subspaces orthogonal to \mathbf{1} and to the other vectors. S^2_0 = NULL (default) takes the residual mean square of the least-squares AMMI fit, so that the prior is equally vague whatever the unit of the trait; with \nu = 1 it carries the information of a single observation. Crossa et al. (2011) report S_0 = 10^8; their published posterior of \sigma^2 is reproduced with the vague default (see the vignette on validation).

Value

A list of class "bam_priors" with these elements.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Hobert, J. P. and Casella, G. (1996). The effect of improper priors on Gibbs sampling in hierarchical linear mixed models. Journal of the American Statistical Association, 91(436), 1461–1473.

Examples

bam_priors(
    mu_mean     = 0
  , mu_var      = 1e8
  , tau_mean    = 0
  , tau_var     = 1e8
  , delta_mean  = 0
  , delta_var   = 1e8
  , lambda_mean = 0
  , lambda_var  = 1e8
  , nu          = 1
  , s2          = NULL
  )

Posterior distribution of the ranks of the genotypes

Description

bam_ranks() gives the posterior probability that each genotype takes each rank, by its mean over the environments (the marginal value) or within one environment, with the expected rank and its HPD interval.

Usage

bam_ranks(fit, higher_better = TRUE, .env = NULL)

Arguments

fit

a result of bayes_ammi() or bayes_stability() (bam_stability() needs bayes_ammi()).

higher_better

TRUE (default) when larger values of the trait are better.

.env

NULL (default: marginal values) or the label of an environment.

Details

In draw s the genotypes are ranked by \bar{\eta}^{(s)}_{i\cdot} (1 = best); the probability of rank r is the share of draws in which genotype i takes it. Ranks summarise the uncertainty of a selection decision better than intervals of the means, since they account for all the genotypes jointly.

Value

An object of class c("bam_ranks", "bam_result", "list"): probabilities (matrix genotypes x ranks) and summary (data.table: Gen, mean_rank, median_rank, lower, upper of the HPD interval of the rank, prob_first). Methods: print(), plot(), bam_interpret().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Dias, K. O. G., dos Santos, J. P. R., Krause, M. D., Piepho, H.-P., Guimaraes, 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. doi:10.1007/s00122-022-04041-y

Examples

data(Maiz)
fit <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 1000
    , .k      = 2
    , .burnin = 200
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )
rk <-
  bam_ranks(
      fit           = fit
    , higher_better = TRUE
    , .env          = NULL
    )
rk$summary
plot(x = rk)

Published values used by the tests and vignettes

Description

Reads inst/extdata/reference.csv: one row per published number, with its source, the R call that recomputes it with this package (call), the expression that extracts the number from the result (extract, in terms of x), the printed value and the digits to compare. The same file drives the test of published values and the vignette that reproduces the published examples. Each distinct call is evaluated once.

Usage

bam_reference(source = NULL, evaluate = FALSE, level = 1)

Arguments

source

optional regular expression; keep rows whose source matches.

evaluate

logical; evaluate the rows and add columns computed and agrees (default FALSE).

level

largest level evaluated (1 = fast rows, 2 = slow rows; default 1); rows above it get NA.

Value

A data.table with columns id, source, quantity, call, extract, value, digits, level, note (and computed, agrees).

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

Examples

ref <-
  bam_reference(
      source   = NULL
    , evaluate = FALSE
    , level    = 1
    )
ref

bam_reference(
    source   = "."
  , evaluate = TRUE
  , level    = 1
  )

Save a figure (vector PDF by default)

Description

A wrapper of ggsave with the defaults of the package: a vector PDF through cairo_pdf (embedded fonts, correct glyphs), A4 landscape (11.69 x 8.27 inches) and 600 dpi for raster formats. The device follows the file extension for other formats (.png, .svg, .tiff, ...).

Usage

bam_save(
  plot = ggplot2::last_plot(),
  filename,
  path = NULL,
  width = 11.69,
  height = 8.27,
  units = c("in", "cm", "mm"),
  dpi = 600,
  scale = 1,
  device = NULL,
  limitsize = TRUE,
  ...
)

Arguments

plot

a ggplot object, for example the result of plot() on a bayesammi result (default: the last plot).

filename

file name; its extension chooses the format.

path

directory to save in (default NULL: the working directory or the path in filename).

width, height

size (default 11.69 x 8.27, A4 landscape, in units).

units

"in" (default), "cm" or "mm".

dpi

resolution for raster formats (default 600).

scale

multiplicative scaling factor (default 1).

device

NULL (default: cairo_pdf for .pdf, otherwise chosen from the extension) or a device function or name, passed to ggsave(). A file name ending in .webp gives a WebP image (drawn as PNG, then converted without loss by the suggested package magick; about half the size of the PNG).

limitsize

logical; refuse images larger than 50 x 50 inches (default TRUE).

...

passed to ggsave.

Value

The file name, invisibly.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

Examples

library(ggplot2)
p <-
  ggplot(data = data.frame(x = 1:10, y = (1:10)^2), mapping = aes(x = x, y = y)) +
  geom_line() +
  bam_theme(
      base_size       = 11
    , base_family     = "serif"
    , ticks_inward    = TRUE
    , legend_position = "bottom"
    )
f <- tempfile(fileext = ".pdf")
bam_save(
    plot      = p
  , filename  = f
  , path      = NULL
  , width     = 11.69
  , height    = 8.27
  , units     = c("in", "cm", "mm")[1]
  , dpi       = 600
  , scale     = 1
  , device    = NULL
  , limitsize = TRUE
  )
bam_save(
    plot      = p
  , filename  = basename(f)
  , path      = tempdir()
  , width     = 11.69
  , height    = 8.27
  , units     = c("in", "cm", "mm")[1]
  , dpi       = 600
  , scale     = 1
  , device    = grDevices::cairo_pdf
  , limitsize = TRUE
  )
bam_save(
    plot      = p
  , filename  = tempfile(fileext = ".png")
  , path      = NULL
  , width     = 6
  , height    = 4
  , units     = c("in", "cm", "mm")[1]
  , dpi       = 300
  , scale     = 1
  , device    = NULL
  , limitsize = TRUE
  )

Priors of the Bayesian stability models

Description

bam_stab_priors() sets the prior distributions of bayes_stability(): normal priors for the overall mean and the contrasts of the genotype effects, environment effects (the environmental index) and slope deviations, and scaled inverse chi-squared priors for the error variance and the deviation variances.

Usage

bam_stab_priors(
  mu_var = NULL,
  g_var = NULL,
  h_var = NULL,
  b_var = 100,
  nu = 1,
  s2 = NULL,
  nu_d = 1,
  s2_d = NULL
)

Arguments

mu_var, g_var, h_var

prior variances of the overall mean and of the contrasts of the genotype and environment effects (NULL: 10^4 times the variance of the response).

b_var

prior variance of the contrasts of the slope deviations b_i (default 100; the slopes are 1 + b_i, without unit).

nu, s2

degrees of freedom and scale of the prior of the error variance (s2 = NULL: the within-cell mean square).

nu_d, s2_d

degrees of freedom and scale of the prior of the deviation variances (s2_d = NULL: the pooled least-squares deviation variance).

Details

\mu \sim N(0, \sigma^2_\mu), \quad \mathbf{K}_G\mathbf{g} \sim N(\mathbf{0}, \sigma^2_g\mathbf{I}), \quad \mathbf{K}_E\mathbf{h} \sim N(\mathbf{0}, \sigma^2_h\mathbf{I}), \quad \mathbf{K}_G\mathbf{b} \sim N(\mathbf{0}, \sigma^2_b\mathbf{I}),

\sigma^2 \sim \mathrm{Inv}\textrm{-}\chi^2(\nu, S^2_0), \qquad \sigma^2_{d,i} \sim \mathrm{Inv}\textrm{-}\chi^2(\nu_d, S^2_d),

where \mathbf{K}_n are the orthonormal contrasts of matrix_k(), so that the effects and the slope deviations sum to zero. Variances left NULL are set from the data: 10^4 times the variance of the response for the effects, and the least-squares estimates for S^2_0 and S^2_d, so that the priors carry the information of \nu and \nu_d observations whatever the unit of the trait.

Value

An object of class c("bam_stab_priors", "list").

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Lian, L. and de los Campos, G. (2016). FW: an R package for Finlay-Wilkinson regression that incorporates genomic/pedigree information and covariance structures between environments. G3: Genes, Genomes, Genetics, 6(3), 589–597. doi:10.1534/g3.115.026328

Examples

bam_stab_priors(
    mu_var = NULL
  , g_var  = NULL
  , h_var  = NULL
  , b_var  = 100
  , nu     = 1
  , s2     = NULL
  , nu_d   = 1
  , s2_d   = NULL
  )

Bayesian stability of the genotypes

Description

bam_stability() gives the posterior distributions of two stability measures of every genotype: the AMMI stability value (ASV) of Purchase et al. (2000) computed draw by draw, and the variance of the interaction effects over the environments (Shukla's stability variance of the model), with posterior means, HPD intervals, the probability of being the most stable genotype and the posterior expected ranks, combined with the mean into a Bayesian yield stability index.

Usage

bam_stability(fit, higher_better = TRUE)

Arguments

fit

a result of bayes_ammi() or bayes_stability() (bam_stability() needs bayes_ammi()).

higher_better

TRUE (default) when larger values of the trait are better.

Details

With the biplot scores u_{ik} = \lambda_k^{1/2}\alpha_{ik} and the sums of squares of the terms proportional to \lambda_k^2,

\mathrm{ASV}_i = \sqrt{\left(\frac{\bar{\lambda}_1^2}{\bar{\lambda}_2^2}u_{i1}\right)^2 + u_{i2}^2}, \qquad V_i = \frac{1}{E - 1}\sum_j ge_{ij}^2,

(\mathrm{ASV}_i = |u_{i1}| with one term), evaluated in every draw. The weight of the first axis uses the posterior means \bar{\lambda}_k of the singular values: the ratio of the draws themselves has no finite mean when the posterior of \lambda_2 reaches zero. The Bayesian yield stability index is the sum of the posterior expected ranks of the mean and of the ASV.

Value

An object of class c("bam_stability", "bam_result", "list") with stability (data.table sorted by the index: Gen, mean, ASV, ASV_lower, ASV_upper, variance, prob_most_stable (smallest ASV), rank_mean, rank_ASV, YSI). Methods: print(), plot(), bam_interpret().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Purchase, J. L., Hatting, H. and van Deventer, C. S. (2000). Genotype x environment interaction of winter wheat in South Africa: II. Stability analysis of yield performance. South African Journal of Plant and Soil, 17(3), 101–107.

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

Examples

data(Maiz)
fit <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 1000
    , .k      = 2
    , .burnin = 200
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )
st <-
  bam_stability(
      fit           = fit
    , higher_better = TRUE
    )
st$stability

Probabilities of superior performance and stability

Description

bam_superiority() computes the probabilities of Dias et al. (2022) for recommending cultivars: the probability that a genotype is among the best selected fraction across environments and within each environment, the probability that it is among the most stable, their joint probability, and pairwise probabilities that one genotype outperforms or is more stable than another.

Usage

bam_superiority(fit, .s = 0.2, higher_better = TRUE)

Arguments

fit

a result of bayes_ammi() or bayes_stability() (bam_stability() needs bayes_ammi()).

.s

selection intensity: the proportion of genotypes selected (default 0.2).

higher_better

TRUE (default) when larger values of the trait are better.

Details

With \Omega^{(s)} the \lceil sG \rceil genotypes with the best marginal values \bar{\eta}^{(s)}_{i\cdot} in draw s (selection intensity s = .s),

P^{\mathrm{perf}}_i = \Pr(i \in \Omega \mid \mathbf{y}), \qquad P^{\mathrm{perf}}_{ij} = \Pr(i \in \Omega_j \mid \mathbf{y}),

where \Omega_j selects on \eta_{ij} within environment j. Stability is the variance over environments of the interaction effects ge_{ij} = \eta_{ij} - \bar{\eta}_{i\cdot} - \bar{\eta}_{\cdot j} + \bar{\eta}_{\cdot\cdot} (Shukla 1972): P^{\mathrm{stab}}_i is the probability that genotype i is among the \lceil sG \rceil genotypes with the smallest variance, and the joint probability is P^{\mathrm{perf}}_i P^{\mathrm{stab}}_i, as in Dias et al. (2022). Pairwise, \Pr(\bar{\eta}_{i\cdot} > \bar{\eta}_{i'\cdot} \mid \mathbf{y}) and \Pr(\mathrm{var}_i < \mathrm{var}_{i'} \mid \mathbf{y}).

Value

An object of class c("bam_superiority", "bam_result", "list"):

marginal: data.table sorted by the joint probability: Gen, mean (posterior mean of the marginal value), prob_performance, prob_stability, prob_joint.
conditional: matrix G \times E: probability of superior performance in every environment.
pairwise_performance, pairwise_stability: matrices G \times G: entry (i, i') is the probability that row genotype i is better (more stable) than i'.

Methods: print(), plot() (which = c("marginal", "conditional", "pairwise")), bam_interpret().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

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

Examples

data(Maiz)
fit <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 1000
    , .k      = 2
    , .burnin = 200
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )
sp <-
  bam_superiority(
      fit           = fit
    , .s            = 0.2
    , higher_better = TRUE
    )
sp$marginal
plot(x = sp, which = "marginal")

ggplot2 theme of bayesammi

Description

The theme used by every plot() method of bayesammi (identical in all packages of the author): theme_classic with a serif font, black axis lines, ticks pointing inwards and no grid, close to the look of TikZ/pgfplots figures in LaTeX documents. The plot() methods also use scales::pretty_breaks(n = 10) on every continuous axis. Save figures with bam_save.

Usage

bam_theme(
  base_size = 11,
  base_family = "serif",
  ticks_inward = TRUE,
  legend_position = "bottom"
)

Arguments

base_size

base font size in points (default 11).

base_family

font family (default "serif"; for example "Latin Modern Roman" when installed, to match LaTeX documents).

ticks_inward

logical; ticks point into the panel, as in pgfplots (default TRUE).

legend_position

position of the legend (default "bottom").

Value

A ggplot2 theme.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

Examples

library(ggplot2)
p <-
  ggplot(data = data.frame(x = 1:20, y = cumsum(sin(1:20))), mapping = aes(x = x, y = y)) +
  geom_line()
p +
  bam_theme(
      base_size       = 11
    , base_family     = "serif"
    , ticks_inward    = TRUE
    , legend_position = "bottom"
    )
p +
  bam_theme(
      base_size       = 12
    , base_family     = "serif"
    , ticks_inward    = FALSE
    , legend_position = "right"
    ) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))

Which genotype wins where: posterior probabilities

Description

bam_winners() gives, for every environment, the posterior probability that each genotype has the best expected response, the most probable winner, and the groups of environments with the same most probable winner (candidate mega-environments). This is the Bayesian counterpart of the which-won-where view of the GGE biplot, with a probability attached to every statement.

Usage

bam_winners(fit, higher_better = TRUE)

Arguments

fit

a result of bayes_ammi() or bayes_stability() (bam_stability() needs bayes_ammi()).

higher_better

TRUE (default) when larger values of the trait are better.

Details

P_{ij} = \Pr\left(\eta_{ij} = \max_{i'}\eta_{i'j} \mid \mathbf{y}\right) \approx \frac{1}{S}\sum_{s=1}^{S} I\left(\eta^{(s)}_{ij} = \max_{i'}\eta^{(s)}_{i'j}\right),

with the minimum instead of the maximum when higher_better = FALSE.

Value

An object of class c("bam_winners", "bam_result", "list"):

probabilities: matrix G \times E of the probabilities P_{ij}.
winners: data.table: Env, Winner (most probable), Probability, Second and Probability2 (the runner-up).
mega: data.table: groups of environments with the same most probable winner.

Methods: print(), plot(), bam_interpret().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Yan, W., Hunt, L. A., Sheng, Q. and Szlavnics, Z. (2000). Cultivar evaluation and mega-environment investigation based on the GGE biplot. Crop Science, 40(3), 597–605. doi:10.2135/cropsci2000.403597x

Examples

data(Maiz)
fit <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 1000
    , .k      = 2
    , .burnin = 200
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "GGE"
    , .het    = FALSE
    , .waic   = TRUE
    )
w <-
  bam_winners(
      fit           = fit
    , higher_better = TRUE
    )
w$winners
plot(x = w)

Bayesian estimation of the AMMI and GGE models

Description

bayes_ammi() estimates the additive main effects and multiplicative interaction (AMMI) model of a multi-environment trial by the Gibbs sampler of Crossa et al. (2011), or the GGE (sites regression, SREG) model of Perez-Elizalde et al. (2012), for balanced or unbalanced trials, with one error variance or one for every environment. The sampler is compiled (C++). Every parameter (overall mean, genotype and environment effects, singular values, genotype and environment scores of the bilinear terms, error variance) has a posterior distribution, so that the biplot carries statistical inference: highest posterior density (HPD) regions around the scores show which genotypes and environments interact significantly and which form homogeneous groups (biplots).

Usage

bayes_ammi(
  .data,
  .y,
  .gen,
  .env,
  .rep,
  .nIter = 1000L,
  .k = 2L,
  .burnin = 0L,
  .thin = 1L,
  .chains = 1L,
  .priors = bam_priors(),
  .level = 0.95,
  .seed = NULL,
  .align = TRUE,
  .model = c("AMMI", "GGE"),
  .het = FALSE,
  .waic = TRUE
)

## Default S3 method:
bayes_ammi(
  .data,
  .y,
  .gen,
  .env,
  .rep,
  .nIter = 1000L,
  .k = 2L,
  .burnin = 0L,
  .thin = 1L,
  .chains = 1L,
  .priors = bam_priors(),
  .level = 0.95,
  .seed = NULL,
  .align = TRUE,
  .model = c("AMMI", "GGE"),
  .het = FALSE,
  .waic = TRUE
)

Arguments

.data

a data frame with one row per plot.

.y, .gen, .env, .rep

the columns of the response, genotypes, environments and replicates: bare names, strings or vectors.

.nIter

number of iterations kept per chain before thinning (after the burn-in).

.k

number of bilinear terms K (default 2, as in Crossa et al. 2011).

.burnin

number of iterations discarded at the start of every chain (default 0).

.thin

keep every .thin-th iteration (default 1).

.chains

number of chains (default 1; two or more give the Gelman-Rubin diagnostic).

.priors

the priors, from bam_priors().

.level

probability of the HPD intervals (default 0.95).

.seed

integer seed for a reproducible run, or NULL; the random number state of the session is restored afterwards.

.align

keep the sign of every bilinear term (see Signs; default TRUE).

.model

"AMMI" (default) or "GGE" (sites regression).

.het

FALSE (default): one error variance; TRUE: one for every environment.

.waic

compute WAIC during sampling (default TRUE; costs one pass over the plots per kept draw).

Details

Model. The response of genotype i in replicate l of environment j is

y_{ijl} = \mu + \tau_i + \delta_j + \sum_{k=1}^{K} \lambda_k \alpha_{ik} \gamma_{jk} + e_{ijl}, \qquad e_{ijl} \sim N(0, \sigma^2),

with \sum_i \tau_i = \sum_j \delta_j = 0, \lambda_1 \ge \cdots \ge \lambda_K > 0, and genotype and environment score vectors \boldsymbol{\alpha}_k, \boldsymbol{\gamma}_k of unit length, orthogonal to \mathbf{1} and to each other (K = .k). The priors are those of bam_priors().

Conditional distributions. With \mathbf{Y} the G \times E table of cell means, B replicates, n = BGE observations, the orthonormal contrasts \mathbf{K}_G of matrix_k() (\boldsymbol{\tau} = \mathbf{K}_G^\top \mathbf{t}) and c_x = \sigma^2 / \sigma^2_x, each iteration draws

\mu \sim N\left(\frac{n\bar{y} + c_\mu \mu_0}{n + c_\mu}, \frac{\sigma^2}{n + c_\mu}\right), \qquad \mathbf{t} \sim N\left(\frac{BE\,\mathbf{K}_G\bar{\mathbf{y}}_{G} + c_\tau \mathbf{K}_G\mathbf{m}_\tau}{BE + c_\tau}, \frac{\sigma^2}{BE + c_\tau}\mathbf{I}\right),

the same for \boldsymbol{\delta} with BG, and for k = 1, \ldots, K

\lambda_k \sim N\left(\frac{B\,\boldsymbol{\alpha}_k^\top\mathbf{Y}\boldsymbol{\gamma}_k + c_\lambda \lambda_0}{B + c_\lambda}, \frac{\sigma^2}{B + c_\lambda}\right) \textrm{ on } (\lambda_{k+1}, \lambda_{k-1}),

p(\boldsymbol{\alpha}_k \mid \cdot) \propto \exp\left(\frac{B\lambda_k}{\sigma^2} \boldsymbol{\alpha}_k^\top\mathbf{Y}\boldsymbol{\gamma}_k\right), \qquad p(\boldsymbol{\gamma}_k \mid \cdot) \propto \exp\left(\frac{B\lambda_k}{\sigma^2}\boldsymbol{\gamma}_k^\top\mathbf{Y}^\top \boldsymbol{\alpha}_k\right),

von Mises-Fisher distributions on the unit sphere of the subspace orthogonal to \mathbf{1} and to the other score vectors (Crossa et al. 2011, eq. 8-9), and finally

\sigma^2 \sim \frac{\nu S^2_0 + \mathrm{SSW} + B\sum_{ij}(\bar{y}_{ij} - \hat{y}_{ij})^2}{\chi^2_{\nu + n}},

with SSW the within-cell sum of squares. The chains start at the least-squares AMMI fit (the first chain) or at dispersed values (the others).

A variance for every environment (.het = TRUE): e_{ijl} \sim N(0, \sigma^2_j) with independent priors \mathrm{Inv}\textrm{-}\chi^2(\nu, S^2_{0j}). With weights w_j = 1/\sigma^2_j the genotype effects keep a closed normal conditional (precision B\sum_j w_j); the overall mean and the environment effects are drawn jointly; \lambda_k has precision B\sum_j w_j\gamma_{jk}^2; the genotype scores stay von Mises-Fisher; the environment scores have the Bingham-von Mises-Fisher conditional

p(\boldsymbol{\gamma}_k \mid \cdot) \propto \exp\left(B\lambda_k\sum_j w_j\gamma_{jk}r_{\cdot jk} - \frac{B\lambda_k^2}{2}\sum_j w_j\gamma_{jk}^2\right),

(r_{\cdot jk} = \sum_i\alpha_{ik}r_{ijk} with r_{ijk} the cell residual without term k), drawn by a Gibbs sweep in the eigenbasis of its quadratic term (Hoff 2009); every \sigma^2_j is a scaled inverse chi-square with \nu + BG degrees of freedom.

GGE model (.model = "GGE"): y_{ijl} = \mu + \delta_j + \sum_k\lambda_k\alpha_{ik}\gamma_{jk} + e_{ijl}: the genotype main effect is part of the multiplicative terms, the genotype scores are orthogonal to \mathbf{1}, the environment scores only to each other (Perez-Elizalde et al. 2012). The conditional distributions are those above without \boldsymbol{\tau} and without the constraint \boldsymbol{\gamma}_k \perp \mathbf{1}.

Unbalanced trials. Missing plots and whole missing genotype-environment cells are completed by data augmentation: with B the largest number of replicates of a cell, every iteration first draws the missing plots from N(\eta_{ij}, \sigma^2_j), then the conditional distributions of the balanced design apply exactly (the plots are assumed missing at random). Every genotype and every environment needs some data.

Model choice. With .waic = TRUE the sampler accumulates the pointwise log-likelihood of the observed plots, giving WAIC; the deviance of every kept draw gives DIC (component fit; bam_compare()).

Signs. The pair (\boldsymbol{\alpha}_k, \boldsymbol{\gamma}_k) and its negative give the same model. With .align = TRUE (default) every draw takes the sign that agrees with the least-squares vectors, so that posterior means and HPD regions refer to one orientation of each axis; with .align = FALSE a weak axis can switch sign during the run and its posterior means shrink towards zero.

Changes from the earlier versions. The singular values were computed once and never sampled; the error variance was not returned (tau1 held the constant precision of the least-squares fit); the residuals for \sigma^2 were computed against the wrong environments; the within-cell sum of squares was a sum of variances; the covariance matrices of the main effects were square roots of variances; the scores were not drawn from their constrained conditional distributions; all \min(G, E) - 1 axes were fitted. tau1 now holds the draws of the error variance \sigma^2.

Value

An object of class c("bam_bayes_ammi", "bam_result", "list"); the components of the earlier versions come first, the draws of all chains one after the other:

mu1: data.table, column mu: draws of the overall mean.
tau1: data.table, column tau: draws of the error variance \sigma^2 (columns tau1, ..., tauE with a variance for every environment).
tao1: data.table, columns tao1, ..., taoG: draws of the genotype effects (no columns for the GGE model).
delta1: data.table, columns delta1, ...: draws of the environment effects.
lambdas1: data.table, columns lambdas1, ...: draws of the singular values.
alphas1, gammas1: matrices with columns .nIter (draw), Gen or Env (number) and alphas1, ... or gammas1, ...: draws of the genotype and environment scores.
alphas0, gammas0: data.tables: least-squares scores (the starting values).
summary: data.table: posterior mean, standard deviation, Monte Carlo standard error, HPD interval, effective sample size and potential scale reduction factor of every parameter (Tables 1 and B1 of Crossa et al. 2011).
scores: data.table: posterior means of the biplot coordinates \lambda_k^{1/2}\alpha_{ik} and \lambda_k^{1/2}\gamma_{jk}, with the mean of every genotype and environment.
draws: coda::mcmc.list of all parameters, one element per chain.
diagnostics: list: gelman (potential scale reduction factors, two or more chains), raftery (Raftery-Lewis dependence factors), heidel (Heidelberger-Welch stationarity and half-width tests), for the overall mean, the error variance, the main effects and the singular values.
ls: list: the least-squares AMMI estimates (mu, tau, delta, lambda, alpha, gamma, sigma2).
fit: list: deviance of every kept draw, Dbar, Dhat, pD, DIC, and with .waic = TRUE lppd, p_waic, WAIC, waic_se and waic_i (the contribution of every observed plot to WAIC, in the order of the data).
settings: list: iterations, burn-in, thinning, chains, number of terms, model, variances, priors (with the scale of the error variance filled in), level, seed, alignment, time used.

Methods: print(), summary(), plot() (which = c("trace", "density", "biplot", "effects")), biplots(), bam_interpret(), the decision tools (bam_predict(), bam_winners(), bam_superiority(), bam_ranks(), bam_stability()), bam_compare(), bam_ppc() and posterior::as_draws_df() (bam-draws).

Author(s)

Muhammad Yaseen (myaseen208@gmail.com), Jose Crossa (j.crossa@cgiar.org), Sergio Perez-Elizalde (sergiop@colpos.mx), Diego Jarquin (diego.jarquin@gmail.com), Jose Miguel Cotes, Kert Viele, Genzhou Liu, Paul L. Cornelius

References

Hoff, P. D. (2009). Simulation of the matrix Bingham-von Mises-Fisher distribution, with applications to multivariate and relational data. Journal of Computational and Graphical Statistics, 18(2), 438–456. doi:10.1198/jcgs.2009.07177

Wood, A. T. A. (1994). Simulation of the von Mises Fisher distribution. Communications in Statistics - Simulation and Computation, 23(1), 157–164.

Examples

data(Maiz)
fm1 <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 200
    , .k      = 2
    , .burnin = 50
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )
fm1
head(fm1$summary)
fm1$fit$WAIC
plot(x = fm1, which = "trace")
plot(x = fm1, which = "biplot")

Bayesian stability analysis: Finlay-Wilkinson, Eberhart-Russell and Shukla models

Description

bayes_stability() fits by Gibbs sampling the joint regression models of Finlay and Wilkinson (1963) and Eberhart and Russell (1966), and the heterogeneous interaction variance model of Shukla (1972), to the plots of a multi-environment trial. Every stability parameter (slope, deviation variance, stability variance, safety-first index) then has a posterior distribution, so genotypes can be compared by probabilities rather than by significance tests. The data may be unbalanced: plots or whole genotype-environment cells may be missing.

Usage

bayes_stability(
  .data,
  .y,
  .gen,
  .env,
  .rep,
  .model = c("ER", "FW", "Shukla"),
  .nIter = 1000L,
  .burnin = 0L,
  .thin = 1L,
  .chains = 1L,
  .priors = bam_stab_priors(),
  .level = 0.95,
  .alpha = 0.05,
  .seed = NULL,
  .waic = TRUE
)

Arguments

.data

data frame with one row per plot.

.y, .gen, .env, .rep

columns of the response, genotype, environment and replicate (bare names or strings).

.model

"ER" (default), "FW" or "Shukla"; see Details.

.nIter

number of iterations kept per chain before thinning (default 1000).

.burnin

number of iterations discarded at the start of every chain (default 0).

.thin

keep every .thin-th iteration (default 1).

.chains

number of chains (default 1); the first starts at the least-squares estimates, the others at dispersed values.

.priors

priors from bam_stab_priors().

.level

probability of the highest posterior density intervals (default 0.95).

.alpha

probability of the safety-first index (default 0.05, as in Eskridge 1990).

.seed

seed of the random numbers (the random number state of the session is restored afterwards).

.waic

compute WAIC during sampling (default TRUE).

Details

For genotype i, environment j and plot l,

y_{ijl} = \mu + g_i + (1 + b_i)h_j + d_{ij} + e_{ijl}, \qquad d_{ij} \sim N(0, \sigma^2_{d,i}), \quad e_{ijl} \sim N(0, \sigma^2),

with \sum_i g_i = \sum_j h_j = \sum_i b_i = 0. The environment effects h_j are the environmental index, estimated jointly with the other parameters as in the Bayesian Finlay-Wilkinson regression of Lian and de los Campos (2016) instead of being fixed at the environment means; the slope of genotype i is 1 + b_i, with average 1.

"ER": Eberhart-Russell: slopes and a deviation variance \sigma^2_{d,i} for every genotype (the Bayesian version of Oliveira et al. 2018). A stable genotype has a slope near 1 and a small deviation variance.
"FW": Finlay-Wilkinson: slopes and one deviation variance for all genotypes.
"Shukla": no slopes (b_i = 0) and an interaction variance \sigma^2_{d,i} for every genotype: Shukla's stability variance, as in the Bayesian heterogeneous variance model of Edwards and Jannink (2006).

The sampler draws the overall mean, the genotype effects, the environmental index and the slopes from their normal conditional distributions with the deviations d_{ij} integrated out (each observed cell mean then has variance \sigma^2_{d,i} + \sigma^2/n_{ij}); then the deviations; then the variances from their scaled inverse chi-squared conditionals. This partially collapsed Gibbs sampler mixes well even when the deviations are large. Empty cells contribute no data; their deviations are drawn from the prior, and bam_predict() predicts them.

stability reports, for every genotype, the posterior mean and HPD interval of its mean \mu + g_i, slope, deviation variance and the share of its variance across environments explained by the regression, R^2_i = (1 + b_i)^2 s^2_h / ((1 + b_i)^2 s^2_h + \sigma^2_{d,i}) with s^2_h the variance of the environmental index; the posterior probability that the slope exceeds 1; and the safety-first index of Eskridge (1990), \mu + g_i - z_{1-\alpha}\sqrt{V_i}, with V_i = b_i^2 s^2_h(1 - 1/E) + \sigma^2_{d,i} (Eberhart-Russell), b_i^2 s^2_h(1 - 1/E) (Finlay-Wilkinson) or s^2_h + \sigma^2_{d,i} (Shukla). When smaller values are better (higher_better = FALSE in bam_interpret() and the decision functions), the index is read the other way: \mu + g_i + z_{1-\alpha}\sqrt{V_i}, smaller better.

The result works with bam_predict(), bam_winners(), bam_superiority(), bam_ranks(), bam_indices(), bam_compare() and bam_ppc(), which use the draws of the cell means \mu + g_i + (1 + b_i)h_j + d_{ij}.

Value

An object of class c("bam_bayes_stability", "bam_result", "list"):

stability: data.table, one row per genotype: Gen; mean, mean_lower, mean_upper; slope, slope_lower, slope_upper, prob_slope_above_1; dev_var, dev_var_lower, dev_var_upper; R2; safety_first, safety_first_lower, safety_first_upper, prob_best_safety.
environments: data.table: Env, posterior mean and HPD interval of the environmental index h_j, and the observed environment mean.
summary: data.table of every parameter (mu, sigma2, g[i], h[j], slope[i], sigma2_d[i]): posterior mean, sd, Monte Carlo error, HPD interval, effective sample size, potential scale reduction factor.
draws: coda::mcmc.list of these parameters.
deviations: array (draws x genotypes x environments) of the deviations d_{ij}.
diagnostics: Gelman-Rubin, Raftery-Lewis and Heidelberger-Welch diagnostics of the overall mean, error variance, slopes and deviation variances.
fit: DIC and WAIC, as in bayes_ammi().
ls: least-squares estimates (the start of the first chain).
settings: iterations, burn-in, thinning, chains, model, priors, level, alpha, seed, time used.

Methods: print(), plot(x, which = c("slopes", "lines", "deviations", "trace", "density")), bam_interpret(), bam_cite().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Finlay, K. W. and Wilkinson, G. N. (1963). The analysis of adaptation in a plant-breeding programme. Australian Journal of Agricultural Research, 14(6), 742–754. doi:10.1071/AR9630742

Eberhart, S. A. and Russell, W. A. (1966). Stability parameters for comparing varieties. Crop Science, 6(1), 36–40. doi:10.2135/cropsci1966.0011183X000600010011x

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

Eskridge, K. M. (1990). Selection of stable cultivars using a safety-first rule. Crop Science, 30, 369–374.

Edwards, J. W. and Jannink, J.-L. (2006). Bayesian modeling of heterogeneous error and genotype x environment interaction variances. Crop Science, 46(2), 820. doi:10.2135/cropsci2005.0164

Oliveira, T. R. A., Carvalho, H. W. L., Nascimento, M., Costa, E. F. N., Amaral Junior, A. T., Gravina, G. A. and Carvalho Filho, J. L. S. (2018). The Eberhart and Russell's Bayesian method used as an instrument to select maize hybrids. Euphytica, 214, 64. doi:10.1007/s10681-018-2146-y

Examples

data(Maiz)
fs <-
  bayes_stability(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .model  = "ER"
    , .nIter  = 1000
    , .burnin = 200
    , .thin   = 1
    , .chains = 1
    , .priors = bam_stab_priors(
          mu_var = NULL
        , g_var  = NULL
        , h_var  = NULL
        , b_var  = 100
        , nu     = 1
        , s2     = NULL
        , nu_d   = 1
        , s2_d   = NULL
        )
    , .level  = 0.95
    , .alpha  = 0.05
    , .seed   = 1
    , .waic   = TRUE
    )
fs
plot(x = fs, which = "slopes")

Biplots with highest posterior density regions

Description

biplots() draws the Bayesian AMMI biplot of the first two bilinear terms with a highest posterior density (HPD) region around every genotype and environment (Crossa et al. 2011, Fig. 3 and 4). A region that contains the origin marks a genotype or environment whose interaction is not different from zero (stable); regions that overlap mark genotypes or environments that cannot be told apart, a homogeneous group.

Usage

biplots(
  model,
  burnin = 0.3,
  thin = 0.2,
  pb = 0.05,
  plot_stable = TRUE,
  plot_unstable = TRUE,
  ncolors = 5,
  .which = c("both", "genotype", "environment"),
  .grid = 100L
)

## Default S3 method:
biplots(
  model,
  burnin = 0.3,
  thin = 0.2,
  pb = 0.05,
  plot_stable = TRUE,
  plot_unstable = TRUE,
  ncolors = 5,
  .which = c("both", "genotype", "environment"),
  .grid = 100L
)

Arguments

model

the result of bayes_ammi() (with .k >= 2).

burnin

proportion of the stored draws of every chain discarded at its start (default 0.3).

thin

proportion of the remaining draws used (default 0.2: every fifth draw).

pb

one or more probabilities p: the regions hold 100(1 - p)% of the posterior (default 0.05, the 95% region; c(0.05, 0.01) draws the 95% and 99% regions of the paper). The classification uses the first value.

plot_stable, plot_unstable

draw the regions that contain the origin, and those that do not (default both).

ncolors

number of colours used for the regions (default 5).

.which

"both" (default), "genotype" or "environment".

.grid

number of grid points per axis of the density estimate (default 100).

Details

The points are the biplot coordinates of the paper,

\left(\lambda_1^{1/2}\alpha_{i1}, \lambda_2^{1/2}\alpha_{i2}\right), \qquad \left(\lambda_1^{1/2}\gamma_{j1}, \lambda_2^{1/2}\gamma_{j2}\right),

computed draw by draw. For each genotype and environment a bivariate kernel density \hat{f} is estimated from the draws (normal kernels, normal reference bandwidths); the 100(1 - p)% HPD region is \{\mathbf{x} : \hat{f}(\mathbf{x}) \ge c_p\} with c_p the p-quantile of \hat{f} at the draws, so that the region holds 100(1 - p)% of them. The region contains the origin when \hat{f}(\mathbf{0}) \ge c_p (stable). The criterion of the paper, the marginal HPD interval of the first coordinate excluding zero, is given as axis1_significant.

Changes from the earlier versions. The regions were drawn for the singular vectors without the factor \lambda_k^{1/2}, with a fixed bandwidth of 0.2 whatever the scale; a genotype counted as stable when the bounding box of its contour covered the origin.

Value

An object of class c("bam_biplots", "bam_result", "list"), with the components of the earlier version:

contorns: data.table of the region boundaries: byplot (Genotype, Environment), ID (number), name, prob (the region, for instance "95%"), piece (a region can have several pieces), x, y, stable.
biplot_data: data.table, one row per genotype and environment: byplot, ID, name, x, y (posterior means of the coordinates), stable (region contains the origin), axis1_lower, axis1_upper (marginal HPD interval of the first coordinate), axis1_significant (it excludes zero), group (sign of the first coordinate when significant: positive, negative, or not significant).
biplot: the ggplot.

Methods: print(), plot(), bam_interpret().

Author(s)

Julian Garcia Abadillo Velasco (garciaabadillo.j@ufl.edu), Diego Jarquin (diego.jarquin@gmail.com), Muhammad Yaseen (myaseen208@gmail.com)

References

Hyndman, R. J. (1996). Computing and graphing highest density regions. The American Statistician, 50(2), 120–126.

Examples

data(Maiz)
fm1 <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 1000
    , .k      = 2
    , .burnin = 200
    , .thin   = 1
    , .chains = 1
    , .priors = bam_priors(
          mu_mean     = 0
        , mu_var      = 1e+8
        , tau_mean    = 0
        , tau_var     = 1e+8
        , delta_mean  = 0
        , delta_var   = 1e+8
        , lambda_mean = 0
        , lambda_var  = 1e+8
        , nu          = 1
        , s2          = NULL
        )
    , .level  = 0.95
    , .seed   = 1
    , .align  = TRUE
    , .model  = c("AMMI", "GGE")
    , .het    = FALSE
    , .waic   = TRUE
    )
bp <-
  biplots(
      model         = fm1
    , burnin        = 0
    , thin          = 1
    , pb            = c(0.05, 0.01)
    , plot_stable   = TRUE
    , plot_unstable = TRUE
    , ncolors       = 5
    , .which        = "both"
    , .grid         = 100
    )
bp$biplot
bp$biplot_data

Environment effects

Description

e_eff() computes the environment main effects, the deviations of the environment means from the grand mean.

Usage

e_eff(.data, .y, .gen, .env)

## Default S3 method:
e_eff(.data, .y, .gen, .env)

Arguments

.data

a data frame with one row per plot.

.y, .gen, .env

the columns of the response, genotypes and environments: bare names, strings or vectors.

Details

\hat{\delta}_j = \bar{y}_{\cdot j} - \bar{y}, \qquad \bar{y}_{\cdot j} = \frac{1}{G}\sum_i \bar{y}_{ij}.

Value

An E \times 1 matrix with column EEffs (row names: environments).

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Examples

data(Maiz)
e_eff(
    .data = Maiz
  , .y    = y
  , .gen  = entry
  , .env  = site
  )

Genotype effects

Description

g_eff() computes the genotype main effects, the deviations of the genotype means from the grand mean.

Usage

g_eff(.data, .y, .gen, .env)

## Default S3 method:
g_eff(.data, .y, .gen, .env)

Arguments

.data

a data frame with one row per plot.

.y, .gen, .env

the columns of the response, genotypes and environments: bare names, strings or vectors.

Details

\hat{\tau}_i = \bar{y}_{i\cdot} - \bar{y}, \qquad \bar{y}_{i\cdot} = \frac{1}{E}\sum_j \bar{y}_{ij},

the least-squares estimates of the AMMI model (equal to the deviations of the plot means when the data are balanced).

Value

A G \times 1 matrix with column GEffs (row names: genotypes).

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Examples

data(Maiz)
g_eff(
    .data = Maiz
  , .y    = y
  , .gen  = entry
  , .env  = site
  )

Least-squares AMMI model

Description

ge_ammi() computes the least-squares estimates of the AMMI model with all interaction axes: grand mean, genotype and environment effects, singular values and vectors of the interaction, and the error variance. They are the starting values of bayes_ammi() and the classical estimates it is compared with.

Usage

ge_ammi(.data, .y, .gen, .env, .rep)

## Default S3 method:
ge_ammi(.data, .y, .gen, .env, .rep)

Arguments

.data

a data frame with one row per plot.

.y, .gen, .env, .rep

the columns of the response, genotypes, environments and replicates: bare names, strings or vectors.

Details

With the cell means \bar{y}_{ij} of B replicates, \hat{\mu} = \bar{y}, \hat{\tau}_i = \bar{y}_{i\cdot} - \bar{y}, \hat{\delta}_j = \bar{y}_{\cdot j} - \bar{y} and the singular value decomposition of the interaction effects (ge_eff()), K = \min(G, E) - 1 axes. The error variance is the residual variance of the model with random replicates within environments (ge_model()), in closed form for balanced data:

\hat{\sigma}^2 = \frac{\mathrm{SSW} - \mathrm{SS}_{B(E)}}{E(B - 1)(G - 1)}

when the replicate mean square is not smaller, otherwise (the replicate variance component estimated as zero) \mathrm{SSW}/(GE(B - 1)); SSW is the within-cell and \mathrm{SS}_{B(E)} the replicate-within-environment sum of squares.

Change from the earlier versions. The error variance came from lme4::lmer(); the closed form gives the same value for balanced data without needing lme4.

Value

A list (class c("bam_ge_ammi", "bam_result", "list")) with the components of the earlier versions:

g, e, Rep, k: numbers of genotypes, environments, replicates and axes.
ge_means: the result of ge_mean().
mu, sigma2, tau: grand mean, error variance and its reciprocal (precision).
tao, delta: genotype and environment effects (g_eff(), e_eff()).
lambdas, alphas, gammas: singular values and genotype and environment singular vectors of the interaction.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Gauch, H. G. (1992). Statistical Analysis of Regional Yield Trials: AMMI Analysis of Factorial Designs. Elsevier, Amsterdam.

Examples

data(Maiz)
fit <-
  ge_ammi(
      .data = Maiz
    , .y    = y
    , .gen  = entry
    , .env  = site
    , .rep  = rep
    )
fit$lambdas
fit$sigma2

Genotype-by-environment interaction effects

Description

ge_eff() computes the interaction effects (the cell means centred on genotype and environment means) and their singular value decomposition, the least-squares estimate of the multiplicative part of the AMMI model.

Usage

ge_eff(.data, .y, .gen, .env)

## Default S3 method:
ge_eff(.data, .y, .gen, .env)

Arguments

.data

a data frame with one row per plot.

.y, .gen, .env

the columns of the response, genotypes and environments: bare names, strings or vectors.

Details

\mathbf{Z} = (z_{ij}), \quad z_{ij} = \bar{y}_{ij} - \bar{y}_{i\cdot} - \bar{y}_{\cdot j} + \bar{y}, \qquad \mathbf{Z} = \sum_k \lambda_k\boldsymbol{\alpha}_k\boldsymbol{\gamma}_k^\top,

with singular values \lambda_1 \ge \lambda_2 \ge \cdots and orthonormal singular vectors orthogonal to \mathbf{1}.

Value

A list (class c("bam_ge_eff", "bam_result", "list")) with ge_effects (the G \times E matrix \mathbf{Z}) and ge_svd (svd() of it: d, u, v).

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Gauch, H. G. (1992). Statistical Analysis of Regional Yield Trials: AMMI Analysis of Factorial Designs. Elsevier, Amsterdam.

Examples

data(Maiz)
fit <-
  ge_eff(
      .data = Maiz
    , .y    = y
    , .gen  = entry
    , .env  = site
    )
fit$ge_svd$d

Genotype-by-environment means

Description

ge_mean() computes the table of genotype-by-environment cell means and the grand mean, the data of the AMMI model.

Usage

ge_mean(.data, .y, .gen, .env)

## Default S3 method:
ge_mean(.data, .y, .gen, .env)

Arguments

.data

a data frame with one row per plot.

.y, .gen, .env

the columns of the response, genotypes and environments: bare names, strings or vectors.

Details

With y_{ijl} the response of genotype i in replicate l of environment j and n_{ij} the number of replicates of the cell,

\bar{y}_{ij} = \frac{1}{n_{ij}}\sum_l y_{ijl}, \qquad \bar{y} = \frac{1}{GE}\sum_{ij}\bar{y}_{ij}.

Genotypes and environments are in the order of their factor levels (sorted values when the columns are not factors).

Value

A list (class c("bam_ge_mean", "bam_result", "list")) with ge_means, the G \times E matrix of means (row and column names: genotypes and environments), and grand_mean.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Examples

data(Maiz)
ge_mean(
    .data = Maiz
  , .y    = y
  , .gen  = entry
  , .env  = site
  )

Linear mixed model of a multi-environment trial

Description

ge_model() fits the two-way model with interaction and random replicates within environments by restricted maximum likelihood (lme4::lmer()). Its residual variance is the error variance of the least-squares AMMI fit (ge_ammi(), which computes it in closed form).

Usage

ge_model(.data, .y, .gen, .env, .rep)

## Default S3 method:
ge_model(.data, .y, .gen, .env, .rep)

Arguments

.data

a data frame with one row per plot.

.y, .gen, .env, .rep

the columns of the response, genotypes, environments and replicates: bare names, strings or vectors.

Details

y_{ijl} = \mu + \delta_j + \tau_i + (\tau\delta)_{ij} + b_{l(j)} + e_{ijl}, \qquad b_{l(j)} \sim N(0, \sigma^2_b), \quad e_{ijl} \sim N(0, \sigma^2).

The package lme4 is needed (it is suggested, not required, since version 1.0.0).

Value

The lmerMod object of lme4::lmer(Y ~ Env + Gen + Env:Gen + (1 | Env:Rep)).

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Bates, D., Maechler, M., Bolker, B. and Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. doi:10.18637/jss.v067.i01

Examples

data(Maiz)
if (requireNamespace("lme4", quietly = TRUE)) {
  fm <-
    ge_model(
        .data = Maiz
      , .y    = y
      , .gen  = entry
      , .env  = site
      , .rep  = rep
      )
  stats::sigma(fm)^2
}

Genotype-by-environment variances

Description

ge_var() computes the variance of the replicates of every genotype-by-environment cell.

Usage

ge_var(.data, .y, .gen, .env)

## Default S3 method:
ge_var(.data, .y, .gen, .env)

Arguments

.data

a data frame with one row per plot.

.y, .gen, .env

the columns of the response, genotypes and environments: bare names, strings or vectors.

Details

s^2_{ij} = \frac{1}{n_{ij} - 1}\sum_l (y_{ijl} - \bar{y}_{ij})^2,

NA for a cell with one replicate. The within-cell sum of squares of the Bayesian AMMI model is \sum_{ij}(n_{ij} - 1)s^2_{ij}.

Change from the earlier versions. bayes_ammi() used the sum of these variances in place of the within-cell sum of squares.

Value

A list (class c("bam_ge_var", "bam_result", "list")) with ge_variances, the G \times E matrix of variances.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Examples

data(Maiz)
ge_var(
    .data = Maiz
  , .y    = y
  , .gen  = entry
  , .env  = site
  )

Orthonormal contrast matrix

Description

matrix_k() returns the (n - 1) \times n matrix \mathbf{K}_n whose rows are orthonormal and orthogonal to \mathbf{1}_n (Helmert contrasts). bayes_ammi() writes the main effects as \boldsymbol{\tau} = \mathbf{K}_G^\top\mathbf{t}, which keeps them summing to zero.

Usage

matrix_k(n)

## Default S3 method:
matrix_k(n)

Arguments

n

number of columns (at least 2).

Details

Row r (r = 1, \ldots, n - 1) has 0 in its first r - 1 positions, (n - r)/c_r in position r and -1/c_r after it, with c_r = \sqrt{(n - r)(n - r + 1)}, so that

\mathbf{K}_n\mathbf{K}_n^\top = \mathbf{I}_{n-1}, \qquad \mathbf{K}_n\mathbf{1}_n = \mathbf{0}, \qquad \mathbf{K}_n^\top\mathbf{K}_n = \mathbf{I}_n - \frac{1}{n}\mathbf{J}_n.

Value

The (n - 1) \times n matrix \mathbf{K}_n.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Examples

K <-
  matrix_k(
    n = 4
  )
K
round(x = K %*% t(K), digits = 10)

Gram-Schmidt orthonormalisation

Description

orthnorm() orthogonalises the columns of a matrix by the Gram-Schmidt process, optionally completes them to a basis of the whole space and scales them to unit length.

Usage

orthnorm(u = NULL, basis = TRUE, norm = TRUE)

## Default S3 method:
orthnorm(u = NULL, basis = TRUE, norm = TRUE)

Arguments

u

a matrix (or vector) whose columns are orthogonalised.

basis

complete the columns to a basis of the whole space (default TRUE).

norm

scale the columns to unit length (default TRUE).

Details

For columns \mathbf{u}_1, \ldots, \mathbf{u}_m,

\mathbf{v}_1 = \mathbf{u}_1, \qquad \mathbf{v}_i = \mathbf{u}_i - \sum_{h < i}\frac{\mathbf{u}_i^\top \mathbf{v}_h}{\mathbf{v}_h^\top\mathbf{v}_h}\mathbf{v}_h,

each divided by its length when norm = TRUE. With basis = TRUE and fewer columns than rows, unit vectors are added (those least explained by the given columns) to complete a basis.

Value

The matrix of orthogonal (orthonormal with norm = TRUE) columns; NULL for u = NULL.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

Examples

V <-
  orthnorm(
      u     = cbind(1, 1:4)
    , basis = TRUE
    , norm  = TRUE
    )
round(x = crossprod(V), digits = 10)

Package {bayesammi}

bayesammi: Bayesian Estimation of the Additive Main Effects and Multiplicative Interaction Model

Description

Details

Notation

Author(s)

References

See Also

Maize multi-environment trial

Description

Usage

Format

Details

Author(s)

Source

References

Examples

The draws of a fit for posterior, bayesplot and coda

Description

Usage

Arguments

Value

Author(s)

References

Examples

Methods for bayesammi results

Description

Usage

Arguments

Value

Author(s)

Posterior distributions of AMMI-based stability indices

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Bayesian AMMI and GGE analysis in a web browser

Description

Usage

Arguments

Details

Value

Author(s)

See Also

Examples

References of a result or of the package

Description

Usage

Arguments

Value

Author(s)

Examples

Compare Bayesian AMMI and GGE models

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Posterior distributions of classical stability indices

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Interpretation of a Bayesian AMMI analysis

Description

Usage