The SREG (GGE) model, the weight of the earlier trial and the number of terms

\[ \newcommand{\vy}{\mathbf{y}} \newcommand{\vx}{\mathbf{x}} \newcommand{\vb}{\mathbf{b}} \newcommand{\ve}{\mathbf{e}} \newcommand{\vone}{\mathbf{1}} \newcommand{\vbeta}{\boldsymbol{\beta}} \newcommand{\vmu}{\boldsymbol{\mu}} \newcommand{\mX}{\mathbf{X}} \newcommand{\mZ}{\mathbf{Z}} \newcommand{\mV}{\mathbf{V}} \newcommand{\mI}{\mathbf{I}} \newcommand{\mJ}{\mathbf{J}} \newcommand{\mSigma}{\boldsymbol{\Sigma}} \newcommand{\mY}{\mathbf{Y}} \newcommand{\mK}{\mathbf{K}} \newcommand{\mH}{\mathbf{H}} \newcommand{\va}{\boldsymbol{\alpha}} \newcommand{\vg}{\boldsymbol{\gamma}} \newcommand{\vt}{\boldsymbol{\tau}} \newcommand{\vd}{\boldsymbol{\delta}} \newcommand{\vtheta}{\boldsymbol{\theta}} \newcommand{\vv}{\mathbf{v}} \newcommand{\E}{\mathrm{E}} \newcommand{\Var}{\mathrm{Var}} \newcommand{\Cov}{\mathrm{Cov}} \newcommand{\tr}{\operatorname{tr}} \]

1 Abstract

The sites regression (SREG) model puts the genotype main effects into the bilinear terms, so its biplot shows which genotypes do best where. This vignette fits it with the first year as prior, then shows how the weight of the earlier trial changes the posterior and how WAIC chooses the number of terms.

2 The SREG model with the first year as prior

fs <-
  bayes_ammi(
      .data       = cultivo2009
    , .y          = y
    , .gen        = entry
    , .env        = site
    , .rep        = rep
    , .nIter      = 3000
    , .prior_data = cultivo2008
    , .n0         = NULL
    , .model      = "SREG"
    , .k          = 2
    , .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
    )
fs

Bayesian SREG with 2 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); 2.0 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    3848.46     30.35    3790.6    3907.58  6000 1.0002
 2:    sigma2 1686438.64 101459.21 1483150.1 1881754.16  5552 1.0007
 3:  delta[1]     128.59    146.43    -142.6     426.81  6000 1.0002
 4:  delta[2]    -533.71    147.60    -829.2    -256.99  6069 1.0000
 5:  delta[3]     866.34    149.15     564.0    1154.05  6203 1.0003
 6:  delta[4]    -380.47    147.99    -666.8     -90.91  5551 0.9999
 7:  delta[5]    1535.74    149.97    1239.4    1824.21  6219 1.0000
 8:  delta[6]     -95.27    149.14    -376.5     196.74  5845 0.9999
 9:  delta[7]   -1758.38    147.83   -2037.1   -1461.39  6313 0.9998
10:  delta[8]     350.94    150.17      64.3     650.00  6000 1.0011
11:  delta[9]   -1471.15    147.61   -1756.5   -1180.45  6000 1.0011
12: delta[10]    -745.08    147.81   -1029.4    -447.52  6167 0.9999
13: delta[11]   -1203.17    149.42   -1494.8    -912.56  5763 1.0002
14: delta[12]   -1897.55    148.98   -2192.3   -1610.74  6000 0.9999
15: delta[13]    -969.30    149.85   -1261.3    -688.02  6000 0.9999
16: delta[14]    -550.54    149.61    -847.5    -266.79  5646 0.9999
17: delta[15]   -2214.80    148.81   -2497.1   -1914.30  6000 1.0005
18: delta[16]      91.07    151.50    -209.2     381.44  6000 0.9999
19: delta[17]     161.26    150.95    -131.9     459.36  6404 1.0002
20: delta[18]    1845.05    148.78    1561.0    2134.90  6000 1.0001
21: delta[19]    -867.08    149.69   -1147.8    -567.86  6000 1.0003
22: delta[20]   -1111.97    150.32   -1387.3    -799.39  6000 1.0000
23: delta[21]   -1625.33    148.66   -1916.0   -1333.87  5454 1.0005
24: delta[22]    4918.80    148.58    4633.1    5218.63  5814 1.0001
25: delta[23]    2864.00    151.26    2561.8    3152.60  5649 1.0003
26: delta[24]     910.64    152.01     609.9    1205.86  6234 1.0002
27: delta[25]    1751.38    154.86    1432.6    2035.94  6000 1.0001
28: lambda[1]    6250.53    579.32    5139.4    7388.08  3888 1.0004
29: lambda[2]    3268.10    577.56    2222.5    4508.37  4130 0.9999
    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); Yan et al. (2000);
Hoff (2009); Spiegelhalter et al. (2002); Watanabe (2010). Full references:
bys_cite(x).

bs <-
  bayes_gge_biplot(
      fit           = fs
    , burnin        = 0
    , thin          = 1
    , pb            = 0.05
    , plot_stable   = FALSE
    , plot_unstable = TRUE
    , ncolors       = 5
    , .which        = "both"
    , .grid         = 100
    )
bs$biplot

On an SREG biplot the first axis reflects mostly the genotype means, so genotypes far out on it are above or below average rather than interacting.

3 The weight of the earlier trial

.n0 sets how many observations per cell the earlier trial counts for. By default it is its number of replicates (3).

w <- c(0.5, 3, 12)
res <-
  data.table::rbindlist(
      l = lapply(
          X   = w
        , FUN = function(n0) {
            f <-
              bayes_ammi(
                  .data       = cultivo2009
                , .y          = y
                , .gen        = entry
                , .env        = site
                , .rep        = rep
                , .nIter      = 1500
                , .prior_data = cultivo2008
                , .n0         = n0
                , .model      = "AMMI"
                , .k          = 2
                , .burnin     = 500
                , .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       = 4
                , .align      = TRUE
                , .waic       = TRUE
                )
            data.table::data.table(
                n0      = n0
              , mu      = f$summary[parameter == "mu", mean]
              , alpha1  = f$summary[block == "tau", mean][1]
              , lambda1 = f$summary[parameter == "lambda[1]", mean]
              )
          }
        )
    )
res

      n0       mu    alpha1  lambda1
   <num>    <num>     <num>    <num>
1:   0.5 3690.593  484.3266 3602.049
2:   3.0 3849.168  136.0094 3578.433
3:  12.0 3980.321 -154.0709 4484.151

As \(n_0\) grows the posterior moves from the second year (grand mean 3627) towards the first (4069). With \(n_0 = 3\) it sits halfway. A small \(n_0\) is the choice when the earlier trial was grown under different conditions.

4 How many terms?

WAIC compares fits with different numbers of terms (Watanabe 2010; Vehtari et al. 2017):

fits <-
  lapply(
      X   = c(1, 2, 3, 5)
    , FUN = function(K) {
        bayes_ammi(
            .data       = cultivo2009
          , .y          = y
          , .gen        = entry
          , .env        = site
          , .rep        = rep
          , .nIter      = 1500
          , .prior_data = cultivo2008
          , .n0         = NULL
          , .model      = "AMMI"
          , .k          = K
          , .burnin     = 500
          , .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       = 5
          , .align      = TRUE
          , .waic       = TRUE
          )
      }
    )
names(fits) <- paste0("K", c(1, 2, 3, 5))
cp <- do.call(what = bys_compare, args = fits)
cp

Comparison of Bayesian models (smaller is better) 

    model terms                  prior      DIC      pD     WAIC  p_waic
   <char> <int>                 <char>    <num>   <num>    <num>   <num>
1:     K1     1 earlier trial (n0 = 3) 15489.01 26.7424 15493.10 29.7240
2:     K2     2 earlier trial (n0 = 3) 15493.61 37.4422 15499.51 41.2472
3:     K3     3 earlier trial (n0 = 3) 15504.81 44.3373 15511.25 48.0848
4:     K5     5 earlier trial (n0 = 3) 15518.21 56.6405 15524.42 58.8980
        se lambda_last_lower lambda_last_mean   dWAIC dWAIC_se
     <num>             <num>            <num>   <num>    <num>
1: 50.1988          1863.926         3023.724  0.0000   0.0000
2: 51.5078          2249.549         3069.854  6.4048   7.3493
3: 51.6190          1456.046         2366.735 18.1503   8.8759
4: 52.0605          1093.993         1900.412 31.3123  10.7548

Sources: Spiegelhalter et al. (2002); Watanabe (2010); Vehtari et al. (2017).
Full references: bys_cite(x).

WAIC prefers K1; K2 lie within two standard errors of the difference (dWAIC_se), so the data do not separate them. With the first year as a strong prior, the singular values of every term are pulled towards that year’s, and few terms suffice for prediction. The paper chose three terms with the Bayes factor, which compares the marginal likelihoods of the models rather than their predictive accuracy.

5 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

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.

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.

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.