Inference on the biplot

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

A least-squares AMMI biplot shows point estimates only. Its Bayesian counterpart gives every genotype and environment a highest posterior density (HPD) region. A region that contains the origin marks a stable genotype or environment, regions that overlap mark a homogeneous group, and the marginal interval on the first term decides the direction of the interaction. This vignette shows how biplots() builds the regions and how to read them.

2 The coordinates

For draw \(s\) the point of genotype \(i\) is \(\left(\lambda_1^{(s)1/2}\alpha^{(s)}_{i1}, \lambda_2^{(s)1/2}\alpha^{(s)}_{i2}\right)\). The draws of a genotype form a cloud whose shape shows the uncertainty of its interaction.

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

3 HPD regions

A bivariate kernel density \(\hat{f}\) is estimated from the cloud. The \(100(1-p)\%\) HPD region is the set where \(\hat{f}\) exceeds the \(p\)-quantile of \(\hat{f}\) at the draws, so it holds \(100(1-p)\%\) of them and is the smallest such set (Hyndman 1996). It contains the origin exactly when \(\hat{f}(\mathbf{0})\) reaches that quantile.

bp <-
  biplots(
      model         = fit
    , burnin        = 0
    , thin          = 0.5
    , pb            = c(0.05, 0.01)
    , plot_stable   = TRUE
    , plot_unstable = TRUE
    , ncolors       = 5
    , .which        = "both"
    , .grid         = 100
    )
bp$biplot

4 Stable genotypes and environments

bp$biplot_data[, list(byplot, name, x, y, stable, axis1_lower, axis1_upper, group)]

         byplot   name            x            y stable  axis1_lower
         <char> <char>        <num>        <num> <lgcl>        <num>
 1:    Genotype      1 -27.39165676   5.76438212  FALSE -38.27946935
 2:    Genotype      2 -17.88925076   8.40300672  FALSE -28.32124518
 3:    Genotype      3 -11.78567779   3.78193269   TRUE -22.73533406
 4:    Genotype      4  37.11735449   0.09401861  FALSE  26.60001160
 5:    Genotype      5  20.43844153  -9.10673615  FALSE   9.03528012
 6:    Genotype      6  19.54334049  -5.78161457  FALSE   8.91289231
 7:    Genotype      7  -8.63141072   9.96597984   TRUE -19.98012491
 8:    Genotype      8 -34.47493736 -18.02866280  FALSE -45.52039819
 9:    Genotype      9  23.07379688   4.90769355  FALSE  12.32697873
10: Environment      1 -15.50657921  -5.55082820  FALSE -26.34872405
11: Environment      2   1.17617282   3.68818862   TRUE  -9.69787599
12: Environment      3 -17.92256587  -3.66236650  FALSE -28.83775413
13: Environment      4  10.29797179  -1.87419076   TRUE  -0.05076488
14: Environment      5   5.90081820  -1.84457566   TRUE  -5.09290073
15: Environment      6  -0.37997269   4.32224532   TRUE -10.77625345
16: Environment      7  -5.64668555  -2.18620563   TRUE -16.65112969
17: Environment      8 -46.94161098   8.79385892  FALSE -56.69780534
18: Environment      9   2.67647433   0.83060165   TRUE  -8.35287890
19: Environment     10 -15.28337774  -1.34734626  FALSE -25.39295576
20: Environment     11  27.11174574   1.51012117  FALSE  16.85410381
21: Environment     12  11.17082006   6.69505350   TRUE  -0.28560818
22: Environment     13  17.28972011   5.07043354  FALSE   5.48043438
23: Environment     14  -0.54832049   0.49671271   TRUE -11.70112079
24: Environment     15  14.55743946  -5.97775500  FALSE   3.28371897
25: Environment     16   5.66452684  -0.03517786   TRUE  -4.62622081
26: Environment     17  -0.01626159  -9.16130318   TRUE -11.30940962
27: Environment     18  15.08382629   7.31282769  FALSE   3.72421775
28: Environment     19  -6.06837963  -1.59997864   TRUE -17.33199141
29: Environment     20  -2.61576189  -5.48031542   TRUE -13.60933486
         byplot   name            x            y stable  axis1_lower
         <char> <char>        <num>        <num> <lgcl>        <num>
    axis1_upper           group
          <num>          <char>
 1:  -16.186961        negative
 2:   -6.547069        negative
 3:   -1.082809        negative
 4:   46.697197        positive
 5:   31.512507        positive
 6:   30.093179        positive
 7:    4.090398 not significant
 8:  -22.302477        negative
 9:   33.828219        positive
10:   -4.743675        negative
11:   12.000181 not significant
12:   -6.725737        negative
13:   21.354530 not significant
14:   16.566620 not significant
15:   10.981398 not significant
16:    4.768966 not significant
17:  -36.925315        negative
18:   13.265884 not significant
19:   -4.598291        negative
20:   37.375407        positive
21:   21.643977 not significant
22:   27.665171        positive
23:    9.279971 not significant
24:   25.392452        positive
25:   16.839005 not significant
26:   10.851369 not significant
27:   25.565660        positive
28:    4.525090 not significant
29:    7.968492 not significant
    axis1_upper           group
          <num>          <char>

A genotype can be stable in the joint sense and still have a first coordinate that differs from zero. This happens when the second term is so uncertain that the region stretches over the origin along the second axis.

5 Only the unstable ones

bp_u <-
  biplots(
      model         = fit
    , burnin        = 0
    , thin          = 0.5
    , pb            = 0.05
    , plot_stable   = FALSE
    , plot_unstable = TRUE
    , ncolors       = 5
    , .which        = "both"
    , .grid         = 100
    )
bp_u$biplot

6 Reading genotype and environment together

A genotype yields more than its main effects predict in environments whose first scores have the same sign as its own, and less where the signs differ. The genotypes in the positive group are therefore adapted to the environments in the positive group.

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

Biplot with highest posterior density regions

7 of 9 genotypes and 8 of 20 environments have a 95% HPD region that excludes
the origin (significant interaction).

Findings:
- Genotypes whose region contains the origin (interaction not different from
  zero, stable): 3 and 7.
- Groups on the first term (marginal HPD interval excludes zero): positive 4,
  5, 6 and 9; negative 1, 2, 3 and 8.
- Environments on the first term: positive 11, 13, 15 and 18; negative 1, 3,
  8 and 10; the others do not differ from zero.
- Overlapping regions mark genotypes (environments) whose interaction
  parameters do not differ: a homogeneous group.

Recommendations:
- Genotypes 4, 5, 6 and 9 are adapted to environments 11, 13, 15 and 18 (both
  positive on the first term).
- Genotypes 1, 2, 3 and 8 are adapted to environments 1, 3, 8 and 10 (both
  negative on the first term).

Cautions:
- The regions are estimated from the draws by a kernel density; with few
  draws their outline is rough.
- The second term is usually estimated with much more uncertainty than the
  first: long regions along the second axis.

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

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