Package {baystability}

Type:	Package
Title:	Bayesian Stability Analysis of Genotype by Environment Interaction (GEI)
Version:	1.0.0
Date:	2026-11-04
Description:	Performs the general Bayesian estimation method of linear-bilinear models for genotype x environment interaction of Perez-Elizalde, Jarquin and Crossa (2012) <doi:10.1007/s13253-011-0063-9>: the additive main effects and multiplicative interaction (AMMI) model and the sites regression (SREG, GGE) model, with the data of an earlier trial as the prior (conjugate normal and von Mises-Fisher priors weighted by its replicates) or with vague priors; a compiled Gibbs sampler; posterior summaries with highest posterior density intervals; convergence diagnostics; biplots with highest posterior density regions; DIC and WAIC; a plain-language interpretation of every result; and a built-in 'shiny' application. For Bayesian stability analysis in general (regression and variance models, stability indices, probabilities of superior performance), see the package 'bayesstability'.
License:	GPL-2
URL:	https://myaseen208.com/baystability/, https://github.com/myaseen208/baystability
BugReports:	https://github.com/myaseen208/baystability/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:	3
Config/skill/prefix:	bys
RoxygenNote:	7.3.1
NeedsCompilation:	yes
Packaged:	2026-11-10 04:14:04 UTC; myaseen208
Author:	Muhammad Yaseen [aut, cre], Diego Jarquin [aut, ctb], Sergio Perez-Elizalde [aut, ctb], Juan Burgueno [aut, ctb], Jose Crossa [aut, ctb]
Maintainer:	Muhammad Yaseen <myaseen208@gmail.com>
Repository:	CRAN
Date/Publication:	2026-11-10 09:40:14 UTC

baystability: Bayesian Stability Analysis of Genotype by Environment Interaction (GEI)

Description

baystability implements the general Bayesian estimation method of linear-bilinear models of Perez-Elizalde, Jarquin and Crossa (2012) for multi-environment trials: the AMMI model and the sites regression (SREG, GGE) model, with the data of an earlier trial as the prior or with vague priors. Every parameter has a posterior distribution, so the biplot carries inference: highest posterior density (HPD) regions show which genotypes and environments interact significantly. For Bayesian stability analysis in general (Finlay-Wilkinson, Eberhart-Russell and Shukla models, stability indices, probabilities of superior performance) see the package bayesstability; for the Bayesian AMMI model of Crossa et al. (2011) see bayesammi.

Details

Functions, in the order of an analysis:

• bys_ls_means, bys_ls_var, bys_ls_gen_effects, bys_ls_env_effects, bys_ls_ge_effects: cell means and variances, main effects, interaction effects and their singular value decomposition.

• bys_ls_ammi, bys_ls_model: least-squares AMMI and the mixed model of the error variance.

• bys_priors, bayes_ammi: vague priors and the Gibbs sampler, with an earlier trial as prior (.prior_data).

• bayes_ammi_biplot: the biplot with HPD regions.

• bys_compare: DIC and WAIC of several fits (for example numbers of terms).

• bys_interpret, bys_cite: a plain-language interpretation and the references of a result.

• bys_contrasts, bys_orthonormal: 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 bys-methods. Utilities: bys_theme, bys_save, bys_notation, bys_reference.

Data. cultivo2008, cultivo2009: grain yield of 12 maize hybrids in 25 environments in two years (Perez-Elizalde et al. 2012).

Vignettes. Start with vignette("baystability_01_getting_started", package = "baystability"); browseVignettes("baystability") 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. The full table is bys_notation().

Author(s)

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

References

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

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

See Also

bys-methods, bys_notation

Bayesian estimation of linear-bilinear models: AMMI and SREG (GGE)

Description

bayes_ammi() estimates the linear-bilinear models of a multi-environment trial by the general Bayesian method of Perez-Elizalde et al. (2012): the additive main effects and multiplicative interaction (AMMI) model and the sites regression (SREG, GGE) model. The data of an earlier trial with the same genotypes and environments can be the prior (.prior_data): its cell means, effects and singular value decomposition give conjugate priors whose weight is its number of replicates, as in the paper's analysis of two years of a maize trial. Without an earlier trial the priors are vague but proper. The sampler is compiled (C++); every parameter has a posterior distribution, and bayes_ammi_biplot() draws the biplot with highest posterior density (HPD) regions.

Usage

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

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

Arguments

Details

Models. For genotype i, replicate l and environment j,

\textrm{AMMI: } y_{ijl} = \mu + \tau_i + \delta_j + \sum_{k=1}^{K} \lambda_k \alpha_{ik} \gamma_{jk} + e_{ijl}, \qquad \textrm{SREG: } y_{ijl} = \mu + \delta_j + \sum_{k=1}^{K} \lambda_k \alpha_{ik} \gamma_{jk} + e_{ijl},

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 orthonormal score matrices \mathbf{U} = (\boldsymbol{\alpha}_k) and \mathbf{V} = (\boldsymbol{\gamma}_k), \mathbf{U}^\top\mathbf{1} = \mathbf{0}, and for AMMI also \mathbf{V}^\top\mathbf{1} = \mathbf{0} (paper, eq. 1.2-1.3). By default all \min(G, E) - 1 terms are fitted (\min(G - 1, E) for SREG), as in the paper and the earlier versions; .k fits fewer.

An earlier trial as prior (.prior_data). With n_0 its number of replicates (.n0), \mathbf{Y}_0 its table of cell means and \mu_0, \boldsymbol{\tau}_0, \boldsymbol{\delta}_0, \mathbf{U}_0\mathbf{D}_0 \mathbf{V}_0^\top its least-squares fit (paper, Section 2.2):

\mu \sim N\left(\mu_0, \frac{\sigma^2}{n_0 GE}\right), \quad \mathbf{K}_G\boldsymbol{\tau} \sim N\left(\mathbf{K}_G\boldsymbol{\tau}_0, \frac{\sigma^2}{n_0 E}\mathbf{I}\right), \quad \mathbf{K}_E\boldsymbol{\delta} \sim N\left(\mathbf{K}_E\boldsymbol{\delta}_0, \frac{\sigma^2}{n_0 G}\mathbf{I}\right),

\lambda_k \sim N(\lambda_{0k}, \sigma^2/n_0) truncated to the order (eq. 2.9), von Mises-Fisher priors p(\mathbf{U}) \propto \mathrm{etr}(n_0\mathbf{Y}_0\mathbf{V}_0 \mathbf{D}_0\mathbf{U}^\top/\sigma^2) and p(\mathbf{V}) \propto \mathrm{etr}(n_0\mathbf{Y}_0^\top\mathbf{U}_0\mathbf{D}_0\mathbf{V}^\top/\sigma^2) (eq. 2.7-2.8), and a S^2_0/\sigma^2 \sim \chi^2_a with a = n_0 - 1 and S^2_0 the within-cell mean square of the earlier trial (eq. 2.11; the values of the earlier versions of the package).

Conditional distributions (paper, Section 2.3). With \mathbf{Y} the cell means of the trial, B its replicates and n = BGE, the linear effects are normal with means (n\bar{y} + n_0GE\mu_0)/(n + n_0GE) and so on; the singular values are truncated normal with mean (B\,\boldsymbol{\alpha}_k^\top\mathbf{Y}\boldsymbol{\gamma}_k + n_0\lambda_{0k})/(B + n_0) and variance \sigma^2/(B + n_0) (eq. 2.15); the score matrices are von Mises-Fisher (eq. 2.13-2.14), drawn one column at a time on the sphere orthogonal to \mathbf{1} and the other columns (Hoff 2009) with parameters (B\lambda_k\mathbf{Y}\boldsymbol{\gamma}_k + n_0\lambda_{0k}\mathbf{Y}_0\boldsymbol{\gamma}_{0k})/\sigma^2 and the same for \mathbf{V}; and

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

As in the paper, the dependence of the priors on \sigma^2 is left out of its conditional distribution (it would involve the normalising constants of the von Mises-Fisher priors). Without an earlier trial the priors of bys_priors() apply, which do not depend on \sigma^2, and all conditionals are exact (Crossa et al. 2011).

Starting values and signs. With an earlier trial the first chain starts at its least-squares fit and every draw keeps the orientation of its score vectors (.align); otherwise the first chain starts at the least-squares fit of the trial. Other chains start at dispersed values. With all terms fitted and a prior trial that disagrees with the data, the posterior of the scores can have several modes; compare chains (.chains) and the potential scale reduction factors.

Unbalanced trials. Missing plots and cells are completed by data augmentation (see bys_priors); an earlier trial with empty cells is completed by the EM algorithm before its fit.

Changes from version 0.2.0. The function no longer reads the earlier trial from objects in the workspace (n0, mu0, tau0, ...; it stopped with "object 'n0' not found" without them): the earlier trial is the argument .prior_data. The singular values were computed once and never sampled; tau1 held a precision drawn from a residual matrix that was never updated; the scores were drawn without their constraints and orthogonalised afterwards; the columns of delta1 were named for the genotypes. tau1 now holds the draws of the error variance \sigma^2.

Value

An object of class c("bys_bayes_ammi", "bys_result", "list"); the components of version 0.2.0 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.

tao1

data.table, columns tao1, ..., taoG: draws of the genotype effects (no columns for SREG).

delta1

data.table, columns delta1, ..., deltaE: draws of the environment effects.

lambdas1

data.table, columns lambdas1, ...: draws of the singular values.

alphas1, gammas1

data.tables, one row per draw, columns alphas1, ..., alphas<GK> (genotype scores, genotypes within terms) and gammas1, ... (environment scores).

Mean_Alphas, Mean_Gammas

data.tables, columns V1, ..., VK: posterior means of the score matrices.

summary

data.table: posterior mean, standard deviation, Monte Carlo standard error, HPD interval, effective sample size and potential scale reduction factor of every parameter.

scores

data.table: posterior means of the biplot coordinates.

draws

coda::mcmc.list of all parameters.

diagnostics

Gelman-Rubin, Raftery-Lewis and Heidelberger-Welch diagnostics.

ls

least-squares fit of the trial; prior: the fit of the earlier trial used as prior.

fit

DIC and WAIC.

settings

iterations, burn-in, thinning, chains, terms, model, priors, level, seed, time used.

Methods: print(), summary(), plot(), bayes_ammi_biplot(), bys_interpret(), bys_compare(), bys_cite().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com), Diego Jarquin (diego.jarquin@gmail.com), Sergio Perez-Elizalde (sergiop@colpos.mx), Juan Burgueno, Jose Crossa (j.crossa@cgiar.org)

References

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

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

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.

See Also

bayes_ammi_biplot, bys_priors, bys_interpret, bys_ls_ammi

Examples

data(cultivo2008)
data(cultivo2009)
fm3 <-
  bayes_ammi(
      .data       = cultivo2009
    , .y          = y
    , .gen        = entry
    , .env        = site
    , .rep        = rep
    , .nIter      = 300
    , .prior_data = cultivo2008
    , .n0         = NULL
    , .model      = "AMMI"
    , .k          = NULL
    , .burnin     = 100
    , .thin       = 1
    , .chains     = 1
    , .priors     = bys_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
    , .waic       = TRUE
    )
fm3
fm3$Mean_Alphas
plot(x = fm3, which = "biplot")

AMMI and GGE biplots with highest posterior density regions

Description

bayes_ammi_biplot() and bayes_gge_biplot() draw the Bayesian AMMI or GGE (SREG) 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

bayes_ammi_biplot(
  fit,
  burnin = 0,
  thin = 0.5,
  pb = 0.05,
  plot_stable = TRUE,
  plot_unstable = TRUE,
  ncolors = 5,
  .which = c("both", "genotype", "environment"),
  .grid = 100L
)

## Default S3 method:
bayes_ammi_biplot(
  fit,
  burnin = 0,
  thin = 0.5,
  pb = 0.05,
  plot_stable = TRUE,
  plot_unstable = TRUE,
  ncolors = 5,
  .which = c("both", "genotype", "environment"),
  .grid = 100L
)

bayes_gge_biplot(
  fit,
  burnin = 0,
  thin = 0.5,
  pb = 0.05,
  plot_stable = TRUE,
  plot_unstable = TRUE,
  ncolors = 5,
  .which = c("both", "genotype", "environment"),
  .grid = 100L
)

## Default S3 method:
bayes_gge_biplot(
  fit,
  burnin = 0,
  thin = 0.5,
  pb = 0.05,
  plot_stable = TRUE,
  plot_unstable = TRUE,
  ncolors = 5,
  .which = c("both", "genotype", "environment"),
  .grid = 100L
)

Arguments

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("bys_biplots", "bys_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(), bys_interpret().

Author(s)

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

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

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

See Also

bayes_ammi, bys_interpret

Examples

data(cultivo2009)
fm1 <-
  bayes_ammi(
      .data       = cultivo2009
    , .y          = y
    , .gen        = entry
    , .env        = site
    , .rep        = rep
    , .nIter      = 1000
    , .prior_data = NULL
    , .n0         = NULL
    , .model      = "AMMI"
    , .k          = 2
    , .burnin     = 200
    , .thin       = 1
    , .chains     = 1
    , .priors     = bys_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
    , .waic       = TRUE
    )
bp <-
  bayes_ammi_biplot(
      fit           = 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

Bayesian GGE (sites regression) model

Description

bayes_gge() fits the sites regression (SREG) model behind the GGE biplot by the method of Perez-Elizalde et al. (2012): it is bayes_ammi(.model = "SREG"), under the name used by the package bayesstability. The genotype main effect is part of the multiplicative terms; an earlier trial enters as prior through .prior_data.

Usage

bayes_gge(
  .data,
  .y,
  .gen,
  .env,
  .rep,
  .nIter = 1000L,
  .prior_data = NULL,
  .n0 = NULL,
  .k = NULL,
  .burnin = 0L,
  .thin = 1L,
  .chains = 1L,
  .priors = bys_priors(),
  .level = 0.95,
  .seed = NULL,
  .align = TRUE,
  .waic = TRUE
)

## Default S3 method:
bayes_gge(
  .data,
  .y,
  .gen,
  .env,
  .rep,
  .nIter = 1000L,
  .prior_data = NULL,
  .n0 = NULL,
  .k = NULL,
  .burnin = 0L,
  .thin = 1L,
  .chains = 1L,
  .priors = bys_priors(),
  .level = 0.95,
  .seed = NULL,
  .align = TRUE,
  .waic = TRUE
)

Arguments

Value

An object of class c("bys_bayes_ammi", "bys_result", "list"); see bayes_ammi().

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

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

See Also

bayes_ammi, bayes_gge_biplot

Examples

data(cultivo2008)
data(cultivo2009)
fg <-
  bayes_gge(
      .data       = cultivo2009
    , .y          = y
    , .gen        = entry
    , .env        = site
    , .rep        = rep
    , .nIter      = 400
    , .prior_data = cultivo2008
    , .n0         = NULL
    , .k          = 2
    , .burnin     = 100
    , .thin       = 1
    , .chains     = 1
    , .priors     = bys_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
    , .waic       = TRUE
    )
fg
bp <-
  bayes_gge_biplot(
      fit           = fg
    , burnin        = 0
    , thin          = 1
    , pb            = 0.05
    , plot_stable   = TRUE
    , plot_unstable = TRUE
    , ncolors       = 5
    , .which        = "both"
    , .grid         = 100
    )
bp$biplot

Deprecated names of baystability 0.2.0

Description

These functions of baystability 0.2.0 still work in 1.0.0, with a warning, and will be removed in a later version. Use the new names, which follow those of the package stability and its Bayesian version bayesstability:

The arguments and results are unchanged.

Usage

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

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

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

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

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

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

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

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

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

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

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

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

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

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

matrix_k(n)

## Default S3 method:
matrix_k(n)

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

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

Arguments

Value

The result of the new function.

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 'bys_bayes_ammi'
as_draws_df(x, ...)

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

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

Arguments

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(cultivo2009)
fit <-
  bayes_ammi(
      .data       = cultivo2009
    , .y          = y
    , .gen        = entry
    , .env        = site
    , .rep        = rep
    , .nIter      = 500
    , .prior_data = NULL
    , .n0         = NULL
    , .model      = "AMMI"
    , .k          = 2
    , .burnin     = 100
    , .thin       = 1
    , .chains     = 2
    , .priors     = bys_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
    , .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 baystability results

Description

The analysis functions of baystability return the named list they have always returned, now with the class c("bys_<analysis>", "bys_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 'bys_result'
print(x, digits = 4L, ...)

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

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

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

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

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

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

Arguments

Value

print() returns x invisibly; summary() an object of class summary.bys_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)

Bayesian AMMI and SREG analysis in a web browser

Description

bys_app() starts the application of baystability: load a trial (the example cultivo2009 or a CSV file with one row per plot) and, optionally, an earlier trial of the same genotypes and environments as the prior (the example cultivo2008 or a CSV file), fit the Bayesian AMMI or SREG (GGE) model of Perez-Elizalde et al. (2012), check convergence, draw the biplot with HPD regions and compare models of different numbers of terms, with and without the earlier trial. Every page shows an interpretation computed from the results (bys_interpret()), the output as R prints it, formatted tables, plots and the R code that repeats the analysis, each with a Copy button.

Usage

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

bys_app_object()

Arguments

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

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

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

References

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)

See Also

bayes_ammi, bayes_ammi_biplot, bys_compare, bys_interpret

Examples

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

References of a result or of the package

Description

bys_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("baystability") cites the package and its main methods.

Usage

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

## Default S3 method:
bys_cite(x = NULL, style = c("text", "bibtex"))

Arguments

Value

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

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

Examples

bys_cite(
    x     = NULL
  , style = "text"
  )
data(cultivo2009)
g <-
  bys_ls_ammi(
      .data = cultivo2009
    , .y    = y
    , .gen  = entry
    , .env  = site
    , .rep  = rep
    )
bys_cite(
    x     = g
  , style = "bibtex"
  )

Compare Bayesian AMMI and GGE models

Description

bys_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

bys_compare(...)

## Default S3 method:
bys_compare(...)

Arguments

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("bys_compare", "bys_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(), bys_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(cultivo2009)
fits <-
  lapply(
      X   = 1:2
    , FUN = function(K) {
        bayes_ammi(
            .data       = cultivo2009
          , .y          = y
          , .gen        = entry
          , .env        = site
          , .rep        = rep
          , .nIter      = 600
          , .prior_data = NULL
          , .n0         = NULL
          , .model      = "AMMI"
          , .k          = K
          , .burnin     = 100
          , .thin       = 1
          , .chains     = 1
          , .priors     = bys_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
          , .waic       = TRUE
          )
      }
    )
bys_compare(
    AMMI1 = fits[[1]]
  , AMMI2 = fits[[2]]
  )

Orthonormal contrast matrix

Description

bys_contrasts() 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

bys_contrasts(n)

## Default S3 method:
bys_contrasts(n)

Arguments

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

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

See Also

bys_orthonormal

Examples

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

Interpretation of a Bayesian AMMI analysis

Description

bys_interpret() turns a result of bayes_ammi(), bayes_ammi_biplot() or bys_ls_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

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

## Default S3 method:
bys_interpret(x, higher_better = TRUE, top = 5L, ...)

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

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

Arguments

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 (bayes_ammi_biplot()): a genotype is stable when its HPD region contains the origin.

Value

An object of class "bys_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; bys_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

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

See Also

bayes_ammi, bayes_ammi_biplot

Examples

data(cultivo2009)
fm1 <-
  bayes_ammi(
      .data       = cultivo2009
    , .y          = y
    , .gen        = entry
    , .env        = site
    , .rep        = rep
    , .nIter      = 500
    , .prior_data = NULL
    , .n0         = NULL
    , .model      = c("AMMI", "GGE")
    , .k          = 2
    , .burnin     = 100
    , .thin       = 1
    , .chains     = 2
    , .priors     = bys_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
    , .waic       = TRUE
    )
txt <-
  bys_interpret(
      x             = fm1
    , higher_better = TRUE
    , top           = 5
    )
txt
cat(format(x = txt, style = "markdown"))

Least-squares AMMI model

Description

bys_ls_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

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

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

Arguments

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 (bys_ls_ge_effects()), K = \min(G, E) - 1 axes. The error variance is the residual variance of the model with random replicates within environments (bys_ls_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("bys_ge_ammi", "bys_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 bys_ls_means().

mu, sigma2, tau

grand mean, error variance and its reciprocal (precision).

tao, delta

genotype and environment effects (bys_ls_gen_effects(), bys_ls_env_effects()).

lambdas, alphas, gammas

singular values and genotype and environment singular vectors of the interaction.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

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

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

See Also

bayes_ammi, bys_ls_ge_effects

Examples

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

Environment effects

Description

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

Usage

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

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

Arguments

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

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

See Also

bys_ls_gen_effects, bys_ls_ge_effects

Examples

data(cultivo2009)
bys_ls_env_effects(
    .data = cultivo2009
  , .y    = y
  , .gen  = entry
  , .env  = site
  )

Genotype-by-environment interaction effects

Description

bys_ls_ge_effects() 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

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

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

Arguments

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("bys_ge_eff", "bys_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

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

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

See Also

bys_ls_ammi

Examples

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

Genotype effects

Description

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

Usage

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

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

Arguments

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

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

See Also

bys_ls_env_effects, bys_ls_ge_effects

Examples

data(cultivo2009)
bys_ls_gen_effects(
    .data = cultivo2009
  , .y    = y
  , .gen  = entry
  , .env  = site
  )

Genotype-by-environment means

Description

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

Usage

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

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

Arguments

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("bys_ge_mean", "bys_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

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

See Also

bys_ls_var, bys_ls_ge_effects, bys_ls_ammi

Examples

data(cultivo2009)
bys_ls_means(
    .data = cultivo2009
  , .y    = y
  , .gen  = entry
  , .env  = site
  )

Linear mixed model of a multi-environment trial

Description

bys_ls_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 (bys_ls_ammi(), which computes it in closed form).

Usage

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

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

Arguments

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

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

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

See Also

bys_ls_ammi

Examples

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

Genotype-by-environment variances

Description

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

Usage

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

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

Arguments

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("bys_ge_var", "bys_result", "list")) with ge_variances, the G \times E matrix of variances.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

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

See Also

bys_ls_means

Examples

data(cultivo2009)
bys_ls_var(
    .data = cultivo2009
  , .y    = y
  , .gen  = entry
  , .env  = site
  )

Notation of baystability

Description

The symbols used in the help pages and vignettes of baystability, 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

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

Arguments

Value

See format.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

See Also

Other baystability utilities: bys_reference(), bys_save(), bys_theme()

Examples

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

Gram-Schmidt orthonormalisation

Description

bys_orthonormal() 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

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

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

Arguments

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

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

See Also

bys_contrasts

Examples

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

Prior distributions of the Bayesian AMMI model

Description

bys_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

bys_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 'bys_priors'
print(x, ...)

Arguments

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 "bys_priors" with these elements.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

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

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.

See Also

bayes_ammi

Examples

bys_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
  )

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

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

Arguments

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)

See Also

Other baystability utilities: bys_notation(), bys_save(), bys_theme()

Examples

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

bys_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

bys_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,
  ...
)

## Default S3 method:
bys_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

Value

The file name, invisibly.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

See Also

Other baystability utilities: bys_notation(), bys_reference(), bys_theme()

Examples

library(ggplot2)
p <-
  ggplot(data = data.frame(x = 1:10, y = (1:10)^2), mapping = aes(x = x, y = y)) +
  geom_line() +
  bys_theme(
      base_size       = 11
    , base_family     = "serif"
    , ticks_inward    = TRUE
    , legend_position = "bottom"
    )
f <- tempfile(fileext = ".pdf")
bys_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
  )
bys_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
  )
bys_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
  )

ggplot2 theme of baystability

Description

The theme used by every plot() method of baystability (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 bys_save.

Usage

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

Arguments

Value

A ggplot2 theme.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

See Also

Other baystability utilities: bys_notation(), bys_reference(), bys_save()

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 +
  bys_theme(
      base_size       = 11
    , base_family     = "serif"
    , ticks_inward    = TRUE
    , legend_position = "bottom"
    )
p +
  bys_theme(
      base_size       = 12
    , base_family     = "serif"
    , ticks_inward    = FALSE
    , legend_position = "right"
    ) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))

Maize trials of 2008 and 2009

Description

Grain yield of 12 maize hybrids in 25 environments in two consecutive years, the multi-environment trial of Perez-Elizalde et al. (2012), Section 3.2: cultivo2008 (the first year, used as the prior) and cultivo2009 (the second year). Randomized complete blocks with 3 replicates in every environment.

Format

Data frames with 900 rows (12 genotypes x 25 environments x 3 replicates) and 4 columns:

site

environment, 1 to 25 (the same environments in both years).

rep

replicate, 1 to 3.

entry

genotype, 1 to 12.

y

grain yield (kg/ha).

Details

The paper describes two replicates; the data have three. The posterior means of its Table B.1 are, to within a posterior standard deviation of a hundredth, the averages of the least-squares estimates of the two years, as the method gives with equal weights for both years (see the vignette on the published analysis). The help page of version 0.2.0 described another data set.

Author(s)

Muhammad Yaseen (myaseen208@gmail.com)

Source

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

Examples

data(cultivo2008)
data(cultivo2009)
str(cultivo2009)
bys_ls_means(
    .data = cultivo2009
  , .y    = y
  , .gen  = entry
  , .env  = site
  )$grand_mean