index.Rmd

---
title: "A retrospective analysis of historical data from multi-environment trial evaluation for dry bean (*Phaseolus vulgaris* L.) in Michigan"
#author: "Written by: Leonardo Volpato"
date: "Last compiled on `r format(Sys.time(), '%d %B, %Y')`"
subtitle: "Leonardo Volpato, Francisco E. Gomez, Evan Wright, and Scott Bales"
output:
  html_document:
    theme: cerulean
    highlight: pygments
    toc: yes
    toc_depth: 4
    toc_float: yes
    number_sections: yes
    df_print: kable
    code_folding: hide
    code_download: yes
  word_document:
    reference_docx: utils/Template_Word.docx
    highlight: pygments
    toc: yes
    toc_depth: 4
    df_print: kable
    number_sections: yes
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 4
  pdf_document:
    toc: yes
    number_sections: yes
    toc_depth: 4
    latex_engine: xelatex
editor_options:
  chunk_output_type: console
# bibliography: utils/bibliography.bib
# csl: utils/apa.csl
header-includes:
  \usepackage{float}
  \usepackage{booktabs}
  \usepackage{colortbl}
  \usepackage{fancyhdr}
  \pagestyle{fancy}
  \fancyhead[C]{}
  \fancyfoot[C]{Multi-environment trials evaluation for dry bean in MI}
  \fancyfoot[L]{Volpato et al. (2023)}
  \fancyfoot[R]{\thepage}
  \renewcommand{\headrulewidth}{0.4pt}
  \renewcommand{\footrulewidth}{0.4pt}
  \usepackage{caption}
  \usepackage{threeparttable}
  \usepackage{threeparttablex}
  \usepackage{longtable}
  \usepackage{fvextra}
  \DefineVerbatimEnvironment{Highlighting}{Verbatim}{breaklines,commandchars=\\\{\}}
always_allow_html: true
tables: yes
nocite: '@*'
---

\captionsetup[table]{labelformat=empty}
\captionsetup[figure]{labelformat=empty}

```{=html}
<!-- 
!!!! IMPORTANT: run `source("utils/render.R")` to publish instead of clicking on 'Knit'
-->
```


```{r directory, message=TRUE, warning=FALSE, include=FALSE}
# Setting up the working directory {-} 
rm(list=ls())  

my.path <- dirname(rstudioapi::getActiveDocumentContext()$path)
setwd(my.path)

#remove.packages("vctrs")
#install.packages("vctrs")
#source("utils/render.R")
```

```{r setup, warning=FALSE, message=TRUE, include=FALSE}
# Set up the environment 
source("utils/config.R")
fast <- FALSE  # Make this false to skip the chunks
```


```{r badges, echo=FALSE, message=TRUE, warning=FALSE, results='asis'}

if (!knitr::is_latex_output() && knitr::is_html_output()) {
  cat(
  "[![Website](https://img.shields.io/badge/visit-website-E91E63)](https://msudrybeanbreeding.github.io/DryBean_MultiEnvTrials/)
[![Website](https://img.shields.io/badge/download-.docx-FF5722)](https://github.com/msudrybeanbreeding/DryBean_MultiEnvTrials/raw/main/word_and_pdf/SupplementaryMaterials.docx)
[![Website](https://img.shields.io/badge/see-.pdf-FF9800)](https://github.com/msudrybeanbreeding/DryBean_MultiEnvTrials/raw/main/word_and_pdf/SupplementaryMaterials_v1.pdf)")
}
```


This is the supplemental material for the manuscript **A retrospective analysis of historical data from multi-environment trial evaluation for dry bean** (*Phaseolus vulgaris* L.) **in Michigan.**

The Supplemental material is divided into six main sections:

1. The first section includes the Supplemental Tables.
2. The second section includes the Supplemental Figures.
3. The third section provides the R code needed to reproduce the descriptive analysis and methods presented in this manuscript.
4. The fourth section offers the R code to replicate the Multi-Environment Trials (MET) analysis.
5. The fifth section contains the R code for reproducing the Multi-Trait Multi-Environment analysis.
6. The sixth and final section lists the references and literature cited throughout this manuscript.  


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r message=FALSE, warning=FALSE, include=FALSE, results='asis'}
#### Packages
      library(rmarkdown)
      library(ggpmisc)
      library(ggpp)
      library(gridExtra)
      library(magrittr)
      library(Matrix)
      library(mapdata)
      library(maps)
      library(spData)
      library(asremlPlus)
      library(tidyquant)
      library(usmap)
      library(data.table)
      library(ggspatial)
      library(ggpattern)
      library(viridisLite)
      library(viridis)
      library(flextable)
      library(lubridate)
      library(raster)
      library(ggcorrplot)
      library(ggdist)
      library(tibble)
      library(RColorBrewer)
      library(rcartocolor)
      library(metan)
      library(ggstatsplot)
      library(sp)
      library(sf)
      library(ggforce)
      library(patchwork)
      library(PerformanceAnalytics)
      library(broom)
      library(quantmod)
      library(xts)
      library(openxlsx)
      library(asreml)
      library(naniar)
      library(TTR)
      library(statgenGxE)
      library(tigris)
      library(plyr)
      library(ggplot2)
      library(stringr)
      library(forcats)
      library(tidyverse)
      library(dplyr)
      library(purrr)
      library(readr)
      library(tidyr)
      library(cowplot)
      library(nadiv)
      library(DT)
      library(zoo)
      library(kableExtra)


```



```{r message=FALSE, warning=FALSE, include=FALSE, results='asis'}
## Data preparation 
data_beans = read.csv("data/DataBean_MET_GYv2.csv",h=T, stringsAsFactors = T)

if (knitr::is_html_output()) {

  print_table(data_beans)
  
}else{
  
  flextable(head(data_beans)) %>%
   add_footer_lines(
     c("Varieties Dry Beans data set from 2017 to 2022",
       "Header data set showing the 6 first entry")) %>%
   autofit() %>%
   add_header_lines("Dry Beans varieties trial") %>%
  theme_design2()
}

```


```{r message=FALSE, warning=FALSE, include=FALSE, results='asis'}

# Data adjustment
# All the effect columns must be as a factor to run in ASReml-r.
cols <- c("rep", "name", "loc","year", "mkt", "year_loc")
data_beans[cols] <- lapply(data_beans[cols], factor)
data_beans <- data.table(data_beans)

```


# Appendix A - Supplemental Tables


```{r loading data S1, echo=FALSE, message=TRUE, warning=FALSE, results='asis'}

data_beans_S1 = read.csv("data/SupplTableS1_v2.csv",h=T, stringsAsFactors = T)

data_beans_S1<- data_beans_S1 %>% 
  dplyr::select(-References,-Pedigree)

#data_beans_S1

if (knitr::is_html_output()) {

#print_table(data_beans) #If the user prefer to use the render table from DT R package.
print_table(data_beans_S1, caption = "Supplemental Table S1: Commercial names (ComName), code, market class (MktCl), source, release year (RelYear), yield (RangeYield) and maturity (RangeMat) range observed in the Michigan Dry Bean Performance Trials analyzed in the study. NA = Missing values")

# flextable(data_beans_S1) %>%
#   autofit() 
  
}else{

kable(data_beans_S1, align = "c", booktabs = T, format = "latex", 
      linesep = "", longtable = TRUE,
       caption = "Supplemental Table S1: Commercial names (ComName), code, market class (MktCl), source, release year (RelYear), yield (RangeYield) and maturity (RangeMat) range observed in the Michigan Dry Bean Performance Trials analyzed in the study.") %>%
      footnote(general = "NA = Missing values", threeparttable = TRUE, fixed_small_size = TRUE) %>%
kable_styling(latex_options = c("striped", "hold_position"), font_size = 8)
}

```

```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r loading data S2, echo=FALSE, message=TRUE, warning=FALSE, results='asis'}

data_beans_S2 = read.csv("data/SupplTableS2_v4.csv",h=T, stringsAsFactors = T)


if (knitr::is_html_output()) {
  
print_table(data_beans_S2, caption = "Supplemental Table S2: Summary of the results for random- and fixed-effects and overall means for each location of the grain yield (GY: 2017 - 2022; GY2: only 2021) and days to maturity (DM), plant height (PH) and lodging (LD) traits evaluated in the study. * Significant at P < 0.05. ** Significant at P < 0.01. *** Significant at P < 0.001. ns, nonsignificant. † LRTg and LRTge, Likelihood ratio tests for genotype and interaction gxe.  § E-MS and E/B-MS, mean squares for environment and and block-within-environment, respectively. ¶ Mean, grand mean; # Standart error. NA = Missing values. TU: Tuscola, SA: Sanilac, HU: Huron, and BA: Bay locations.")

} else {
kable(data_beans_S2, align = "c", booktabs = T, format = "latex", linesep = "", escape = FALSE,
      caption = "Supplemental Table S2: Summary of the results for random- and fixed-effects and overall means for each location of the grain yield (GY: 2017 - 2022; GY2: only 2021), days to maturity (DM), plant height (PH) and lodging (LD) traits evaluated in the study.") %>% 
      footnote(general = "* Significant at P < 0.05. ** Significant at P < 0.01. *** Significant at P < 0.001. ns, nonsignificant. † LRTg and LRTge, Likelihood ratio tests for genotype and interaction gxe. § E-MS and E/B-MS, mean squares for environment and and block-within-environment, respectively. ¶ Mean, grand mean; # Standart error. NA = Missing values. TU: Tuscola, SA: Sanilac, HU: Huron, and BA: Bay locations.", threeparttable = TRUE, fixed_small_size = TRUE) %>% 
      row_spec(c(9,18), hline_after = T) %>%
kable_styling(latex_options = c("striped", "hold_position"),font_size = 8)

}
```


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


Supplemental Table S3: Estimated variance components and genetic parameters for grain yield (GY), days to maturity (DM), plant height (PH) and lodging (LD) for 37 Black, 37 Navy and 14 Red Small beans genotypes evaluated in BA, SA, and TU environments during 2021 growing season only.

**Black Beans**

| REML                                               |       GY2       |      DM      |      PH      |      LD       |
|:--------------|:-------------:|:-------------:|:-------------:|:-------------:|
| $\mathop {\hat\sigma} \nolimits_g^2$               |  68572 (32%)^‡^  | 2.86 (39.5%) | 0.63 (10.5%) | 0.063 (17.6%) |
| $\mathop {\hat\sigma} \nolimits_{ge}^{2}$          | 17115 (8%) |  0.8 (11%)   | 0.18 (3.0%) | 0.02 (5.32%)  |
| $\mathop {\hat\sigma} \nolimits_{\varepsilon}^{2}$ | 128582 (60%) | 3.58 (49.5%) | 5.17 (86.4%) | 0.27 (77.1%)  |
| $\mathop {\hat\sigma} \nolimits_{p}^{2}$           |     214269     |     7.24     |     5.98     |     0.356     |
| $\mathop h\nolimits_g^2$                           |      0.32       |    0.395     |    0.105     |     0.176     |
| $\mathop R\nolimits_{gei}^2$                       |      0.08      |     0.11     |    0.0305    |    0.0532     |
| $\mathop h\nolimits_{gm}^2$                        |      0.81      |    0.835     |    0.561     |     0.682     |
| $Acc$                                              |      0.90       |    0.914     |    0.749     |     0.826     |
| $\mathop r\nolimits_{ge}$                          |      0.12      |    0.182     |    0.0341    |    0.0645     |
| $CVg$                                              |      8.19      |     1.84     |     4.37     |     14.5      |
| $CVr$                                              |      11.21      |     2.06     |     12.5     |     30.3      |
| $CVg/CVr$                                          |      0.73      |    0.894     |    0.349     |     0.478     |

**Navy Beans**

| REML                                               |       GY2       |      DM      |      PH      |      LD       |
|:--------------|:-------------:|:-------------:|:-------------:|:-------------:|
| $\mathop {\hat\sigma} \nolimits_g^2$               | 50195 (25.72%)^‡^  |  1.97 (32%)  | 1.36 (19.1%) | 0.195 (39.3%) |
| $\mathop {\hat\sigma} \nolimits_{ge}^{2}$          | 22779 (11.67%)  | 0.63 (10.2%) | 0.81 (11.4%) | 0.04 (7.19%)  |
| $\mathop {\hat\sigma} \nolimits_{\varepsilon}^{2}$ | 122206 (62.61%) | 3.56 (57.8%) | 4.95 (69.5%) | 0.27 (53.5%)  |
| $\mathop {\hat\sigma} \nolimits_{p}^{2}$           |     195181     |     6.17     |     7.12     |     0.496     |
| $\mathop h\nolimits_g^2$                           |      0.26       |     0.32     |    0.191     |     0.393     |
| $\mathop R\nolimits_{gei}^2$                       |      0.12      |    0.102     |    0.114     |    0.0719     |
| $\mathop h\nolimits_{gm}^2$                        |      0.74      |    0.796     |    0.666     |     0.851     |
| $Acc$                                              |      0.86      |    0.892     |    0.816     |     0.923     |
| $\mathop r\nolimits_{ge}$                          |      0.16      |     0.15     |    0.141     |     0.118     |
| $CVg$                                              |      7.45      |     1.54     |     6.35     |     23.6      |
| $CVr$                                              |      11.62      |     2.07     |     12.1     |     27.6      |
| $CVg/CVr$                                          |      0.64      |    0.744     |    0.524     |     0.856     |

```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```

**Red Beans**

| REML^§^                                            |        GY2         |      DM       |      PH       |      LD      |
|:--------------|:-------------:|:-------------:|:-------------:|:-------------:|
| $\mathop {\hat\sigma} \nolimits_g^2$               | 0.1 (0.01%)^‡^ | 0.61 (12.2%) | 1.81 (25.3%) | 0.18 (34.6%) |
| $\mathop {\hat\sigma} \nolimits_{ge}^{2}$          | 75950 (31.4%) | 0.66 (13.2%)  | 0.35 (4.9%)  | 0.03 (5.0%) |
| $\mathop {\hat\sigma} \nolimits_{\varepsilon}^{2}$ | 165822 (68.59%)  | 3.72 (74.6%) | 4.99 (69.8%) | 0.32 (60.4%) |
| $\mathop {\hat\sigma} \nolimits_{p}^{2}$           |     241772     |     4.98      |     7.142     |     0.53     |
| $\mathop h\nolimits_G^2$                           |       0.00        |     0.12      |    0.2531     |     0.35     |
| $\mathop R\nolimits_{gei}^2$                       |       0.31        |     0.13      |    0.04876    |     0.05     |
| $\mathop h\nolimits_{gm}^2$                        |       0.00        |     0.54      |    0.7728     |     0.84     |
| $Acc$                                              |       0.00        |     0.73      |    0.8791     |     0.91     |
| $\mathop r\nolimits_{ge}$                          |       0.31        |     0.15      |    0.06528    |     0.07     |
| $CVg$                                              |       0.00        |     0.86      |     7.307     |     22.9     |
| $CVr$                                              |       12.90       |     2.12      |     12.13     |     30.2     |
| $CVg/CVr$                                          |       0.00        |      0.4      |    0.6022     |     0.76     |

$\mathop {\hat\sigma} \nolimits_G^2$, genotypic variance; $\mathop {\hat\sigma} \nolimits_{GE}^{2}$ , variance of G × E interaction; $\mathop {\hat\sigma} \nolimits_{\varepsilon}^{2}$ residual variance; $\mathop {\hat\sigma} \nolimits_{P}^{2}$ phenotypic variance; $\mathop h\nolimits_G^2$ broad-sense heritability; $\mathop r\nolimits_{gei}^2$ coefficient of determination for the genotype-vs-environment interaction effects; $\mathop h\nolimits_{gm}^2$ heritability of the genotypic mean; $Acc$ accuracy of genotype selection; $\mathop r\nolimits_{GE}$ correlation between genotypic values across environments; $CVg$ genotypic coefficient of variation; $CVr$ residual coefficient of variation; $CVg/CVr$ and is the ratio between genotypic and residual coefficient of variation.^‡^ Parenthetical values indicate the percentage of the observed phenotypic variance $\mathop {\hat\sigma} \nolimits_{P}^{2}$. The genetic parameters were estimated by $\mathop h\nolimits_G^2 = \frac{\mathop {\hat\sigma} \nolimits_G^2} {\mathop {\hat\sigma} \nolimits_G^2 + \mathop {\hat\sigma} \nolimits_{GE}^2 + \mathop {\hat\sigma} \nolimits_e^2 }$ ; where $\mathop {\hat\sigma} \nolimits_G^2$ is the genotypic variance; $\mathop {\hat\sigma} \nolimits_{GE}^2$ is the genotype-by-environment interaction variance; and $\mathop {\hat\sigma} \nolimits_{\varepsilon}^{2}$ is the residual variance. Coefficient of determination of the interaction effects: $R_{gei}^2=\ \frac{{\hat{\sigma}}_{GE}^2}{{\hat{\sigma}}_G^2\ +\ {\hat{\sigma}}_{GE}^2\ +{\hat{\sigma}}_\varepsilon^2}$; heribability on the mean basis: $\mathop h\nolimits_{gm}^2 = \frac{\mathop {\hat\sigma} \nolimits_G^2}{[\mathop {\hat\sigma} \nolimits_G^2 + \mathop {\hat\sigma} \nolimits_{GE}^2 /e + \mathop {\hat\sigma} \nolimits_e^2 /\left( {eb} \right)]}$ where *e* and *b* are the number of environments and blocks, respectively; The accuracy of selection: $Acc=\sqrt{h_{GM}^2}$; genotype-environment correlation $r_{GE}=\ \frac{{\hat{\sigma}}_g^2}{{\hat{\sigma}}_G^2\ +\ {\hat{\sigma}}_{GE}^2}$; genotypic coefficient of variation $CVg = \left( {\sqrt {\mathop {\hat \sigma }\nolimits_g^2 } /\mu } \right) \times 100$ residual coefficient of variation estimated: $CVr = \left( {\sqrt {\mathop {\hat \sigma }\nolimits_e^2 } /\mu } \right) \times 100$ where $\mu$ is the grand mean; ratio between genotypic and residual coefficient of variation: ${CV}_g/{CV}_r$.


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r sd geno sel bb, echo=FALSE, message=FALSE, warning=FALSE, results='asis'}

sel_SD_BB<- readxl::read_xlsx(here::here("data", "sel_SD_bb.xlsx")) %>% 
  subset(SELECTED == "yes") %>% 
  dplyr::select(-SELECTED)  

colnames(sel_SD_BB) = c('Geno','BLUPs','Overhall','SD (%)','WAASB','Overhall','SD (%)') 

if (knitr::is_html_output()) {
  
  print_table(sel_SD_BB, caption = "Supplemental Table S4: Selection differential of selected genotypes for mean performance and stability (WAASB index) at black beans. Selected genotypes:  B14, B19, B21, B22, B25, B26, B51, B52, B53, B10, B17, B23, B61, B63, B66, B68, B2, B40, B46, B70, B24, B27, B30, B33, B60, B7.") 
  
} else{
  
kable(sel_SD_BB, align = "c", booktabs = T, format = "latex", linesep = "", 
      caption = "Supplemental Table S4: Selection differential of selected genotypes for mean performance and stability (WAASB index) at black beans. Selected genotypes:  B14, B19, B21, B22, B25, B26, B51, B52, B53, B10, B17, B23, B61, B63, B66, B68, B2, B40, B46, B70, B24, B27, B30, B33, B60, B7.") %>% 
         add_header_above(c(" "=1, "Mean performance" = 3, "Stability (WAASB)" = 3)) %>% 
         kable_styling(latex_options = c("striped", "hold_position"),font_size = 10)

}
  

```


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r sd geno sel nb, echo=FALSE, message=FALSE, warning=FALSE, results='asis'}

sel_SD_NB<- readxl::read_xlsx(here::here("data", "sel_SD_nb.xlsx")) %>% 
  subset(SELECTED == "yes") %>% 
  dplyr::select(-SELECTED)  

colnames(sel_SD_NB) = c('Geno','BLUPs','Overhall','SD (%)','WAASB','Overhall','SD (%)') 

if (knitr::is_html_output()) {
  
  print_table(sel_SD_NB, caption = "Table S5: Selection differential of selected genotypes for mean performance and stability (WAASB index) at navy beans. Selected genotypes: N37, N39, N41, N58, N64, N10, N21, N33, N3, N57, N68, N54, N69, N70, N23, N4, N44, N47, N5, N6, N60, N7, N27, N29, N59.")
  
}else{

kable(sel_SD_NB, align = "c", booktabs = T, format = "latex", linesep = "", caption = "Supplemental Table S5: Selection differential of selected genotypes for mean performance and stability (WAASB index) at navy beans. Selected genotypes: N37, N39, N41, N58, N64, N10, N21, N33, N3, N57, N68, N54, N69, N70, N23, N4, N44, N47, N5, N6, N60, N7, N27, N29, N59.") %>% 
         add_header_above(c(" "=1, "Mean performance" = 3, "Stability (WAASB)" = 3)) %>% 
         kable_styling(latex_options = c("striped", "hold_position"),font_size = 10)

}

```


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r sd geno sel sr, echo=FALSE, message=FALSE, warning=FALSE, results='asis'}

sel_SD_SR<- readxl::read_xlsx(here::here("data", "sel_SD_sr.xlsx")) %>% 
  subset(SELECTED == "yes") %>% 
  dplyr::select(-SELECTED)  

colnames(sel_SD_SR) = c('Geno','BLUPs','Overhall','SD (%)','WAASB','Overhall','SD (%)') 

if (knitr::is_html_output()) {
  
  print_table(sel_SD_SR, caption = "Supplemental Table S6: Selection differential of selected genotypes for mean performance and stability (WAASB index) at red beans. Selected genotypes:R10, R13, R18, R1, R2, R21, R11.")
  
}else{

kable(sel_SD_SR, align = "c", booktabs = T, format = "latex", linesep = "",
      caption = "Supplemental Table S6: Selection differential of selected genotypes for mean performance and stability (WAASB index) at red beans. Selected genotypes:R10, R13, R18, R1, R2, R21, R11.") %>% 
         add_header_above(c(" "=1, "Mean performance" = 3, "Stability (WAASB)" = 3)) %>% 
         kable_styling(latex_options = c("striped", "hold_position"),font_size = 10)
}

```


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r sd geno blup bb, echo=FALSE, message=FALSE, warning=FALSE, results='asis'}

blups_BB<- readxl::read_xlsx(here::here("data", "blups_bb.xlsx")) 
blups_BB$GEN <- as.character(blups_BB$GEN)
blups_BB[, -1] <- round(blups_BB[, -1], 2)

if (knitr::is_html_output()) {
  
  print_table(blups_BB, caption = "Supplemental Table S7: BLUP-based indexes for selecting genotypes with performance and stability at black beans")
  
}else{

kable(blups_BB, align = "c", booktabs = T, format = "latex", linesep = "",longtable = T,
      caption = "Supplemental Table S7: BLUP-based indexes for selecting genotypes with performance and stability at black beans.") %>% 
         kable_styling(latex_options = c("striped", "hold_position"), font_size = 6)

}

```


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r sd geno blup nb, echo=FALSE, message=FALSE, warning=FALSE, results='asis'}

blups_NB<- readxl::read_xlsx(here::here("data", "blups_nb.xlsx")) 
blups_NB$GEN <- as.character(blups_NB$GEN)
blups_NB[, -1] <- round(blups_NB[, -1], 2)

if (knitr::is_html_output()) {
  
  print_table(blups_NB, caption = "Supplemental Table S8: BLUP-based indexes for selecting genotypes with performance and stability at navy beans.")
  
}else{

kable(blups_NB, align = "c", booktabs = T, format = "latex", linesep = "",longtable = T,
      caption = "Supplemental Table S8: BLUP-based indexes for selecting genotypes with performance and stability at navy beans.") %>% 
         kable_styling(latex_options = c("striped", "hold_position"),font_size = 6)

}

```


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r sd geno blups sr, echo=FALSE, message=FALSE, warning=FALSE, results='asis'}

blups_SR<- readxl::read_xlsx(here::here("data", "blups_sr.xlsx")) 
blups_SR$GEN <- as.character(blups_SR$GEN)
blups_SR[, -1] <- round(blups_SR[, -1], 2)

if (knitr::is_html_output()) {
  
print_table(blups_SR, caption = "Supplemental Table S9: BLUP-based indexes for selecting genotypes with performance and stability at red beans.")
  
}else{
  
kable(blups_SR, align = "c", booktabs = T, format = "latex", linesep = "",longtable = T,
      caption = "Supplemental Table S9: BLUP-based indexes for selecting genotypes with performance and stability at red beans.") %>% 
         kable_styling(latex_options = c("striped", "hold_position"),font_size = 6)

}

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



Supplemental Table S10: Selection differential of the WAASBY index for 27 Black, 29 Nayv and 12 Red beans.

**Black Beans**

| Factor  |	Trait	  |  Xo^‡^	 |   Xs^§^	  |   SD (%)^¶^    |   Sense  |
|:--------------|:-------------:|:-------------:|:-------------:|:--------------:|:-------------:|
| FA 1   |   GY  	|   3,201	  | 3,258	  | 57.92 (1.81)  | Increase |
| FA 1	 |   LD	  |   1.73	  |  1.63	  | -0.1 (-5.556) | Decrease |
| FA 2   |   DM	  |   91.78  	| 91.90	  | 0.12 (0.1324) | Decrease |
| FA 2   |   PH	  |   18.15	  | 18.34	  | 0.19 (1.063)  | Increase |

**Navy Beans**		

| Factor  |	Trait	  |  Xo^‡^	 |   Xs^§^	  |   SD (%)^¶^    |   Sense  |
|:--------------|:-------------:|:-------------:|:-------------:|:--------------:|:-------------:|
| FA 1     |  	GY	 |  3,006	  |  3,145	 |  138.90 (4.619)  | Increase |
| FA 1	   |    DM	 |  91.25	  |  91.54	 |  0.29 (0.3184)   | Decrease |
| FA 1	   |    PH	 |  18.34	  |  19.69	 |  1.35 (7.371)    | Increase |
| FA 1	   |    LD	 |  1.88	  |    1.65	 |  -0.22 (-11.85)  | Decrease |

**Red Beans**	

| Factor  |	Trait	  |  Xo^‡^	 |   Xs^§^	  |   SD (%)^¶^   |   Sense  |
|:------------|:-----------:|:-------------:|:-------------:|:--------------:|:-------------:|
| FA 1	 |    PH	  |   18.25	  |  20.07	  |  1.82 (9.966)   |  Increase |
| FA 1   |    LD	  |    1.87	  |  1.91	    |  0.04 (2.104)   |  Decrease |
| FA 2	 |    GY	  |   3,156	  |  3,162	  |  5.92 (0.1875)  |  Increase |
| FA 2	 |    DM	  |   90.85	  |  90.88	  |  0.03 (0.02776) |  Decrease |

GY, grain yield; PH, plant height; DM, days to maturity; LD, lodging;^‡^ Xo, mean for WAASBY index of the original population; ^§^ Xs, mean for WAASBY index of the selected genotypes (BB: B55, B1, B29, B20, B28; NB: N38, N6, N61, N35, N52, N22; SR: R2, R13); SD, Selection differential; ^¶^% of SD in parenthesis.


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


# Appendix B - Supplemental Figures


```{r weather data1, echo=FALSE, fig.align="center", fig.height=figheight, fig.width=figwidth, fig.pos= "H", out.width="100%", message=FALSE, warning=FALSE, dpi=600, fig.cap="Supplemental Figure S1: Box plot distribution from the weather data collected by location from 2017 to 2022 cultivation years. Data were scaled to better plot visualization. Weather data obtained using  daily average value from planting to harvesting.TU: Tuscola, SA: Sanilac, HU: Huron, and BA: Bay locations."}

data_beans_w = read.csv("data/Wdata_scale_beans.csv",h=T, stringsAsFactors = T)

data_beans_w$year<- as.factor(data_beans_w$year)

plot1<- ggplot(data=data_beans_w, aes( y=data, fill=loccode )) +
  geom_boxplot()+
  theme_bw()+
  facet_wrap("variable") +
  scale_fill_grey(start = 0.9, end = 0.4) +
  #facet_grid("variable") +
  theme(axis.text.x=element_blank(),
        strip.text=element_text(face="bold"), legend.position = c(0.8, 0.08),
        legend.direction = "horizontal")+
  labs(title="", 
       subtitle="", 
       caption="Data source: https://mawn.geo.msu.edu/", 
       x=NULL, y=NULL) + 
  labs(fill = "Loc")

plot1
```
atmp_max: Max Air Temperature (1.5m), atmp_min: Min Air Temperature (1.5m), mstr0_max: Max Soil Moisture (0-30cm), mstr0_min: Min Soil Moisture (0-30cm), mstr1_max: Max Soil Moisture (30-60cm), mstr1_min: Min Soil Moisture (30-60cm), pcpn: Precipitation, relh_max: Max Relative Humidity (1.5m), relh_min: Min Relative Humidity (1.5m), rpet: Reference Potential Evapotranspiration, sden_max: Max Solar Flux, soil0_max: Max Soil Temperature (5cm), soil0_min: Min Soil Temperature (5cm), soil1_max: Max Soil Temperature (10cm), soil1_min: Min Soil Temperature (10cm), srad: Total Solar Flux, wspd_max: Max Wind Speed (3m). GGDtotal: Growing Degree-Day Calculations in Celsius (C) given by: $(max-min)/2-10$



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r weather data2, echo=FALSE, fig.align="center", fig.height=figheight, fig.width=figwidth, fig.pos= "H", out.width="100%", message=FALSE, warning=FALSE, dpi=600, fig.cap="Supplemental Figure S2: Box plot distribution from the weather data collected by year from 2017 to 2022 cultivation locations. Data were scaled to better plot visualization. Weather data obtained using  daily average value from planting to harvesting.TU: Tuscola, SA: Sanilac, HU: Huron, and BA: Bay locations."}

data_beans_w = read.csv("data/Wdata_scale_beans.csv",h=T, stringsAsFactors = T)

data_beans_w$year<- as.factor(data_beans_w$year)

plot2<- ggplot(data=data_beans_w, aes( y=data, fill=year )) +
  geom_boxplot()+
  theme_bw()+
  facet_wrap("variable") +
  scale_fill_grey(start = 0.9, end = 0.4) +
  #facet_grid("variable") +
  theme(axis.text.x=element_blank(),
        strip.text=element_text(face="bold"), legend.position = c(0.8, 0.08),
        legend.direction = "horizontal")+
  labs(title="", 
       subtitle="", 
       caption="Data source: https://mawn.geo.msu.edu/", 
       x=NULL, y=NULL) + 
  labs(fill = "Years")

plot2
```
atmp_max: Max Air Temperature (1.5m), atmp_min: Min Air Temperature (1.5m), mstr0_max: Max Soil Moisture (0-30cm), mstr0_min: Min Soil Moisture (0-30cm), mstr1_max: Max Soil Moisture (30-60cm), mstr1_min: Min Soil Moisture (30-60cm), pcpn: Precipitation, relh_max: Max Relative Humidity (1.5m), relh_min: Min Relative Humidity (1.5m), rpet: Reference Potential Evapotranspiration, sden_max: Max Solar Flux, soil0_max: Max Soil Temperature (5cm), soil0_min: Min Soil Temperature (5cm), soil1_max: Max Soil Temperature (10cm), soil1_min: Min Soil Temperature (10cm), srad: Total Solar Flux, wspd_max: Max Wind Speed (3m). GGDtotal: Growing Degree-Day Calculations in Celsius (C) given by: $(max-min)/2-10$

```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r plot_ggstatsplot1, echo=FALSE, fig.height=figheight, fig.width=figwidth, message=FALSE, fig.pos= "H", out.width="100%", warning=FALSE, dpi=600, fig.cap="Supplemental Figure S3: Combination of box and density plots for between subjects comparisons by locations (left) and year (right) of grain yield (GY in Kg/ha) for black (BB), Navy (NB) and Small Red (SR) beans. BA: Bay, HU: Huron, SA:Sanilac, TU: Tuscola locations."}

## Figure showing the results corrected across years
knitr::include_graphics(here::here("main-figures", "Plot_dist_GY.png"))

```




```{r TDscatter, eval=FALSE, fig.align="center", fig.cap="Supplemental Figure S0: Scatterplot matrix of trials for grain yield (GY). The diagonal shows the distribution and the bottom left shows the scatter points for GY in each location.", fig.height=figheight, fig.pos="H", fig.width=figwidth, message=FALSE, warning=FALSE, dpi=600, include=FALSE, out.width="100%"}

dropsTD <- statgenSTA::createTD(data = data_beans, genotype = "name", trial = "loc")
options("statgen.genoColors" = c("black", "blue", "red"))
plot(dropsTD, plotType = "scatter", traits = "gy_kg_ha", colorGenoBy = "mkt", 
     colorTrialBy = "trial", title = "Scatterplots of trials for grain yield (Kg/ha)" )

```


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r selected geno bb, echo=FALSE, warning=FALSE, message=FALSE, out.width="100%", fig.pos= "H", fig.align="center", dpi = 600, fig.cap="Supplemental Figure S4:Proportion of the phenotypic variance for 71 Black (BB), 72 Navy (NB), and 21 Red (SR) beans for grain yield in Kg/ha (GY: 2017 - 2022; GY2: only 2021) days to maturity (DM, days), plant height (PH, cm) and lodging (LD,scale) traits evaluated in the study. * Significant at P < 0.05. ** Significant at P < 0.01. *** Significant at P < 0.001. ns, nonsignificant."}

knitr::include_graphics(here::here("main-figures", "Rplot_VarCompFinal.png"))

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r plot_ggstatsplot_met_fig, echo=FALSE, fig.height=figheight, fig.width=figwidth, message=FALSE, fig.pos= "H", out.width="100%", warning=FALSE, dpi=600, fig.cap="Supplemental Figure S5: Combination of box and violin plots along with jittered data points and grand mean values for between subjects comparisons by locations of grain yield (GY) for black (BB), Navy (NB) and Small Red (SR) beans. Pairwise Games-Howell test used. Comparisons showing only significant values. BA: Bay, HU: Huron, SA: Sanilac, TU: Tuscola. Pairwise Games-Howell test used. Comparisons showing only significant between the pairs of environments."}

## Figure showing the results corrected across years
knitr::include_graphics(here::here("main-figures", "Boxviolin_blues.png"))
```


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```


```{r plot_ggstatsplot2, echo=FALSE, fig.height=figheight, fig.width=figwidth, message=FALSE, fig.pos= "H", out.width="100%", warning=FALSE, dpi=600, fig.cap="Supplemental Figure S6: Genotype’s performance across the environments for Black (BB), Navy (NB), and Small Red (SR) beans using the estimated means (BLUEs) values. BA: Bay, HU: Huron, SA:Sanilac, TU: Tuscola locations"}

## Figure showing the results corrected across years
knitr::include_graphics(here::here("main-figures", "BLUES_means.png"))

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r S7 ok,echo=FALSE, warning=FALSE, message=FALSE, out.width="100%", fig.pos= "H", fig.align="center", dpi = 600, fig.cap="Supplemental Figure S7: Genotype x environment interaction characterization of all beans market classes combined dataset from 2017 to 2022. The left panel (a) shows a GGE biplot using the adjusted grain yield (GY) phenotypic means for each location averaged across years. Mega-environments (MEs, in the divided sector by the perpendicular lines to the sides of the polygon) were built using the winning genotypes from the GGE biplot. The right top panel (b) shows the discriminativeness and representativeness among environments and in relation to (c) the ideal environment. The right bottom panel (d) shows the GEI correlation among locations. PC, principal component; BA, Bay; TU, Tuscola; HU, Huron; SA, Sanilac."}

knitr::include_graphics(here::here("main-figures", "S12_allMKT2.png"))

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r selected geno bb2, echo=FALSE, warning=FALSE, message=FALSE, out.width="100%", fig.pos= "H", fig.align="center", dpi = 600, fig.cap="Supplemental Figure S8: Mean performance (a) and stability (b) for grain yield (GY) of 72 black beans genotypes. The vertical dashed and solid lines show, respectively, the  mean of the selected genotype and the overall mean for both mean performance and WAASB index."}

knitr::include_graphics(here::here("main-figures", "Sel_GY_BB.png"))

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r selected geno nb, echo=FALSE, warning=FALSE, message=FALSE, out.width="100%", fig.pos= "H", fig.align="center", dpi = 600, fig.cap="Supplemental Figure S9: Mean performance (a) and stability (b) for grain yield (GY) of 71 navy beans genotypes. The vertical dashed and solid lines show, respectively, the  mean of the selected genotype and the overall mean for both mean performance and WAASB index."}

knitr::include_graphics(here::here("main-figures", "Sel_GY_NB.png"))

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r selected geno sr, echo=FALSE, warning=FALSE, message=FALSE, out.width="100%", fig.pos= "H", fig.align="center", dpi = 600, fig.cap="Supplemental Figure S10: Mean performance (a) and stability (b) for grain yield (GY) of 21 red beans genotypes. The vertical dashed and solid lines show, respectively, the  mean of the selected genotype and the overall mean for both mean performance and WAASB index."}

knitr::include_graphics(here::here("main-figures", "Sel_GY_SR.png"))

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r blups CI2,echo=FALSE, warning=FALSE, message=FALSE, out.width="100%", fig.align="center", fig.pos= "H", dpi = 600, fig.cap= "Supplemental Figure S11: BLUP-based stability indexes coincidence index (CI) using selection intensity of 20 top genotypes. BB, black; NB, navy; SR, red beans."}

knitr::include_graphics(here::here("main-figures", "CI_blups.png"))

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r selected waasb,echo=FALSE, warning=FALSE, message=FALSE, out.width="100%", fig.pos= "H", fig.align="center", dpi = 600, fig.cap="Supplemental Figure S12: Ranks of dry beans genotypes (a: 72 black, b:71 navys and c: 21 small red/pink) considering different weights for stability and yielding. The most-left ranks were obtained considering the stability only. The most right-ranks were obtained considering the grain yield only. Between the extremes, ranks were obtained different weights for stability and yielding. The four clusters represent four classes of genotypes: (1) Poorly productive and unstable genotypes; (2) productive but unstable genotypes; (3) stable but poorly productive genotypes; and (4), highly productive and stable genotypes. The ranks highlighted by a black rectangle are the same as those BLUPs predicted to WAASBY index and the red rectangle box are the selected genotypes by the Multi-Trait Stability Index (MTSI)."}

knitr::include_graphics(here::here("main-figures", "waasb_ratio_new.png"))

```




```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r s13,echo=FALSE, warning=FALSE, message=FALSE, out.width="100%", fig.pos= "H", fig.align="center", dpi = 600, fig.cap="Supplemental Figure S13: Joint interpretation for mean performance and stability for grain yield (GY) in Kg/ha (a) and mean performance (b) for grain yield (GY) analysis. The vertical dashed and solid lines show, respectively, the mean of the selected genotype and the overall mean for both mean performance and WAASB index. Black, Navy and Small Red beans are represented by bars colors black, blue and red. The x axis shows the arithmetic mean for each genotype × environment interaction and y axis the weighted average of absolute scores from the singular value decomposition of the matrix of best linear unbiased predictions for the genotype × environment interaction effects generated by a combined linear mixed-effect model (WAASB) of 72 Black, 71 Navy and 21 Small Red beans evaluated in 4 environments (BA: Bay, HU: Huron, SA: Sanilac, TU: Tuscola). The blue shaded squared in the quadrant IV are the potential selected genotypes considering combined analysis."}

knitr::include_graphics(here::here("main-figures", "S13.png"))

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r s14,echo=FALSE, warning=FALSE, message=FALSE, out.width="100%", fig.pos= "H", fig.align="center", dpi = 600, fig.cap="Supplemental Figure S14: MTME model: Genotypic correlation between GY, DPM, and PH across genotypes and locations for all market classes (BB, NB, and SR) combined analysis (CombMKT). The genotype (G) and genotype × environment interactions (GEI) correlation coefficients between the GY, DM, PH, and LD traits."}

knitr::include_graphics(here::here("main-figures", "Rplot_MTME-allmkt.png"))

```



```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```




```{r figure pred year 2, echo=FALSE, include=FALSE, fig.height=figheight, fig.width=figwidth, message=FALSE, fig.pos= "H", out.width="100%", warning=FALSE, dpi=600, fig.cap="Genetic yield gains for black (BB), navy (NB), and (SR) market classes of dry bean varieties regressed against year. For BB yields have improved by 57.7 kg ha−1 annually, while NB and SR yields have decreased by, respectively, 4.8 kg and 13.3 ha−1 annually. Best linear unbiased predictors (BLUPs) of varieties yields in each market class were based on trials in four environments (BA: Bay, HU: Huron, SA: Sanilac, TU: Tuscola) across the state of Michigan from 2017 to 2022."}

knitr::include_graphics(here::here("main-figures", "plot_year_gains-1.png"))

```


```{r echo=FALSE, message=TRUE, warning=FALSE, results='asis'}
cat("\n\n\\pagebreak\n")
```



```{r full_code, echo=F, ref.label=knitr::all_labels(),eval=FALSE}
# Full Code

# The full script of executive code contained in this document is reproduced here.
```


# Appendix C - R codes

# Getting started {.unnumbered}

The present analysis aims to dissect the genotype by environment interaction study (aka GEI) using a data set from the Dry Beans breeding program at Michigan State University - MSU.

The trait in study is the grain yield (GY) per plot (lb/plot) adjusted to the international measurements (Kg/ha). A previous data analysis (not shown here) was done to perform the historical data mining and adjust of raw data. Besides the GY, plant height (PH), date of maturity (DM) and lodging (LD) were investigate using a subset from 2021, which contains the all the data available.

The main focus of this manuscript, as describe in the [published paper](link%20here), is to investigate the varieties performance of GY across four locations at different Michigan counties. However, the MTME (three locations and 4 traits) also was studied in this manuscript when available (only in 2021 in Bay, Tuscola and Sanilac sites). Thus, different types of analysis will be performed in order to study the GEI in the Multi-Environment-Trials (MET) data to provide better varieties recommendations to the Dry Beans breeding program, which:

-   [x]  Multi-Environment Trials -- Genotype x Environment Interaction to grain yield (GY)
-   [x]  Mean performance and stability of multiple traits

## Packages

R packages version from: `r format(Sys.Date())`.

```{r message=TRUE, warning=FALSE, results='asis',echo = F}
          library(rmarkdown) 
          library(ggpmisc) 
          library(ggpp) 
          library(gridExtra) 
          library(magrittr)
          library(Matrix)
          library(mapdata) 
          library(maps) 
          library(spData) 
          library(asremlPlus) 
          library(tidyquant) 
          library(usmap) 
          library(data.table) 
          library(ggspatial)
          library(ggpattern) 
          library(viridisLite) 
          library(viridis) 
          library(flextable) 
          library(lubridate) 
          library(raster) 
          library(ggcorrplot) 
          library(ggdist) 
          library(tibble) 
          library(RColorBrewer) 
          library(rcartocolor) 
          library(metan) 
          library(ggstatsplot) 
          library(sp) 
          library(sf) 
          library(ggforce) 
          library(patchwork) 
          library(PerformanceAnalytics) 
          library(broom) 
          library(quantmod) 
          library(xts) 
          library(openxlsx) 
          library(asreml) 
          library(naniar) 
          library(TTR) 
          library(statgenGxE) 
          library(tigris) 
          library(plyr) 
          library(ggplot2) 
          library(stringr) 
          library(forcats) 
          library(tidyverse) 
          library(dplyr) 
          library(purrr) 
          library(readr) 
          library(tidyr) 
          library(cowplot) 
          library(nadiv) 
          library(DT) 
          library(zoo) 
          library(kableExtra) 


summary(report::report(sessionInfo()))

```

## Data preparation

```{r loading data, warning=FALSE, message=TRUE, results='asis'}
data_beans = read.csv("data/DataBean_MET_GYv2.csv",h=T, stringsAsFactors = T)

if (knitr::is_html_output()) {

  print_table(data_beans)
  
}else{
  
  flextable(head(data_beans)) %>%
   add_footer_lines(
     c("Varieties Dry Beans data set from 2017 to 2022",
       "Header data set showing the 6 first entry")) %>%
   autofit() %>%
   add_header_lines("Dry Beans varieties trial") %>%
  theme_design2()
}

```


```{r dataadjs1,warning=FALSE, message=TRUE, results='asis'}

# Data adjustment
# All the effect columns must be as a factor to run in ASReml-r.
cols <- c("rep", "name", "loc","year", "mkt", "year_loc")
data_beans[cols] <- lapply(data_beans[cols], factor)
data_beans <- data.table(data_beans)

```

## Descriptive Stats - Raw data {.tabset}

Data set distribution, checking data and locations of study.

```{r eval=T, child=if (fast == FALSE) '1_descriptive.Rmd'}
```



## Predicted by year mixed model analysis

To predict the BLUPs per year and plot it versus year: BLUPs prediction for the vector of the variable GY in the ith genotype, and jth loc within `year`. 

```{r eval=T, child=if (fast == FALSE) '2.1_PredByYear_MET_desc.Rmd'}
```



## Two Stage mixed model analysis

The following pipeline will be used: 1-Stage = BLUEs estimation for the vector of the variable GY in the ith genotype, and jth year within `loc`. Then the BLUEs from the 1-Stage (Yik) will be used to predict the BLUPs (Yijl) of the ith genotype in the lth location and jth rep in the 2-Stage in which this second model have `name` and `loc` effects as random.

```{r eval=T, child=if (fast == FALSE) '2_two_stage_MET_desc.Rmd'}
```

# Appendix D - R codes


## MET analysis - **Market class combined**

```{r eval=T,child=if (fast == FALSE) '2.2_comb_MKT_MET.Rmd'}
```


## MET analysis - **Black beans**

```{r eval=T,child=if (fast == FALSE) '3_BB_MET.Rmd'}
```

## MET analysis - **Navy beans**

```{r eval=T,child=if (fast == FALSE) '4_NB_MET.Rmd'}
```

## MET analysis - **Red beans**

```{r eval=T,child=if (fast == FALSE) '5_SR_MET.Rmd'}
```

# Appendix E - R codes

## Multi-trait Multi-env analysis


- GY data set: were obtained from the adjusted data in the 2nd-stage mixed models analysis.
- PH, PM, and LD were obtained from the coincident locations (TU, BA and SA) in 2021 only.

```{r eval=T,child=if (fast == FALSE) '6_MultTrait.Rmd'}
```

```{r warning=FALSE, message=TRUE, results='asis'}
cat("\n\n\\pagebreak\n")
```

# Package References

```{r echo=FALSE, message=FALSE, warning=FALSE, results='asis'}
report::cite_packages(sessionInfo())
```

<!-- # References -->