4_NB_MET.Rmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---


### MET analysis - ASReml

Running MET using `ASReml` - only to comparison of variance components with `metan` outputs

```{r mod met asreml nb, warning=FALSE, message=TRUE, results='asis'}
mod.met.asreml.nb1 <- asreml(fixed       = yield ~  loc + loc:rep,
                     random      = ~  name + name:loc,
                     data        = blues_stage.I_NB,
                     predict     = predict.asreml(classify = "name"),
                     trace       = F,
                     maxit       = 500)

summFix.nb.met.asreml<- data.frame(wald(mod.met.asreml.nb1))
summFix.nb.met.asreml

summ.nb.met.asreml<- data.frame(summary.asreml(mod.met.asreml.nb1)$varcomp)
summ.nb.met.asreml

#print(summary.asreml(mod.met.asreml.nb1)$bic)
mod.met.asreml.nb<- as.data.frame((mod.met.asreml.nb1$predictions$pvals[1:3]))
names(mod.met.asreml.nb) <- c("name", "yield_BLUPS_MET", "SE")

###
```


### MET analysis - lme4

Running MET using `metan` R package [@olivotoMeanPerformanceStability2019](Olivoto et al., 2019a).

```{r mod metan nb, warning=FALSE, message=TRUE, results='asis'}
#str(blues_stage.I_NB)
mixed_mod.nb<- 
   gamem_met(blues_stage.I_NB,
             env = loc,
             gen = name,
             rep = rep,
             resp = yield,
             random = "gen", #Default
             verbose = TRUE) #Default

```

### Printing the model outputs {.tabset}

#### Likelihood Ratio Tests

The output `LRT` contains the Likelihood Ratio Tests for genotype and genotype-vs-environment random effects. 

```{r  warning=FALSE, message=TRUE, results='asis'}
data_mod_nb_test <- get_model_data(mixed_mod.nb, "lrt")
data_mod_nb_test

#customize the display of numbers and other data in a tibble
# old <- options(pillar.sigfig = 6)
# 
# blues_stage.I_NB %>%
#   group_by(loc) %>%
#   dplyr::summarise(Mean = mean(yield, na.rm = TRUE))
```


#### Detailed parameters

```{r  warning=FALSE, message=TRUE, results='asis'}
data_mod_nb_det <- get_model_data(mixed_mod.nb, "details")
data_mod_nb_det

```



#### Random effects 

The output `LRT` contains the Likelihood Ratio Tests for genotype and genotype-vs-environment random effects. 

```{r  warning=FALSE, message=TRUE, results='asis'}
old <- options(pillar.sigfig = 8)

data_mod_nb_var <- get_model_data(mixed_mod.nb, "vcomp")
data_mod_nb_var

```


#### Variance components and genetic parameters 


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

old <- options(pillar.sigfig = 4)

data_mod_nb_comp <- get_model_data(mixed_mod.nb)
data_mod_nb_comp
```


### MET - GGE biplot {.tabset}

Genotype plus Genotype-vs-Environment interaction (GGE). Mega-environment identification in multi-environment trials (MET) according to [@W. Yan et al. 2007](link here)

#### GGE ENV biplot

GGE biplot done using: 

- **sd**: each value is divided by the standard deviation of its corresponding environment.
- **environment**: environment-centered (G+GE)
- **environment**: singular value is entirely partitioned into the environment eigenvectors, also called column metric preserving


```{r nb gge fig1, warning=FALSE, message=FALSE, fig.height=figheight, fig.width=figwidth, fig.align="center", dpi = 600}
gge_model.nb <- gge(blues_stage.I_NB, loc, name, yield, 
                    centering = "environment", #1
                    scaling = "sd", #2
                    svp = "environment")#2

a <- plot(gge_model.nb, type=4,
          size.text.env = 4.5,
          plot_theme = theme_metan(grid =  "both",color.background = transparent_color()),
         axis_expand = 1.5,
         col.alpha.circle = 0.8,
         shape.gen = NA,
         col.gen = NA,
         size.text.lab = NA,
         size.text.gen = NA,
         leg.lab=c('Env')
         #title = FALSE
         )


gge_model.nb <- gge(blues_stage.I_NB, loc, name, yield,
                    centering = "environment", #1
                    scaling = "sd", #2Y
                    svp = "environment")#2)

b <- plot(gge_model.nb, type = 6,
             size.text.env = 5,
          plot_theme = theme_metan(grid =  "both",color.background = transparent_color()),
         axis_expand = 1.5,
        # col.alpha.circle = 100,
          col.alpha.circle = 0.8,
         size.text.lab = 13
        #title = FALSE
        )

 arrange_ggplot(a, b,
                guides = "collect",
  tag_levels = "a",
  tag_prefix = "(",
  tag_suffix = ")")
 
```


#### Biplot type 3: Which-won-where

GGE biplot done using: 

- **sd**: each value is divided by the standard deviation of its corresponding environment.
- **environment**: environment-centered (G+GE)
- **genotype**: singular value is entirely partitioned into the environment eigenvectors, also called column metric preserving


```{r nb gge fig2, warning=FALSE, message=FALSE, fig.height=figheight, fig.width=figwidth, fig.align="center", dpi = 600}
gge_model.nb <- gge(blues_stage.I_NB, loc, name, yield,
                     centering = "environment", #2
                    scaling = "sd", #1
                    svp = "genotype")#2)

e <- plot(gge_model.nb, type = 3,
          size.text.env = 5,
          plot_theme = theme_metan(grid =  "both",color.background = transparent_color()),
         axis_expand = 1.2,
         size.line = 0.7,
         size.text.gen = 4,
         size.text.win = 4.5
         #title = FALSE

         )

print(e)

```


### Mean performance and stability analysis {.tabset}

WAASP index and BLUPs to estimate stability analysis.

```{r mod metan nb stab, warning=FALSE, message=FALSE, fig.height=figheight, fig.width=figwidth, fig.align="center", dpi = 600}
waasb_model_nb <- 
  waasb(blues_stage.I_NB,
        env = loc,
        gen = name,
        rep = rep,
        resp = yield,
        random = "gen", #Default
        verbose = TRUE,
        wresp = 60) #weight for response variable 60 and 40 for yielding and stability, respectively)

waasb_model<- waasb_model_nb$yield$model

waasp_plot <- plot_scores(waasb_model_nb, type = 3,
          title = FALSE,
          size.tex.gen = 4,
          size.tex.env = 4,
          size.tex.lab = 13,
         # highlight = c("N38", "N6" , "N61", "N35" ,"N52", "N22"),
         plot_theme = theme_metan(grid =  "both",color.background = transparent_color())
        ) +
  
  geom_mark_rect(aes(filter =  Code  %in% c("N70", "N37", "N69", "N60"),
                     ),
               label.fontsize = 10,
               show.legend = F,
               con.cap = 0,
               con.colour = "red",
               color = "red",
               expand = 0.005,
               label.buffer = unit(10, "cm"))+
#theme_gray()+
theme(legend.position = c(0.1, 0.9),
      legend.background = element_blank(),
      legend.title = element_blank(),
      aspect.ratio = 1) +
  labs(x = "GY")
print(waasp_plot)

```

```{r nb selected waasb, warning=FALSE, message=TRUE, results='asis'}

waasb_model_meanWaasb<-mean(waasb_model$WAASB)
waasb_model_meanY<-mean(waasb_model$Y)

selected <- waasb_model %>%
  dplyr::filter(Y >= waasb_model_meanY & WAASB <= waasb_model_meanWaasb) 

selected_table <- selected 

if (knitr::is_html_output()) {
  
  print_table(selected_table)
  
}else{
  
selected_table[,1:8]
}



#selected$Code
 
```



#### Selection differentials

```{r mod metan nb stab2, warning=FALSE, message=FALSE, fig.height=figheight, fig.width=figwidth, fig.align="center", dpi = 600, fig.cap="Mean performance (a) and stability (b) for grain yield (GY) of 71 Navy beans genotypes. The vertical dashed and solid lines shows, respectivelly, the  mean of the selected genotype and the overall mean for both mean performance and WAASB index"}

#Create a data frame with BLUPS - selected and non-selected
blups_sel <-
  gmd(waasb_model_nb, "blupge") %>%
  add_cols(SELECTED = ifelse(GEN %in% selected$Code, "yes", "no")) %>% 
    dplyr::rename(BLUPs_sel = yield) %>% 
  droplevels()

blups_sel_mean<-
  gmd(waasb_model_nb, "blupge") %>%
  add_cols(SELECTED = ifelse(GEN %in% selected$Code, "yes", "no")) %>% 
  filter(SELECTED == "yes") %>% 
  dplyr::summarise(mean_GY = mean(yield,na.rm = TRUE), n = n()) 

# Create a data frame with the waasb index - selected and non-selected
waasb_sel <-
  gmd(waasb_model_nb, "WAASB") %>%
  add_cols(SELECTED = ifelse(GEN %in% selected$Code, "yes", "no")) %>% 
  dplyr::rename(WAASB_sel = yield) %>% 
  droplevels()
#str(waasb_sel)

waasb_sel_mean<-
  gmd(waasb_model_nb, "WAASB") %>%
  add_cols(SELECTED = ifelse(GEN %in% selected$Code, "yes", "no")) %>% 
  filter(SELECTED == "yes") %>% 
  dplyr::summarise(mean_GY = mean(yield,na.rm = TRUE), n = n()) 

p1<- plot_selected(blups_sel, GEN, BLUPs_sel, mean_sel = blups_sel_mean$mean_GY) +
  labs(y = "GY")

p3<- plot_selected(waasb_sel, GEN, WAASB_sel, mean_sel = waasb_sel_mean$mean_GY) +
  labs(y = "WAASB index")

arrange_ggplot(p1, p3,
  guides = "collect",
  tag_levels = "a",
  tag_prefix = "(",
  tag_suffix = ")")
```

Percentage (SD_gain in %) gain from the selected genotypes compared to the general mean.


```{r mod metan nb stab3_gain, warning=FALSE, message=TRUE, results='asis'}

blups_sel2 <-
  gmd(waasb_model_nb, "blupg") %>%
  add_cols(SELECTED = ifelse(GEN %in% selected$Code, "yes", "no")) %>% 
    dplyr::rename(BLUPs_sel = yield) %>% 
  droplevels()

blups_sel_mean2<-
  gmd(waasb_model_nb, "blupg") %>%
  add_cols(SELECTED = ifelse(GEN %in% selected$Code, "yes", "no")) %>% 
  filter(SELECTED == "yes") %>% 
  dplyr::summarise(mean_GY = mean(yield,na.rm = TRUE), n = n()) 

SD_blups<- as_tibble((blups_sel_mean2$mean_GY/mean(blups_sel2$BLUPs_sel, na.rm = T)) -1)*100
SD_WAASP<- as_tibble((waasb_sel_mean$mean_GY /mean(waasb_sel$WAASB_sel, na.rm = T)) -1)*100

SD_comb<- full_join(SD_blups, SD_WAASP, by = "value") %>% 
  dplyr::rename(SD_gain = value) %>% 
  tibble::add_column(Comp_name = c('BLUPs', 'WAASB')) %>% 
  relocate(Comp_name)

SD_comb$n_selected<- blups_sel_mean2$n
SD_comb
```


```{r mod metan nb stab4, warning=FALSE, message=TRUE, results='asis'}

blups_sel2$mean_blup <- mean(blups_sel2$BLUPs_sel, na.rm = T)
waasb_sel$mean_waasb <- mean(waasb_sel$WAASB_sel, na.rm = T)

#str(waasb_sel)
data_comb<- merge(blups_sel2, waasb_sel, by = c("GEN", "SELECTED"))
#names(data_comb)
## SD for each genotype
data_sel_perc <- data_comb %>%
 rowwise %>%
  mutate(Perc_blup_gain = ((BLUPs_sel/mean_blup)*100)-100) %>% 
  mutate(Perc_WAASB_gain = ((WAASB_sel/mean_waasb)*100)-100) %>% 
  as_tibble()

# data_sel_perc_mean <- data_sel_perc %>% 
#   dplyr::filter(SELECTED  == "yes")
# 
# mean(data_sel_perc_mean$Perc_blup_gain)

if (knitr::is_html_output()) {
  
  print_table(data_sel_perc)
  
}else{
  
data_sel_perc[,1:7]
}


data_sel_perc<- data_sel_perc %>% 
  dplyr::relocate(GEN,SELECTED,BLUPs_sel,mean_blup,Perc_blup_gain,
                 WAASB_sel,mean_waasb ,Perc_WAASB_gain)

#write.xlsx(data_sel_perc, "./data/sel_SD_nb_2.xlsx")


data_sel_perc2 <- data_sel_perc %>% 
  dplyr::select(GEN,SELECTED, BLUPs_sel, WAASB_sel, Perc_blup_gain, Perc_WAASB_gain)

data_sel_perc2

```

```{r mod metan nb stab5, warning=FALSE, message=FALSE, fig.height=figheight, fig.width=figwidth, fig.align="center", dpi = 600}
##BLUPs indexes
stab_blups_nb<- blup_indexes(waasb_model_nb) 
stab_blups_nb<- as_tibble(stab_blups_nb$yield)

data_waasby <- waasb_model_nb$yield$model %>% 
  dplyr::filter(type != "ENV") %>% 
  dplyr::select("Code", "WAASBY", "OrWAASBY") %>% 
  dplyr::rename(GEN = Code)

stab_blups_nb<- stab_blups_nb %>% 
  full_join(data_waasby, by = "GEN")

if (knitr::is_html_output()) {
  
  print_table(stab_blups_nb)
  
}else{
  
stab_blups_nb[,1:8]
}



# library(openxlsx)
# write.xlsx(stab_blups_nb, "./data/blups_nb_2.xlsx")

```

#### Scenarios of waasby estimation


Planning different scenarios of waasby estimation by changing the weights assigned to the stability and the mean performance according to [@olivotoMeanPerformanceStability2019](Olivoto et al., 2019a).

```{r mod metan nb stab6, warning=FALSE, message=FALSE, fig.height=figheight, fig.width=figwidth, fig.align="center", dpi = 600}
scenarios <- wsmp(waasb_model_nb,progbar = F)

scen1<- plot(scenarios, type = 1) + 
  theme(axis.text.y = element_text(size=8))

scen2<-plot(scenarios, type = 2) +
  theme(axis.text.y = element_text(size=8))

arrange_ggplot(scen1, scen2,
  guides = "collect",
  legend.position = "right",
  tag_levels = "a",
  tag_prefix = "(",
  tag_suffix = ")")
```


#### Coincidence index of genotype selection

Computes the coincidence index (Hamblin and Zimmermann, 1986) as follows:

$$
{CI = \frac{A-C}{M-C}\times 100}
$$

where *A* is the number of selected genotypes common to different methods; 
*C* is the number of expected genotypes selected by chance; 
and *M* is the number of genotypes selected according to the selection intensity.



```{r mod metan nb stab7, warning=FALSE, message=FALSE, fig.height=figheight, fig.width=figwidth, fig.align="center", dpi = 600}
coinc_1 <- stab_blups_nb %>% dplyr::select(GEN,HMRPGV_R) %>% arrange(HMRPGV_R)
coinc_2 <- stab_blups_nb %>% dplyr::select(GEN,RPGV_R) %>% arrange(RPGV_R)
coinc_3 <- stab_blups_nb %>% dplyr::select(GEN,HMGV_R) %>% arrange(HMGV_R)
coinc_4 <- stab_blups_nb %>% dplyr::select(GEN,OrWAASBY) %>% arrange(OrWAASBY)
coinc_5 <- stab_blups_nb %>% dplyr::select(GEN,WAASB_R) %>% arrange(WAASB_R)

selc_perc<- round(nrow(stab_blups_nb)*0.2)

coinc_1.1 <-1
coinc_1.2 <- coincidence_index(sel1 = coinc_1$GEN[1:selc_perc], 
                                        sel2 = coinc_2$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_1.3 <- coincidence_index(sel1 = coinc_1$GEN[1:selc_perc], 
                                        sel2 = coinc_3$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_1.4 <- coincidence_index(sel1 = coinc_1$GEN[1:selc_perc], 
                                        sel2 = coinc_4$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_1.5 <- coincidence_index(sel1 = coinc_1$GEN[1:selc_perc], 
                                        sel2 = coinc_5$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_2.2 <-1
coinc_2.3 <- coincidence_index(sel1 = coinc_2$GEN[1:selc_perc], 
                                        sel2 = coinc_3$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_2.4 <- coincidence_index(sel1 = coinc_2$GEN[1:selc_perc], 
                                        sel2 = coinc_4$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_2.5 <- coincidence_index(sel1 = coinc_2$GEN[1:selc_perc], 
                                        sel2 = coinc_5$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_3.3<- 1
coinc_3.4 <- coincidence_index(sel1 = coinc_3$GEN[1:selc_perc], 
                                        sel2 = coinc_4$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_3.5 <- coincidence_index(sel1 = coinc_3$GEN[1:selc_perc], 
                                        sel2 = coinc_5$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_4.4 <- 1
coinc_4.5 <- coincidence_index(sel1 = coinc_4$GEN[1:selc_perc], 
                                        sel2 = coinc_5$GEN[1:selc_perc], 
                                        total = 71)/100
coinc_5.5 <- 1


coinc<- c(coinc_1.1,coinc_1.2,coinc_2.2,coinc_1.3,coinc_2.3,
          coinc_3.3,coinc_1.4, coinc_2.4, coinc_3.4,
          coinc_4.4, coinc_1.5, coinc_2.5,
          coinc_3.5, coinc_4.5,
          coinc_5.5)
  
z=matrix(0,5,5)
z[upper.tri(z)| row(z)==col(z)] <- coinc

rownames(z)=c(
"HMRPGV",
"RPGV",
'HMGV',
'WAASBY',
'WAASB')

colnames(z)=rownames(z)

plotNB<- ggcorrplot(z, colors = c("#6D9EC1", "gray" ,"#E46726"),  
           show.legend = T,
legend.title = "CI" ,lab_size=5,tl.srt = 90,type = c("upper"), lab = T,digits = 4,
outline.color = "white",pch.col = "white", tl.col = "blue",show.diag = F) +
  labs(title = "BLUP-based stability indexes coincidence at NB",
           subtitle = "Selection intensity of 20% top genotypes")

print(plotNB)

```