Posterior distributions of stability indices

\[ \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

Dozens of stability indices are in use, each computed once from the cell means of a trial. Computed on every posterior draw of the cell means instead, each index gets a posterior mean, an HPD interval and, for every genotype, a probability of being the best. bam_indices() does this for five classical indices and bam_ammi_indices() for six AMMI-based ones. The cell means come from any fit: AMMI or GGE (bayes_ammi()) or a stability model (bayes_stability()).

fit <-
  bayes_ammi(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .nIter  = 3000
    , .k      = 2
    , .burnin = 500
    , .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   = 11
    , .align  = TRUE
    , .model  = "AMMI"
    , .het    = FALSE
    , .waic   = TRUE
    )

2 Classical indices

For one draw of the cell means \(\eta_{ij}\), with interaction \(\epsilon_{ij} = \eta_{ij} - \bar{\eta}_{i\cdot} - \bar{\eta}_{\cdot j} + \bar{\eta}_{\cdot\cdot}\):

Index	Formula	Better
Wricke’s ecovalence	\(W_i = \sum_j \epsilon_{ij}^2\)	smaller
Shukla’s stability variance	\(\frac{G}{(G-2)(E-1)}W_i - \frac{\sum_k W_k}{(G-1)(G-2)(E-1)}\)	smaller
Lin and Binns’ superiority	\(\sum_j (\eta_{ij} - M_j)^2 / (2E)\), \(M_j\) the best in environment \(j\)	smaller
Annicchiarico’s reliability	\(\bar{Z}_i - z_{1-p} s_{Z_i}\), \(Z_{ij} = 100\,\eta_{ij}/\bar{\eta}_{\cdot j}\)	larger
Eskridge’s safety-first	\(\bar{\eta}_{i\cdot} - z_{1-\alpha} s_i\), \(s_i^2\) the variance across environments	larger

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

Posterior stability indices (larger values of the trait are better) 

Posterior means; rank_overall is the mean of the expected ranks over the indices (1 = best)
      Gen  wricke shukla lin_binns annicchiarico safety_first rank_overall
   <char>   <num>  <num>     <num>         <num>        <num>        <num>
1:      6 2464000 135300    118900        104.50         3049        2.570
2:      5 2855000 161800    152200        102.30         2942        3.192
3:      3 1129000  44930    392500         96.20         2930        3.213
4:      7 1112000  43790    508100         92.99         2736        4.241
5:      2 2225000 119100    612900         90.25         2732        5.370
6:      9 3395000 198200    308900         94.29         2552        5.383
7:      4 7851000 499800    245700         96.98         2579        6.125
8:      1 4611000 280600    672000         88.58         2775        6.532
9:      8 7841000 499100   1159000         80.39         2550        8.374

Most probably the best on each index:
           index    Gen prob_best      mean     lower     upper
          <fctr> <char>     <num>     <num>     <num>     <num>
1:        wricke      7     0.453 1112000.0     66.68 2938000.0
2:        shukla      7     0.453   43790.0 -36900.00  167000.0
3:     lin_binns      6     0.705  118900.0  27980.00  226900.0
4: annicchiarico      6     0.709     104.5     99.48     109.5
5:  safety_first      6     0.542    3049.0   2742.00    3354.0

All HPD intervals and probabilities: x$indices

Sources: Wricke (1962); Shukla (1972); Lin and Binns (1988); Annicchiarico
(1992); Eskridge (1990). Full references: bam_cite(x).

The indices measure different things. Wricke’s and Shukla’s indices measure stability alone and rank genotypes in the same order. Lin and Binns’ index measures distance from the best, so it mixes mean and stability. The rank correlation between Wricke’s and Lin and Binns’ posterior means is 0.08. Shukla’s estimator can be negative, as in the classical analysis; 5 of its 9 HPD intervals reach below zero, a sign that the stability variance of those genotypes cannot be told from zero.

Because the indices are computed from the cell means, without plot error, they describe the stability of the expected performance. The HPD intervals show how well 20 environments determine them; they are wide.

3 AMMI-based indices

With coordinates \(\mathrm{PC}_{ik} = \lambda_k^{1/2}\alpha_{ik}\) and the share \(\theta_k = \lambda_k^2/\sum_m \lambda_m^2\) of each fitted term:

Index	Formula	Source
ASV	\(\sqrt{(\lambda_1^2/\lambda_2^2\,\mathrm{PC}_{i1})^2 + \mathrm{PC}_{i2}^2}\)	Purchase et al. (2000)
SIPC	\(\sum_k \lvert\mathrm{PC}_{ik}\rvert\)	Sneller et al. (1997)
EV	\(\sum_k \alpha_{ik}^2 / K\)	Zobel (1994)
Za	\(\sum_k \theta_k \lvert\alpha_{ik}\rvert\)	Zali et al. (2012)
WAAS	\(\sum_k \theta_k \lvert\mathrm{PC}_{ik}\rvert\)	Olivoto et al. (2019)
WAASY	mean and WAAS rescaled to 0-100 and weighted	Olivoto et al. (2019)

ai <-
  bam_ammi_indices(
      fit           = fit
    , higher_better = TRUE
    , .weight_y     = 50
    )
ai

Posterior AMMI-based stability indices (2 term(s); WAASY with weight 50 on the mean) 

Posterior means; rank_overall is the mean of the expected ranks over the indices (1 = best)
      Gen    ASV  SIPC      EV     Za   WAAS  WAASY rank_overall
   <char>  <num> <num>   <num>  <num>  <num>  <num>        <num>
1:      3 130.10 19.92 0.05095 0.1646 11.650 67.600        2.504
2:      7  99.51 22.02 0.08482 0.1422  9.614 63.650        2.658
3:      6 210.90 27.92 0.07374 0.2580 18.510 79.290        3.737
4:      2 192.50 28.22 0.08125 0.2425 17.210 48.010        4.273
5:      5 219.60 32.11 0.10290 0.2768 19.640 73.780        4.364
6:      9 252.50 32.34 0.09467 0.3057 22.000 54.950        5.142
7:      1 296.40 37.70 0.12270 0.3590 25.850 33.960        6.476
8:      4 402.40 43.50 0.15210 0.4679 34.200 48.210        7.565
9:      8 376.60 53.51 0.23690 0.4737 33.680  5.992        8.281

Most probably the best on each index:
    index    Gen prob_best     mean     lower    upper
   <fctr> <char>     <num>    <num>     <num>    <num>
1:    ASV      7     0.583 99.51000 5.678e+00 204.1000
2:   SIPC      3     0.395 19.92000 5.479e+00  36.0500
3:     EV      3     0.385  0.05095 3.757e-06   0.1475
4:     Za      7     0.526  0.14220 1.018e-02   0.2696
5:   WAAS      7     0.551  9.61400 4.464e-01  18.5300
6:  WAASY      6     0.566 79.29000 6.224e+01 100.0000

All HPD intervals and probabilities: x$indices

Sources: Purchase et al. (2000); Sneller et al. (1997); Zobel (1994); Zali et
al. (2012); Olivoto et al. (2019). Full references: bam_cite(x).

WAASY gives the mean and the stability equal weights (.weight_y = 50). A breeder who values the mean more raises the weight:

data.table::rbindlist(
    l = lapply(
        X   = c(30, 50, 70)
      , FUN = function(wy) {
          a <-
            bam_ammi_indices(
                fit           = fit
              , higher_better = TRUE
              , .weight_y     = wy
              )
          a$indices[index == "WAASY"][order(-mean)][1:3, list(weight = wy, Gen, mean, prob_best)]
        }
      )
  )

   weight    Gen     mean   prob_best
    <num> <char>    <num>       <num>
1:     30      7 74.80483 0.319666667
2:     30      3 74.70340 0.296166667
3:     30      6 72.28332 0.216500000
4:     50      6 79.29195 0.565500000
5:     50      5 73.78358 0.297333333
6:     50      3 67.60028 0.087833333
7:     70      6 86.30058 0.696333333
8:     70      5 80.16992 0.292833333
9:     70      4 63.66496 0.006666667

4 Indices from a stability model

The same classical indices can be computed from the cell means of a bayes_stability() fit:

fs <-
  bayes_stability(
      .data   = Maiz
    , .y      = y
    , .gen    = entry
    , .env    = site
    , .rep    = rep
    , .model  = "ER"
    , .nIter  = 2000
    , .burnin = 500
    , .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   = 11
    , .waic   = TRUE
    )
bam_indices(
    fit           = fs
  , higher_better = TRUE
  , .prob         = 0.25
  , .alpha        = 0.05
  )$means

      Gen   wricke    shukla lin_binns annicchiarico safety_first rank_overall
   <char>    <num>     <num>     <num>         <num>        <num>        <num>
1:      6  1822522  86285.70  156276.3     105.29231     3155.480       1.8692
2:      3  1607810  71756.31  489402.1      95.20404     2910.646       3.4322
3:      5  3901674 226980.19  226178.4     100.68529     2741.897       4.0675
4:      9  2996713 165742.18  409454.0      94.96091     2612.972       4.6945
5:      7  2905543 159572.83  674857.3      91.22647     2895.690       4.9025
6:      2  2475026 130440.08  719596.2      89.47825     2649.955       5.3462
7:      4  6580280 408239.21  274473.6      97.27475     2493.226       5.6621
8:      1  6387406 395187.58  772477.1      85.52653     2488.491       7.3905
9:      8 10736649 689497.31 1399778.8      80.73232     3017.037       7.6353

5 Interpretation

bam_interpret(
    x             = ind
  , higher_better = TRUE
  , top           = 5
  )

Posterior stability indices

Best on average over the indices (mean expected rank): 6, 5, 3, 7 and 2.

Findings:
- Most probably the best on wricke: 7 (probability 0.45).
- Most probably the best on shukla: 7 (probability 0.45).
- Most probably the best on lin_binns: 6 (probability 0.71).
- Most probably the best on annicchiarico: 6 (probability 0.71).
- Most probably the best on safety_first: 6 (probability 0.54).
- The indices disagree on the best genotype: they weigh mean and stability
  differently.

Recommendations:
- Prefer genotypes that rank well on several indices and whose HPD intervals
  separate from the others; overlapping intervals mean the trial cannot tell
  them apart.

Cautions:
- Stability indices describe the environments of this trial; confirm over
  years before recommending.

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

6 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.

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.

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.

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

Farshadfar, E. (2008). Incorporation of AMMI stability value and grain yield in a single non-parametric index (GSI) in bread wheat. Pakistan Journal of Biological Sciences, 11(14), 1791–1796.