Getting started with baystability

1 Abstract

baystability fits the linear-bilinear models of a multi-environment trial, AMMI and SREG (GGE), by the general Bayesian method of Perez-Elizalde et al. (2012). Its distinctive feature is that the data of an earlier trial can be the prior. This vignette analyses the second year of a maize trial with the first year as prior, from the data to the biplot with highest posterior density (HPD) regions and a written interpretation.

2 The data

cultivo2008 and cultivo2009 hold the grain yield (kg/ha) of 12 maize hybrids in the same 25 environments in two consecutive years, three replicates per environment.

data(cultivo2008)
data(cultivo2009)
str(cultivo2009)
Classes 'tbl_df', 'tbl' and 'data.frame':   900 obs. of  4 variables:
 $ site : int  1 1 1 1 1 1 1 1 1 1 ...
 $ entry: int  2 2 2 10 10 10 12 12 12 13 ...
 $ rep  : int  1 2 3 1 2 3 1 2 3 1 ...
 $ y    : num  5079 4899 5296 5353 4353 ...
 - attr(*, "spec")=List of 2
  ..$ cols   :List of 4
  .. ..$ site : list()
  .. .. ..- attr(*, "class")= chr [1:2] "collector_integer" "collector"
  .. ..$ entry: list()
  .. .. ..- attr(*, "class")= chr [1:2] "collector_integer" "collector"
  .. ..$ rep  : list()
  .. .. ..- attr(*, "class")= chr [1:2] "collector_integer" "collector"
  .. ..$ y    : list()
  .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
  ..$ default: list()
  .. ..- attr(*, "class")= chr [1:2] "collector_guess" "collector"
  ..- attr(*, "class")= chr "col_spec"
m08 <-
  bys_ls_means(
      .data = cultivo2008
    , .y    = y
    , .gen  = entry
    , .env  = site
    )
m09 <-
  bys_ls_means(
      .data = cultivo2009
    , .y    = y
    , .gen  = entry
    , .env  = site
    )
c(year_2008 = m08$grand_mean, year_2009 = m09$grand_mean)
year_2008 year_2009 
  4069.25   3627.30 

3 The least-squares AMMI fit of the second year

ls09 <-
  bys_ls_ammi(
      .data = cultivo2009
    , .y    = y
    , .gen  = entry
    , .env  = site
    , .rep  = rep
    )
round(ls09$lambdas)
 [1] 5204 3886 3740 2777 2704 2106 1882 1875 1363  974  828

4 The Bayesian fit with the first year as prior

.prior_data makes the first year the prior. Its cell means, effects and singular value decomposition become the prior means, and its replicates set the prior’s weight. All 11 bilinear terms are fitted, as in the paper.

fm <-
  bayes_ammi(
      .data       = cultivo2009
    , .y          = y
    , .gen        = entry
    , .env        = site
    , .rep        = rep
    , .nIter      = 3000
    , .prior_data = cultivo2008
    , .n0         = NULL
    , .model      = "AMMI"
    , .k          = NULL
    , .burnin     = 1000
    , .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       = 2012
    , .align      = TRUE
    , .waic       = TRUE
    )
fm
Bayesian AMMI with 11 multiplicative term(s), an earlier trial as prior (n0 = 3): 12 genotypes, 25 environments, 3 replicates 
2 chain(s) x 3000 kept draws (iterations 3000 after a burn-in of 1000, thinning 1); 11.4 s

Posterior means, standard deviations and 95% HPD intervals
     parameter       mean        sd      lower      upper   ess   rhat
        <char>      <num>     <num>      <num>      <num> <num>  <num>
 1:         mu  3.848e+03     30.93  3.789e+03    3910.31  6000 1.0004
 2:     sigma2  1.705e+06 101590.85  1.510e+06 1902454.99  5104 1.0001
 3:     tau[2]  1.370e+02    103.33 -6.953e+01     337.71  6000 0.9999
 4:    tau[10]  3.182e+02    102.04  1.224e+02     514.97  5763 0.9999
 5:    tau[12] -6.377e+01    101.81 -2.608e+02     141.86  6000 1.0000
 6:    tau[13]  3.917e+00    102.06 -1.912e+02     205.12  6164 1.0001
 7:    tau[14]  9.089e+01    101.98 -1.153e+02     284.47  6286 1.0000
 8:    tau[18]  1.820e+02    102.97 -1.753e+01     387.84  6000 1.0002
 9:    tau[19]  4.988e+02    102.13  3.000e+02     703.98  5514 0.9999
10:    tau[20]  6.379e+01    100.38 -1.272e+02     264.90  6000 0.9998
11:    tau[23]  1.255e+02    104.03 -7.588e+01     330.28  5382 1.0000
12:    tau[38] -4.785e+02    102.31 -6.797e+02    -279.61  5639 1.0001
13:    tau[40] -2.697e+02    101.76 -4.716e+02     -68.79  6981 1.0001
14:    tau[42] -6.080e+02    102.07 -7.997e+02    -405.36  6000 1.0004
15:   delta[1]  1.315e+02    153.33 -1.800e+02     421.88  6171 0.9999
16:   delta[2] -5.364e+02    149.50 -8.325e+02    -250.83  5713 0.9999
17:   delta[3]  8.658e+02    151.62  5.748e+02    1166.47  5385 1.0005
18:   delta[4] -3.809e+02    153.16 -6.794e+02     -88.67  6000 1.0005
19:   delta[5]  1.536e+03    147.78  1.263e+03    1836.33  6000 0.9999
20:   delta[6] -9.399e+01    150.80 -3.937e+02     191.12  6254 1.0001
21:   delta[7] -1.758e+03    150.99 -2.032e+03   -1444.18  5843 1.0002
22:   delta[8]  3.495e+02    151.95  5.839e+01     658.30  6130 1.0000
23:   delta[9] -1.466e+03    149.35 -1.750e+03   -1161.58  6000 1.0007
24:  delta[10] -7.463e+02    151.99 -1.031e+03    -433.34  6173 1.0001
25:  delta[11] -1.204e+03    150.09 -1.482e+03    -893.10  6202 1.0002
26:  delta[12] -1.899e+03    149.21 -2.184e+03   -1607.04  5741 1.0005
27:  delta[13] -9.696e+02    150.25 -1.266e+03    -683.74  5986 1.0008
28:  delta[14] -5.460e+02    150.42 -8.333e+02    -248.70  6000 0.9999
29:  delta[15] -2.217e+03    152.99 -2.515e+03   -1919.88  5609 0.9999
30:  delta[16]  8.578e+01    152.21 -2.110e+02     386.78  6000 0.9999
31:  delta[17]  1.605e+02    149.05 -1.443e+02     438.51  6000 1.0006
32:  delta[18]  1.844e+03    148.30  1.554e+03    2135.29  6000 1.0002
33:  delta[19] -8.668e+02    151.53 -1.154e+03    -558.32  6250 1.0002
34:  delta[20] -1.112e+03    150.02 -1.410e+03    -820.77  5576 1.0002
35:  delta[21] -1.622e+03    151.07 -1.922e+03   -1330.08  6000 0.9999
36:  delta[22]  4.917e+03    150.81  4.613e+03    5208.71  6232 1.0008
37:  delta[23]  2.864e+03    152.44  2.566e+03    3164.43  6000 1.0002
38:  delta[24]  9.098e+02    150.05  6.153e+02    1204.77  5683 1.0006
39:  delta[25]  1.753e+03    150.35  1.462e+03    2048.28  6001 1.0001
40:  lambda[1]  3.665e+03    399.09  2.886e+03    4435.63  1952 1.0018
41:  lambda[2]  3.282e+03    353.10  2.601e+03    3971.58  1609 1.0025
42:  lambda[3]  2.819e+03    314.38  2.210e+03    3436.75  1376 1.0017
43:  lambda[4]  2.500e+03    298.48  1.925e+03    3075.36  1041 1.0010
44:  lambda[5]  2.221e+03    302.41  1.665e+03    2853.30  1063 1.0015
45:  lambda[6]  1.852e+03    305.39  1.258e+03    2453.89  1303 1.0011
46:  lambda[7]  1.529e+03    294.92  9.447e+02    2104.87  1222 1.0033
47:  lambda[8]  1.233e+03    286.55  6.697e+02    1784.19  1179 1.0053
48:  lambda[9]  9.331e+02    283.57  3.983e+02    1496.26  1192 1.0036
49: lambda[10]  6.427e+02    272.62  1.472e+02    1171.91  1245 1.0026
50: lambda[11]  3.308e+02    232.07  2.213e-01     768.85  1980 1.0039
     parameter       mean        sd      lower      upper   ess   rhat
        <char>      <num>     <num>      <num>      <num> <num>  <num>

Scores of the bilinear terms: x$summary (blocks alpha, gamma); biplot coordinates: x$scores

Sources: Perez-Elizalde et al. (2012); Crossa et al. (2011); Hoff (2009);
Spiegelhalter et al. (2002); Watanabe (2010). Full references: bys_cite(x).

The posterior mean of the overall mean, 3848, lies between the two yearly means (4069 and 3627). The prior counts as much as the data, because both years have three replicates.

plot(x = fm, which = "trace", pars = c("mu", "sigma2", "lambda[1]", "lambda[2]"))

5 The biplot with HPD regions

bp <-
  bayes_ammi_biplot(
      fit           = fm
    , burnin        = 0
    , thin          = 1
    , pb            = c(0.10, 0.05)
    , plot_stable   = FALSE
    , plot_unstable = TRUE
    , ncolors       = 5
    , .which        = "both"
    , .grid         = 100
    )
bp$biplot

Only genotypes and environments whose HPD region excludes the origin are drawn with their regions: they interact significantly on the first two terms.

6 Interpretation

bys_interpret(
    x             = fm
  , higher_better = TRUE
  , top           = 5
  )
Bayesian AMMI analysis

11 bilinear term(s): lambda_1 = 3,665, lambda_2 = 3,282; the chains have
converged.

Findings:
- The chains have converged by the criteria of Crossa et al. (2011) (largest
  potential scale reduction factor 1.007); the smallest effective sample size
  of the main parameters is 1,040.
- 2 of 200 Heidelberger-Welch stationarity and half-width tests failed (about
  5% fail by chance at level 0.05).
- Overall mean 3,848 (95% HPD 3,789 to 3,910); error variance 1,705,000, a
  coefficient of variation of 33.9%.
- Genotypes significantly above average (HPD of the effect excludes zero): 10
  and 19; below average: 38, 40 and 42.
- Environments with significantly high means: 3, 5, 8, 18, 22 and 3 more;
  low: 2, 4, 7, 9, 10 and 8 more.
- Singular values with 95% HPD intervals: lambda_1 = 3,665 (2,886 to 4,436);
  lambda_2 = 3,282 (2,601 to 3,972); lambda_3 = 2,819 (2,210 to 3,437);
  lambda_4 = 2,500 (1,925 to 3,075); lambda_5 = 2,221 (1,665 to 2,853);
  lambda_6 = 1,852 (1,258 to 2,454); lambda_7 = 1,529 (944.7 to 2,105);
  lambda_8 = 1,233 (669.7 to 1,784); lambda_9 = 933.1 (398.3 to 1,496);
  lambda_10 = 642.7 (147.2 to 1,172); lambda_11 = 330.8 (0.2213 to 768.8).
- On the first bilinear term, genotypes with significantly positive scores:
  10, 20 and 42; negative: 2, 12 and 19.
- Environments with significantly positive first scores: 5 and 25; negative:
  4 and 9.
- Genotypes and environments with significant first scores of the same sign
  interact positively: the genotype yields more there than its main effect
  predicts; of opposite sign, less.

Recommendations:
- Recommend the genotypes with significant first scores for the environments
  whose scores have the same sign (bayes_ammi_biplot() shows the regions).
- Draw bayes_ammi_biplot() to see which genotypes and environments form
  homogeneous groups (overlapping HPD regions).

Cautions:
- The HPD interval of lambda_11 reaches down to near zero: the term is weakly
  supported and its scores are uncertain.

Sources:
- 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
- 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.

7 Where to go next

  • Vignette 2 gives the model, the priors from an earlier trial and the sampler.
  • Vignette 3 reproduces the published analysis of these data.
  • Vignette 4 fits the SREG (GGE) model and chooses the number of terms.
  • Vignette 5 lists the corrections to version 0.2.0 and the checks of the sampler.

8 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, 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