---
title: "Empirical inputs for simulations - IITA"
site: workflowr::wflow_site
date: "2021-July-01"
output: 
  workflowr::wflow_html:
    code_folding: hide
editor_options:
  chunk_output_type: inline
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = F, 
                      tidy='styler', tidy.opts=list(strict=FALSE,width.cutoff=100), highlight=TRUE)
```

```{bash set-up remote R environment, eval=F}
# pull a singularity image and save in the file rocker.sif
# singularity pull ~/rocker2.sif docker://rocker/tidyverse;

# 1) start a screen shell 
screen; # or screen -r if re-attaching...
# 2) start the singularity Linux shell inside that
singularity shell ~/rocker2.sif; 
#singularity shell /workdir/$USER/rocker.sif; 
# Project directory, so R will use as working dir.
cd /home/mw489/implementGMSinCassava/;
# 3) Start R
R
```

# Concept

We want to simulate the advantages of different breeding pipelines
(variety development pipelines; VDPs) based on proposed differences in
plot-size (number of plants), number of reps, locations and overall
trial-stage size.

**Question:** How do we empirically know what difference varying these
parameters is associated with?

We can't compute the selection index *per se* for each trial because of
variation in the traits scored. Then again, neither can breeders.

We want:

-   an objective approach to estimate the selection error for an SI on a
    trial-by-trail basis,
-   ideally using existing GS pipeline / data,
-   dealing with heterogeneity in traits-observed in given trial.

Proposed approach:

-   Generate a best estimate of "true" net merit using the *SELECTION
    INDEX* GEBV/GETGV from full genomic prediction, *i.e.* the entire
    available training population and all of the latest available data.
    The predictions made
    [here](https://wolfemd.github.io/implementGMSinCassava/06-GenomicPredictions.html)
    are the most up-to-date.

    -   Alternative best estimate of "true" net merit: actual national
        performance trial data (e.g. from NCRPs in Nigeria?)

-   For each trial, and each trait scored:

    -   Fit mixed-model

    -   Extract trial-specific BLUPs

    -   Compute the SELIND for the current trial using BLUPs for
        whatever component traits were scored ($SI_{TrialBLUP}$).

    -   Regress $SI_{GETGV}$ on the $SI_{TrialBLUP}$

        -   Extract the $\sigma^2_e$ of the regression as the
            trial-specific estimate of the selection error

-   Finally, the relationship between SI error and plot-size can be
    plotted and a function for SI error given plot-size can be developed
    that can be plugged into simulations of any VDP of interest.

# Selection Index GETGV: proxy for true net merit

-   Generate a best estimate of "true" net merit using the *SELECTION
    INDEX* GETGV from full genomic prediction, *i.e.* the entire
    available training population and all of the latest available data.
    The predictions made
    [here](https://wolfemd.github.io/implementGMSinCassava/06-GenomicPredictions.html)
    are the most up-to-date.


```{r load getgv and support data}
library(tidyverse); library(magrittr); 
# GBLUPs
### Two models AD and DirDom
#gpreds_ad<-readRDS(file = here::here("output","genomicPredictions_ModelAD.rds"))
gpreds_dirdom<-readRDS(file = here::here("output","genomicPredictions_ModelDirDom.rds"))
si_getgvs<-gpreds_dirdom$gblups[[1]] %>% 
#si_getgvs<-gpreds_ad$gblups[[1]] %>% 
  filter(predOf=="GETGV") %>% 
  select(GID,SELIND)

## IITA Germplasm Ages
ggcycletime<-readxl::read_xls(here::here("data","PedigreeGeneticGainCycleTime_aafolabi_01122020.xls")) %>% 
  mutate(Year_Accession=as.numeric(Year_Accession))

# Need germplasmName field from raw trial data to match GEBV and cycle time
dbdata<-readRDS(here::here("output","IITA_ExptDesignsDetected_2021May10.rds"))
si_getgvs %<>% 
  left_join(dbdata %>% 
              distinct(germplasmName,GID)) %>% 
  group_by(GID) %>% 
  slice(1) %>% 
  ungroup()
# table(ggcycletime$Accession %in% si_getgvs$germplasmName)
# FALSE  TRUE 
#   193   614 

si_getgvs %<>% 
  left_join(.,ggcycletime %>% 
              rename(germplasmName=Accession) %>% 
              mutate(Year_Accession=as.numeric(Year_Accession))) %>% 
  mutate(Year_Accession=case_when(grepl("2013_|TMS13",germplasmName)~2013,
                                  grepl("TMS14",germplasmName)~2014,
                                  grepl("TMS15",germplasmName)~2015,
                                  grepl("TMS18",germplasmName)~2018,
                                  !grepl("2013_|TMS13|TMS14|TMS15|TMS18",germplasmName)~Year_Accession))

# Declare the "eras" as PreGS\<2012 and GS\>=2013.
si_getgvs %<>% 
  filter(Year_Accession>2012 | Year_Accession<2005)
si_getgvs %<>% 
  mutate(GeneticGroup=ifelse(Year_Accession>=2013,"GS","PreGS"))

# SELECTION INDEX WEIGHTS
## from IYR+IK
## note that not ALL predicted traits are on index
SIwts<-c(logFYLD=20,
         HI=10,
         DM=15,
         MCMDS=-10,
         logRTNO=12,
         logDYLD=20,
         logTOPYLD=15,
         PLTHT=10) 

# SELIND GETGVS (for input to estimateSelectionError func below)
getgvs<-si_getgvs %>% 
  select(GID,SELIND); 
```

```{r show rate of genetic gain - getgv vs. accession age}
# , fig.height=10, fig.width=12
si_getgvs %>% 
  select(GeneticGroup,GID,Year_Accession,SELIND) %>% 
  ggplot(.,aes(x=Year_Accession,y=SELIND,color=GeneticGroup)) + 
  geom_point(size=1.25) + geom_smooth(method=lm, se=TRUE, size=1.5) + 
  # facet_wrap(~Trait,scales='free_y', ncol=2) + 
  theme_bw() +
  theme(axis.text = element_text(face = 'bold',angle = 0, size=14),
        axis.title = element_text(face = 'bold',size=16),
        strip.background.x = element_blank(),
        strip.text = element_text(face='bold',size=18)) + 
  scale_color_viridis_d() + 
  labs(title = "Selection Index GETGV vs. Accession Age by 'era' [GS vs. PreGS]",
       subtitle = "SI GETGV from modelType='DirDom'")
```

```{r barplot of mean GEBV vs. Cycle}
si_gebvs<-gpreds_dirdom$gblups[[1]] %>% 
  filter(predOf=="GEBV") %>% 
  select(GID,SELIND) %>% 
  mutate(GeneticGroup=case_when(grepl("2013_|TMS13",GID)~"C1",
                                  grepl("TMS14",GID)~"C2",
                                  grepl("TMS15",GID)~"C3",
                                  grepl("TMS18",GID)~"C4",
                                  !grepl("2013_|TMS13|TMS14|TMS15|TMS18",GID)~"PreGS"),
         GeneticGroup=factor(GeneticGroup,levels = c("PreGS","C1","C2","C3","C4")))
si_gebvs %>% 
  group_by(GeneticGroup) %>% 
  summarize(meanGEBV=mean(SELIND),
            stdErr=sd(SELIND)/sqrt(n()),
            upperSE=meanGEBV+stdErr,
            lowerSE=meanGEBV-stdErr) %>% 
  ggplot(.,aes(x=GeneticGroup,y=meanGEBV,fill=GeneticGroup)) + 
  geom_bar(stat = 'identity', color='gray60', size=1.25) + 
  geom_linerange(aes(ymax=upperSE,
                     ymin=lowerSE), color='gray60', size=1.25) + 
  #facet_wrap(~Trait,scales='free_y', ncol=1) + 
  theme_bw() +
  geom_hline(yintercept = 0, size=1.15, color='black') + 
  theme(axis.text.x = element_text(face = 'bold',angle = 0, size=12),
        axis.title.y = element_text(face = 'bold',size=14),
        legend.position = 'none',
        strip.background.x = element_blank(),
        strip.text = element_text(face='bold',size=14)) + 
  scale_fill_viridis_d() + 
  labs(x=NULL,y="Mean GEBVs")
```

# Estimate selection error for each trial

-   For each trial, and each trait scored:

    -   Fit mixed-model

    -   Extract trial-specific BLUPs

    -   Compute the SELIND for the current trial using BLUPs for
        whatever component traits were scored ($SI_{TrialBLUP}$).

    -   Regress $SI_{GETGV}$ on the $SI_{TrialBLUP}$

        -   Extract the $\sigma^2_e$ of the regression as the
            trial-specific estimate of the selection error

## Available plot-level data

See [here](01-cleanTPdata.html) for details about the cleaned trial data most recently downloaded from IITA/Cassavabase and used below.

Quick summary of the number of unique plots, locations, years, etc. in the cleaned plot-basis data.

```{r summarize trial data}
dbdata %>% 
  summarise(Nplots=nrow(.),
            across(c(locationName,studyYear,studyName,TrialType,GID), ~length(unique(.)),.names = "N_{.col}")) %>% 
  rmarkdown::paged_table()
```

```{r}
dbdata %>% 
  count(TrialType,CompleteBlocks,IncompleteBlocks) %>% 
  spread(TrialType,n) %>% 
  rmarkdown::paged_table()
```

Not much NCRP data, but there is some.

```{r plot MaxNOHAV vs. TrialType, cols.print=16}
dbdata %>% 
  filter(studyYear>=2013) %>% 
  distinct(studyYear,locationName,studyName,TrialType,CompleteBlocks,IncompleteBlocks,MaxNOHAV) %>% 
  filter(!is.na(MaxNOHAV)) %>% #count(TrialType)
  mutate(TrialType=factor(TrialType,levels=c("CrossingBlock","GeneticGain","CET","ExpCET","PYT","AYT","UYT","NCRP"))) %>% 
  ggplot(.,aes(x=TrialType,y=MaxNOHAV,fill=TrialType)) + 
  geom_boxplot(notch = T) + 
  theme_bw() + theme(axis.text.x = element_text(angle=45,vjust=.5)) +
  labs(title = "Max number harvested as a proxy for planned plot size",
       subtitle="MaxNOHAV = The maximum number stands harvested per trial")
```

## Trial-by-trial mixed-models

Decided it makes sense to restrict consideration to >2012 to measure the selection error during the current "era" at IITA.

```{r subset to 2013 trial data}
trials2keep<-dbdata %>% 
  filter(studyYear>=2013) %>% 
  distinct(studyYear,locationName,studyName,TrialType,CompleteBlocks,IncompleteBlocks,MaxNOHAV) %>% 
  filter(!is.na(MaxNOHAV)) %$% 
  unique(studyName)
# 633 trials, 94830 plots

trialdata<-dbdata %>% 
  filter(studyYear>=2013,
         studyName %in% trials2keep) %>% 
  nest(TrialData=-c(studyYear,locationName,studyName,TrialType,CompleteBlocks,IncompleteBlocks,MaxNOHAV))

trialdata %<>% 
  mutate(propGenotyped=map_dbl(TrialData,
                              ~length(which(!is.na(unique(.$FullSampleName))))/length(unique(.$GID))))

trialdata %>% head
```
```{r unit test inputs for estimateSelectionError, eval=F}
# TrialData<-trialdata$TrialData[[1]]
# CompleteBlocks<-trialdata$CompleteBlocks[[1]]
# IncompleteBlocks<-trialdata$IncompleteBlocks[[1]]
# 
# TrialData<-trialdata %>% filter(propGenotyped>0.75) %>% slice(4) %$% TrialData[[1]]
# CompleteBlocks<-trialdata %>% filter(propGenotyped>0.75) %>% slice(4) %$% CompleteBlocks[[1]]
# IncompleteBlocks<-trialdata %>% filter(propGenotyped>0.75) %>% slice(4) %$% IncompleteBlocks[[1]]
# ncores=4

# rm(TrialData,CompleteBlocks,IncompleteBlocks)
```

Function to estimation selection error for each trial: `estimateSelectionError()`.
```{r estimateSelectionError function, eval=F}
estimateSelectionError<-possibly(function(TrialData,CompleteBlocks,IncompleteBlocks,
                                 SIwts,getgvs,...){
  
  # SET-UP THE DATA TRAIT-BY-TRIAL~~~~~~~~~~~~~~~~~~~~~
  blups<-TrialData %>% 
    select(observationUnitDbId,GID,plotNumber,repInTrial,blockInRep,PropNOHAV,
           all_of(names(SIwts))) %>% 
    pivot_longer(cols = c(all_of(names(SIwts))),
                 names_to = "Trait",
                 values_to = "TraitValue") %>% 
    nest(TraitByTrialData=c(-Trait))
  
  # FIT MIXED-MODELS TRAIT-BY-TRIAL~~~~~~~~~~~~~~~~~~~~~~
  
  ## if model fails, by design, returns NULL for a given trait or trait-trial
  ## Output will simply be absent. 
  fit_model<-possibly(function(Trait,TraitByTrialData,CompleteBlocks,IncompleteBlocks){
    # debug settings ~~~~~~~~~~
    # TraitByTrialData<-blups$TraitByTrialData[[1]]
    # Trait<-blups$Trait[[1]]
    # TraitByTrialData<-blups$TraitByTrialData[[8]]
    # Trait<-blups$Trait[[8]]
    # rm(TraitByTrialData,Trait)
    
    # Model formula based on trial design
    modelFormula<-paste0("TraitValue ~ (1|GID)")
    modelFormula<-ifelse(CompleteBlocks,
                         paste0(modelFormula,"+(1|repInTrial)"),modelFormula)
    modelFormula<-ifelse(IncompleteBlocks,
                         paste0(modelFormula,"+(1|blockInRep)"),modelFormula)
    modelFormula<-ifelse(Trait %in% c("logRTNO","logFYLD","logTOPYLD","logDYLD"),
                         paste0(modelFormula,"+PropNOHAV"),modelFormula)
    require(lme4); require(lme4)
    model_out<-lmer(as.formula(modelFormula),data=TraitByTrialData)
    propMiss<-length(which(is.na(TraitByTrialData$TraitValue))) / length(TraitByTrialData$TraitValue)
    
    if(is.na(model_out) | !exists("model_out")){
      out <-tibble(H2=NA,VarComps=list(NULL),BLUPs=list(NULL),Model=modelFormula,propMiss=propMiss)
    } else {
      varcomps<-as.data.frame(VarCorr(model_out))[,c("grp","vcov")] %>%
        spread(grp,vcov)
      Vg<-varcomps$GID
      H2<-Vg/(Vg+varcomps$Residual)
      BLUP<-ranef(model_out, condVar=TRUE)[["GID"]]
      PEV <- c(attr(BLUP, "postVar"))
      blup<-tibble(GID=rownames(BLUP),BLUP=BLUP$`(Intercept)`,PEV=PEV) %>%
        mutate(REL=1-(PEV/Vg),
               drgBLUP=BLUP/REL,
               WT=(1-H2)/((0.1 + (1-REL)/REL)*H2))
      out <- tibble(H2=H2,
                    VarComps=list(varcomps),
                    BLUPs=list(blup),
                    Model=modelFormula,
                    propMiss=propMiss) }
    return(out) },
    otherwise = NULL)
  
  blups %<>%
    mutate(modelOut=pmap(.,fit_model,
                         CompleteBlocks=CompleteBlocks,
                         IncompleteBlocks=IncompleteBlocks)) %>% 
    select(-TraitByTrialData) %>% 
    unnest(modelOut) %>% 
    unnest(VarComps)
  
  # COMPUTE SELIND FROM BLUPs~~~~~~~~~~~~~~~~~~~~~~
  si_blups<-blups %>% 
    select(Trait,BLUPs) %>% 
    unnest(BLUPs) %>% 
    select(GID,Trait,BLUP) %>% 
    spread(Trait,BLUP) %>% 
    select(GID,any_of(names(SIwts))) %>% 
    column_to_rownames(var = "GID") %>% 
    as.matrix
  si_blups<-si_blups%*%SIwts[colnames(si_blups)] %>% 
    as.data.frame %>% 
    rownames_to_column(var = "GID") %>% 
    rename(SI_BLUP=V1) %>% 
    left_join(getgvs)
  
  # Correlation between SELIND (SI GETGV) and the SI of BLUPs for current trial
  cor2si<-si_blups %$% cor(SI_BLUP,SELIND,use = 'complete.obs')
  
  # Regress SELIND on SI_BLUPs
  regSIonTrialBLUP<-lm(SELIND~SI_BLUP,data = si_blups)
  
  # TWO MEASURES OF SELECTION ERROR
  ## 1) regression r-squared [or 1 - r2_si actually]
  ## 2) mean squared error, 
  ##### not sure if the scaling will make sense across trials 
  ##### because of differences in traits included
  r2_si<-regSIonTrialBLUP %>% summary %$% r.squared
  
  mse<-regSIonTrialBLUP %>% anova %>% as.data.frame %>% .["Residuals","Mean Sq"]
  
  NcloneForReg<-si_blups %>% filter(!is.na(SI_BLUP),
                                    !is.na(SELIND)) %>% nrow()
  # Collect and return outputs for current trial
  trial_out<-tibble(cor2si=cor2si,
                    r2_si=r2_si,
                    TrialMSE=mse,
                    NcloneForReg=NcloneForReg,
                    SI_BLUPs=list(si_blups),
                    BLUPs=list(blups))

  return(trial_out) 
},otherwise = NA)
```

Run `estimateSelectionError()` across 633 trials. 

```{r, eval=F}
#require(furrr); require(future.callr); plan(callr, workers = 10)
require(furrr); plan(multicore, workers = 20)
options(future.globals.onReference = "error")
options(future.globals.maxSize=+Inf); options(future.rng.onMisuse="ignore")

trialdata %<>% 
  mutate(SelectionError=future_pmap(.,estimateSelectionError,
                                    SIwts=SIwts,getgvs=getgvs))

saveRDS(trialdata,here::here("output","estimateSelectionError.rds"))
```

# RESULTS


## Plot selection error vs. MaxNOHAV

Out of 633 trials, 472 produced successful model fits and subsequent estimates of TrialMSE (selection index error). 
```{r}
estSelError<-readRDS(here::here("output","estimateSelectionError.rds"))
estSelError %<>% 
  select(-TrialData) %>% 
  unnest(SelectionError) %>% 
  select(-SI_BLUPs,-BLUPs,-SelectionError) %>% 
  filter(!is.na(TrialMSE))
```

Here's the `str()` of the estimates I made.
```{r}
estSelError %>% str
```

  - cor2si = correlation between SI computed from each trial's BLUPs and the SI computed from GETGV (all training data and traits used)
  - r2_si = r-squared, regression of SI_GETGV on SI_TrialBLUP
  - TrialMSE = mean squared error from that regression
  - NcloneForReg = the number of clones with estimates of both SI_TrialBLUP and SI_GETGV for a given trial. Avoid considering trials with too few data points.
  
  
```{r}
estSelError %>% rmarkdown::paged_table()
```


```{r cor2si vs MaxNOHAV}
estSelError %>% 
  ggplot(.,aes(x=MaxNOHAV,y=cor2si,color=NcloneForReg,size=NcloneForReg)) + 
  geom_point() + geom_smooth() + theme_bw() + theme(panel.grid = element_blank()) +
  labs(title = "The cor(SI_GETGV,SI_TrialBLUPs) vs. MaxNOHAV as a proxy for plot size")
```

```{r TrialMSE vs MaxNOHAV - NcloneForReg more than 5}
estSelError %>% 
  filter(NcloneForReg>=5) %>% 
  ggplot(.,aes(x=MaxNOHAV,y=TrialMSE,color=NcloneForReg,size=NcloneForReg)) + 
  geom_point() + geom_smooth() + theme_bw() + theme(panel.grid = element_blank()) +
  labs(title = "Trial-specific selection error ests. vs. MaxNOHAV",
       subtitle = "NcloneForReg>=5")
```

```{r TrialMSE vs MaxNOHAV - pos cor only}
estSelError %>% 
  filter(NcloneForReg>=10,
         cor2si>=0) %>% 
  ggplot(.,aes(x=MaxNOHAV,y=TrialMSE)) + 
  geom_point() + geom_smooth() + theme_bw() + theme(panel.grid = element_blank()) +
  labs(title = "Trial-specific selection error ests. vs. MaxNOHAV",
       subtitle = "Trials with Positive 'cor2si'. NcloneForReg>=10 and cor2si>=0")
```

```{r TrialMSE vs MaxNOHAV - binner - Nclone more than 5}
estSelError %>% 
  filter(NcloneForReg>=5) %>% 
  ggplot(.,aes(x=cut_width(MaxNOHAV,5),y=TrialMSE)) + 
  geom_boxplot() + theme_bw() + theme(panel.grid = element_blank()) +
  labs(title = "Trial-specific selection error ests. vs. MaxNOHAV",
       subtitle = "NcloneForReg>=5",
       xlab = "Binned MaxNOHAV")
```

```{r TrialMSE vs MaxNOHAV - big and pos cor only}
estSelError %>% 
  filter(NcloneForReg>=10,
         cor2si>=0.2) %>% 
  ggplot(.,aes(x=cut_width(MaxNOHAV,5),y=TrialMSE)) + 
  geom_boxplot() + theme_bw() + theme(panel.grid = element_blank()) +
  labs(title = "Trial-specific selection error ests. vs. MaxNOHAV",
       subtitle = "NcloneForReg>=10 and cor2si>=0.2",
       xlab = "Binned MaxNOHAV")
```

```{r 1-r2 vs MaxNOHAV}
estSelError %>% 
  filter(NcloneForReg>=5) %>% 
  ggplot(.,aes(x=MaxNOHAV,y=1-r2_si,size=NcloneForReg,color=NcloneForReg)) + 
  geom_point() + geom_smooth(method = 'lm') + geom_smooth(color='red')
```

```{r cor2si vs r2_si}
estSelError %>% 
  filter(NcloneForReg>5) %>% 
  ggplot(.,aes(x=r2_si,y=cor2si,size=NcloneForReg,color=NcloneForReg)) + 
  geom_point()
```

```{r TrialMSE vs cor2si}
estSelError %>% 
  filter(NcloneForReg>=5) %>% 
  ggplot(.,aes(x=cor2si,y=TrialMSE,size=NcloneForReg,color=NcloneForReg)) + 
  geom_point()
```

## Regression analysis

Try to measure an effect size of increasing the number of stands per plot.

Regression 1: NcloneForReg>=5
```{r}
lm(TrialMSE~MaxNOHAV+NcloneForReg,
   data = estSelError %>% filter(NcloneForReg>=5)) %>% summary
```

Regression 2: NcloneForReg>=10, cor2si>=0.2
```{r}
lm(TrialMSE~MaxNOHAV+NcloneForReg,
   data = estSelError %>% filter(NcloneForReg>=10,
                                 cor2si>=0.2)) %>% summary
```


