Code_for_Release.Rmd

---
title: "Code for Release"
author: "Zach Winn"
date: "`r Sys.Date()`"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r, include=FALSE}
#global options
knitr::opts_chunk$set(echo = TRUE)

#library
library(tidyverse) #data organization
library(readxl) #reading in data
library(gaston) #reading in genomic data
library(rrBLUP) #calculations
library(asreml) #gBLUP models
library(caret) #phenomic models
library(knitr) #display tables
library(fuzzyjoin) #match names
library(Hmisc) #calculate correlations
library(reshape2) #make new data frames
library(ggcorrplot) #visualize correlations
library(parallel) #parallel computing
library(foreach) #parallel computing
library(iterators) #parallel computing
library(doParallel) #parallel computing
library(stringr) #string managment 

#set wd
setwd("C:/Users/zwinn/OneDrive/Publication/Manuscripts In Progress/PS in CSU Material/Scripts/Scripts_for_Release")
```

# Read in data

Here, I read in the relevant data and make sure everything is formatted correctly...I am also coding the data to maintain the privacy of the CSU program. 

```{r}
####coding data####
#read in data
files<-list.files()

if(length(unique(c("pheno_coded.txt.gz","geno_coded.txt.gz","geno_amat_coded.txt.gz") %in% files))==1 &
   unique(c("pheno_coded.txt.gz","geno_coded.txt.gz","geno_amat_coded.txt.gz") %in% files)[1]==TRUE){
  
  #print
  print("Data has already been privatized, reading in privatized data")
  
  #read in files
  yp_dat<-read.table(gzfile("pheno_coded.txt.gz"),
                     header = TRUE)
  gbs_dat<-read.table(gzfile("geno_coded.txt.gz"),
                      header = TRUE)
  gbs_amat<-read.table(gzfile("geno_amat_coded.txt.gz"),
                       header = TRUE)
  
}else{
  
  #print
  print("Private data is not found. Attempting to convert data.")
  
  #read in yield plot data
  yp_dat<-read.csv("agronomic_drone_data_aligned.csv")
  
  #read in headrow data
  hr_dat_20<-read.csv("headrow_formatted_file_2020.csv")
  hr_dat_21<-read.csv("headrow_formatted_file_2021.csv")
  
  #format data
  yp_dat<-yp_dat %>%
    select(-thermal_mean_flight_1,
           -thermal_mean_flight_2,
           -thermal_mean_flight_3,
           -thermal_mean_flight_4) %>%
    mutate(grain_yield=grain_yield/0.0149) %>%
    filter(!nursery=="cl_pyn")
  
  #read in gbs data
  #gbs_dat<-as.matrix(read.vcf("POPVAR_CALL_2022_filt_imp.vcf.gz", convert.chr = F))-1
  
  #select only individuals in the data
  #gbs_dat<-gbs_dat[rownames(gbs_dat) %in% unique(yp_dat$fullsamplename),]
  
  #calculate relationship matrix (this will take a bit...)
  #gbs_amat<-as.matrix(A.mat(gbs_dat))
  
  #try troubleshooting with another GBS matrix
  gbs_dat<-read.vcf("POPVAR_CALL_2022_filt_imp.vcf.gz", convert.chr = F)
  gbs_dat<-gbs_dat[gbs_dat@ped$id %in% yp_dat$fullsamplename,]
  gbs_amat<-GRM(gbs_dat, autosome.only = F)
  gbs_amat<-gbs_amat[rownames(gbs_amat) %in% unique(yp_dat$fullsamplename),
                     colnames(gbs_amat) %in% unique(yp_dat$fullsamplename)]
  
  #code data and remove pedigrees
  #convert matrix
  gbs_dat_code<-as.matrix(gbs_dat)
  
  #pull names and assign dummy
  names<-data.frame(fullsamplename=rownames(gbs_dat_code))
  names$Dummy=paste("Geno_",seq(1:nrow(names)),sep = "")
  
  #recode pheno data
  pheno_code<-yp_dat %>%
    left_join(names, by="fullsamplename") %>%
    select(-fullsamplename) %>%
    rename(fullsamplename=Dummy) %>%
    select(1:5, fullsamplename, everything()) %>%
    mutate(id="Not_Available",
           pedigree="Not_Available") %>%
    drop_na(fullsamplename)
  
  #recode gbs data
  gbs_dat_code<-as.data.frame(gbs_dat_code) %>%
    rownames_to_column(var = "fullsamplename") %>%
    left_join(names, by="fullsamplename") %>%
    select(-fullsamplename) %>%
    rename(fullsamplename=Dummy) %>%
    select(fullsamplename, everything())
  rownames(gbs_dat_code)=gbs_dat_code$fullsamplename
  gbs_dat_code<-gbs_dat_code %>%
    select(-fullsamplename)
  colnames(gbs_dat_code)=paste("Marker_", seq(1:ncol(gbs_dat_code)), sep = "")
  
  #recode grm
  gbs_amat_coded<-as.data.frame(gbs_amat) %>%
    rownames_to_column(var="fullsamplename") %>%
    left_join(names, by="fullsamplename") %>%
    select(-fullsamplename) %>%
    rename(fullsamplename=Dummy) %>%
    select(fullsamplename, everything())
  rownames(gbs_amat_coded)=gbs_amat_coded$fullsamplename
  gbs_amat_coded<-gbs_amat_coded %>%
    select(-fullsamplename)
  colnames(gbs_amat_coded)=rownames(gbs_amat_coded)

  
  #write out results
  write.table(pheno_code,
              gzfile("pheno_coded.txt.gz"),
              sep = "\t",
              row.names = FALSE,
              quote = FALSE)
  write.table(gbs_dat_code,
              gzfile("geno_coded.txt.gz"),
              sep = "\t",
              row.names = TRUE,
              quote = FALSE)
  write.table(gbs_amat_coded,
              gzfile("geno_amat_coded.txt.gz"),
              sep = "\t",
              row.names = TRUE,
              quote = FALSE)
  
  #read in files
  yp_dat<-read.table(gzfile("pheno_coded.txt.gz"),
                     header = TRUE)
  gbs_dat<-read.table(gzfile("geno_coded.txt.gz"),
                      header = TRUE)
  gbs_amat<-read.table(gzfile("geno_amat_coded.txt.gz"),
                       header = TRUE)
}
```

# Calculate per plot, narrow sense heritability of all traits in each field  

```{r}
#check files
file<-list.files()[grep("h2",list.files(), ignore.case = T)]

if(length(file)>=1){
  
  print("Heritability within environment per trait has already been calculated; moving on!")
  h2<-read.csv(file)
  remove(file)
  
}else{

  #remove
  remove(file)
  
  #empty vector
  h2<-c()
  
  
  for (i in unique(yp_dat$environment)){
    
    #print location
    print(i)
    
    #filter data
    a<-yp_dat %>%
      filter(environment==i)
    
    #filter GRM
    grm<-gbs_amat[rownames(gbs_amat) %in% unique(a$fullsamplename),
                  colnames(gbs_amat) %in% unique(a$fullsamplename)]
    
    #make sure the GRM and pheno data have same fullsamplenames
    a<-a %>%
      filter(fullsamplename %in% rownames(grm))
    
    for(j in colnames(a)[14:ncol(a)]){
      
      #print trait
      print(j)
      
      #pull data
      b<-a %>%
        select(fullsamplename, j) %>%
        drop_na()
      
      if(nrow(b)==0){
        
        #say something
        print(paste(j, "was not taken in", i))
        remove(b)
        
      }else{
        
        #format
        b[,1]=as.factor(b[,1])
        b[,2]=as.numeric(b[,2])
        
        #rename
        colnames(b)[2]="y"
        
        #run model
        fit<-try(asreml(fixed = y~1,
                        random = ~vm(fullsamplename, grm),
                        residual = ~idv(units),
                        data = b,
                        trace = FALSE,
                        maxit = 75))
        
        if(class(fit)=="try-error"){
          
          #say something
          print(paste(j, "caused a singlualtiry in the model"))
          remove(b, fit)
          
        }else{
          
          #calculate heritability
          h<-vpredict(fit, ~(V1/(V1+V2))) %>%
            mutate(Environment=i,
                   Trait=j) %>%
            select(Environment,
                   Trait,
                   everything())
          
          #bind in
          h2<-rbind(h2, h)
          
          #remove
          remove(b,fit,h)
          
        }
        
      }
      
    }
    
    #remove
    remove(a, grm)
    
  }
  
  #write out results
  write.csv(h2,
            "h2_within_environment.csv",
            row.names = F)
    
}
```

# Correlation with environment

```{r}
corrs<-c()

for(i in unique(yp_dat$environment)){
  
  a<-yp_dat %>%
    filter(environment==i)
  a<-a[,!colSums(is.na(a))>=nrow(a)-10]
  a<-a %>%
    drop_na()
  a<-rcorr(as.matrix(a[,14:ncol(a)]))
  b<-melt(a$r)
  colnames(b)[3]="r"
  c<-melt(a$P)
  colnames(c)[3]="P-Value"
  d<-b %>%
    left_join(c, by=c("Var1", "Var2"))
  d$`P-Value`<-round(d$`P-Value`, 4)
  Year=ifelse(length(grep("2019", i))>0, "2019",
              ifelse(length(grep("2020", i))>0, "2020", 
                     ifelse(length(grep("2021", i))>0, "2021", "error")))
  Location="Fort Collins"
  Trial=ifelse(length(grep("aynd1", i))>0, "AYN-D1",
               ifelse(length(grep("aynd2", i))>0, "AYN-D2", 
                      ifelse(length(grep("csu_elite", i))>0, "ELITE", 
                             ifelse(length(grep("ayn", i))>0, "AYN", 
                                    ifelse(length(grep("cl_pyn", i))>0, "CL-PYN", 
                                           ifelse(length(grep("_pyn", i))>0, "PYN", "error"))))))
  d<-data.frame(Environment=i,
                Year=Year,
                Location=Location,
                Trial=Trial,
                `Variable 1`=d$Var1,
                `Variable 2`=d$Var2,
                r=d$r,
                `P-Value`=d$`P-Value`,
                check.names = F)
  d<-d %>%
    drop_na(`P-Value`)
  
  corrs<-rbind(corrs,d)
  remove(a,b,c,d, Year, Location, Trial)

}


#look at yield within environment

for(i in unique(yp_dat$environment)){
  
  a<-corrs %>%
    filter(Environment==i) %>%
    filter(`Variable 2`=="grain_yield") %>%
    filter(`Variable 1` %in% unique(corrs$`Variable 1`)[grep("mean_", unique(corrs$`Variable 1`))]) %>%
    mutate(Significance=ifelse(`P-Value`<0.001, "P(t)<0.001",
                               ifelse(`P-Value`<0.01, "P(t)<0.01",
                                      ifelse(`P-Value`<0.05, "P(t)<0.05", "NS"))))
  a$Significance<-factor(a$Significance, levels = c("P(t)<0.001", "P(t)<0.01", "P(t)<0.05", "NS"))
  a$`Variable 1`=gsub("_mean_", "_", a$`Variable 1`)
  a$`Variable 1`=gsub("_flight", " (flight", a$`Variable 1`)
  a$`Variable 1`=paste(a$`Variable 1`, ")", sep = "")
  a$`Variable 1`=gsub("_", " ", a$`Variable 1`)
  a$`Variable 1`=str_to_title(a$`Variable 1`)
  a$`Variable 1`=gsub("Rededge", "RE", a$`Variable 1`)
  a$`Variable 1`=gsub("Nir", "NIR", a$`Variable 1`)
  a$`Variable 1`=gsub("Ndvi", "NDVI", a$`Variable 1`)
  a<-a[order(a$`Variable 1`),]
  b<-ggplot(data = a, aes(x=`Variable 1`, y=r, fill=Significance))+
    geom_col()+
    scale_fill_manual(values = c("NS" = "gray",
                                 "P(t)<0.05" = "pink",
                                 "P(t)<0.01" = "lightyellow",
                                 "P(t)<0.001" = "lightgreen"))+
    theme_bw()+
    theme(axis.text.x = element_text(angle = 90))+
    labs(title = paste("Correlation of Visual Indicies to Yield in the",
                       unique(a$Trial),
                       "Nursery in",
                       unique(a$Year)),
         x="Index",
         y="Correlation (r)")+
    geom_label(aes(label=round(r, 2)), nudge_y = ifelse(a$r>0, -.05, .05))
  print(b)
  ggsave(filename = paste("yield_corr_plot_",i, ".png", sep = ""),
         plot = b,
         width = 14,
         height = 4,
         units = "in",
         dpi = 320)
  remove(a,b)
    
}
```

# Calculate BLUEs within and among environments

```{r}
#format dataset
dat_ml<-yp_dat %>%
  select(environment, fullsamplename, row, col, grain_yield, 21:44) %>%
  drop_na() %>%
  filter(fullsamplename %in% rownames(gbs_amat))

#check files
file<-list.files()[grep("BLUEs", list.files(), ignore.case = T)]

if(length(file)==3){
  
  print("BLUEs have already been calculated across environments; moving on!")
  dat_ml_blues_envs<-read.csv(file[1])
  dat_ml_blues_year<-read.csv(file[2])
  dat_ml_blues<-read.csv(file[3])
  remove(file)
  
}else{
  
  dat_ml_blues_loc<-c()
  
  for(i in unique(dat_ml$environment)){
    
    #announce
    print(paste("Analyzing", i))
    
    #pull data
    a<-dat_ml %>%
      filter(environment==i)
    
    #place to put the means
    blues<-a %>%
      distinct(fullsamplename)
    
    for(j in colnames(a)[5:ncol(a)]){
      
      #announce
      print(paste("Analyzing", j))
      
      #pull data
      b<-a %>%
        select(fullsamplename, row, col, all_of(j))
      
      #format
      b[,1:3]<-lapply(b[,1:3], as.factor)
      b[,4]<-as.numeric(b[,4])
      
      #rename
      colnames(b)[4]<-"y"
      
      #fit 4 models
      fit1<-try(asreml(fixed = y ~ fullsamplename,
                       random = ~idv(row)+idv(col),
                       residual = ~idv(units),
                       data = b,
                       workspace = "8gb",
                       trace = F))
      
      pred<-as.data.frame(fit1$coefficients$fixed) %>%
        rownames_to_column(var = "fullsamplename") %>%
        mutate(fullsamplename=gsub("fullsamplename_", "", fullsamplename),
               fullsamplename=gsub("[(]", "", fullsamplename),
               fullsamplename=gsub("[)]", "", fullsamplename))
      bu<-pred %>%
        filter(fullsamplename=="Intercept")
      pred<-pred %>%
        filter(fullsamplename!="Intercept") %>%
        mutate(effect=effect+bu$effect)
      colnames(pred)[2]<-j
      blues<-blues %>%
        left_join(pred, by = "fullsamplename")
      
      remove(fit1,b,pred,bu)
    
    }
    
    blues<-blues %>%
      mutate(environment=i) %>%
      select(environment, everything())
    
    #put them on the big object
    dat_ml_blues_loc<-rbind(dat_ml_blues_loc, blues)
    
    #remove
    remove(a, blues)
    
  }
  
  #write out results
  write.csv(dat_ml_blues_loc,
            "BLUEs_within_environment.csv",
            row.names = F)
  
  dat_ml_blues_envs<-read.csv("BLUEs_within_environment.csv")
  
  #make dataset
  temp<-dat_ml_blues_envs %>%
    mutate(year=ifelse(environment %in% dat_ml_blues_envs$environment[grep("2019", dat_ml_blues_envs$environment)], 2019,
                       ifelse(environment %in% dat_ml_blues_envs$environment[grep("2020", dat_ml_blues_envs$environment)], 2020,
                              ifelse(environment %in% dat_ml_blues_envs$environment[grep("2021", dat_ml_blues_envs$environment)], 2021, "ERROR")))) %>%
    select(environment, year, fullsamplename, grain_yield, 4:ncol(dat_ml_blues_envs)) %>%
    drop_na()
  
  dat_ml_blues_year<-c()
  
  for(i in unique(temp$year)){
    
    print(paste("Analying all data in year =", i))
    
    a<-temp %>%
      filter(year==i)
    
    results<-a %>%
      distinct(fullsamplename, year)
    
    for(j in colnames(a)[4:ncol(a)]){
    
      print(paste("Analying", j, "in year =", i))
        
      b<-a %>%
        select(environment, fullsamplename, all_of(j)) 
      
      b[,1:2]<-lapply(b[,1:2], as.factor)
      b[,3]<-as.numeric(b[,3])
      colnames(b)[3]="y"
      
      fit<-asreml(fixed = y ~ fullsamplename,
                  random = ~idv(environment),
                  residual = ~idv(units),
                  data = b,
                  trace = T,
                  workspace = "8gb")
      
      pred<-data.frame(fit$coefficients$fixed) %>%
        rownames_to_column("id") %>%
        separate(id, into = c("factor", "fullsamplename1", "fullsamplename2"), sep = "_") %>%
        mutate(fullsamplename=paste(fullsamplename1, fullsamplename2, sep = "_")) %>%
        select(factor, fullsamplename, effect)
      u<-pred %>%
        filter(factor=="(Intercept)")
      pred<-pred %>%
        filter(!factor=="(Intercept)") %>%
        select(fullsamplename, effect) %>%
        rename(bu=effect) %>%
        mutate(bu=bu+u$effect)
      colnames(pred)[2]<-j
      
      results<-results %>%
        left_join(pred, by = "fullsamplename")
      
      remove(b,fit,pred)
      
    }
    
    dat_ml_blues_year<-rbind(dat_ml_blues_year, results)
    
    remove(a, results)
    
  }

  write.csv(dat_ml_blues_year,
            "BLUEs_within_year.csv",
            row.names = F)

  dat_ml_blues<-dat_ml_blues_envs %>%
    distinct(fullsamplename)

  for(i in colnames(dat_ml_blues_envs)[3:ncol(dat_ml_blues_envs)]){

    #announce
    print(paste("Calculating adjusted means across environments for", i))
    
    #pull data
    a<-dat_ml_blues_envs %>%
       select(environment, fullsamplename, all_of(i))
    colnames(a)[3]="y"
    
    #format columns
    a[,1:2]<-lapply(a[,1:2], as.factor)
    a[,3]<-as.numeric(a[,3])  
    
    #announce
    print("Fitting model...")
    
    #calculate blues across environments
    fit<-asreml(fixed = y ~ fullsamplename,
                random = ~idv(environment),
                residual = ~idv(units),
                data = a,
                workspace = "12gb",
                trace = T)
      
    #announce
    print("Estimating BLUEs...")
    
    pred<-predict(fit, 
                  classify = "fullsamplename", 
                  pworkspace="12gb",
                  trace = F)
    
    #format BLUEs
    pred<-pred$pvals[,1:2]
    colnames(pred)[2]=i
    
    #bind in the data
    dat_ml_blues<-dat_ml_blues %>%
      left_join(pred, by="fullsamplename")
    
    #remove
    remove(a,fit,pred)
    
    #announce
    print(paste("Finshed at:", Sys.time()))
      
  }
  
  #write results
  write.csv(dat_ml_blues,
            "ME_BLUEs.csv",
            row.names = F)
  
}
```

# PPE Methods

```{r}
##################################
###80:20 split cross validation###
##################################

#check files
file<-list.files()[grep("MLM_FFCV", list.files(), ignore.case = T)]

if("MLM_FFCV_raw_output.csv" %in% file &
   "MLM_FFCV_summary.csv" %in% file){
  
  print("Machine learning, 80:20 split cross validation is complete; Moving on!")
  
}else{
    
  #set number permutation
  perms=30
  cv_results=c()
  dat_ml<-dat_ml_blues %>%
    select(grain_yield, 3:ncol(dat_ml_blues)) %>%
    drop_na()
  
  #80:20 split
  cv_results=c()
  
  for(i in 1:perms){
    
    #announce
    print(paste("Commencing Permutation =", i))
    
    #select train and test
    train=sample(rownames(dat_ml), size = round(nrow(dat_ml)*0.8, 0))
    train=dat_ml[rownames(dat_ml) %in% train,]  
    test=dat_ml[!rownames(dat_ml) %in% rownames(train),]
    
    #make folds of the data
    myFolds <- createFolds(train[,1], k = 5)
    
    #tune parameters
    myControl <- trainControl(
      classProb = FALSE,
      verboseIter = FALSE,
      savePredictions = TRUE,
      index = myFolds,
      number = 5,
      repeats = 1000,
      method = "repeatedcv")
    
    #create model statement
    q<-paste(colnames(train)[1], " ~ .", sep = "")
    
    #create k
    k_knn=seq(0, 50, by=5)
    k_knn=c(1, k_knn[2:length(k_knn)])
    k_rf=c(1,5,10,15,20)
    
    #linear regression
    lm_model <- train(formula(q), 
                      data = train,
                      method = "lm",
                      trControl = myControl)
    
    #knn model
    knn_model <- train(formula(q), 
                       data = train,
                       method = "knn",
                       tuneGrid = expand.grid(k = k_knn),
                       trControl = myControl)
    
    
    #rf model
    rf_model <- train(formula(q),
                      train,
                      method = "rf",
                      tuneGrid = expand.grid(mtry = k_rf),
                      trControl = myControl)
    
    #print summary
    print(lm_model)
    print(knn_model)
    print(rf_model)
    
    #make prediction
    pred_lm<-predict(lm_model, test[,-1])
    pred_knn<-predict(knn_model, test[,-1])
    pred_rf<-predict(rf_model, test[,-1])
    
    for(j in c("pred_lm", "pred_knn", "pred_rf")){
    
      if(j=="pred_lm"){
        a<-pred_lm
        name<-"Linear Regression"
      }else if (j=="pred_knn"){
        a<-pred_knn
        name<-"K-Nearest Neighbors"
      }else if (j=="pred_rf"){
        a<-pred_rf
        name<-"Random Forest"
      }else{
        print("Something went wrong...")
      }
      
      #get regression
      fit<-lm(a~test[,1])
      r=summary(fit)$r.squared
      Fstat=summary(fit)$fstatistic
      rse=summary(fit)$sigma
      
      #correlate
      plot(x=a, 
           y=test[,1], 
           xlab = "Predicted Grain Yield (Bu/A)", 
           ylab = "Observed Grain Yield (Bu/A)",
           main = paste("(Perm = ",
                        i,
                        ") Scatterplot of Predicted vs Observed Grain Yield (Bu/A) (", 
                        name,
                        ")",
                        sep = ""))
      
      legend("bottomright" ,legend=c(paste("r =", round(sqrt(r), 2)),
                                     paste("r2 =", round(r, 2)),
                                     paste("F-Statistic =", round(Fstat[1], digits = 1)),
                                     paste("DFn =", round(Fstat[2], digits = 1)),
                                     paste("DFd =", round(Fstat[3], digits = 1))))  
      
      #make dataframe
      a<-data.frame(Permutation=i,
                    Model=name,
                    r = sqrt(r),
                    r2 = r,
                    RSE = rse,
                    `F-Statistic`=Fstat[1],
                    DFn=Fstat[2],
                    DFd=Fstat[3],
                    check.names = F)
      
      #bind
      cv_results<-rbind(cv_results, a)
      
      #remove
      remove(a, 
             name,
             fit,
             r,
             Fstat,
             rse)
      
    }
    
    #remove
    remove(train, 
           test, 
           myFolds,
           myControl,
           q,
           k_knn,
           k_rf,
           lm_model,
           knn_model,
           rf_model,
           pred_lm,
           pred_knn,
           pred_rf)
    
  }
  
  #summarize
  ff_cv_sum<-cv_results %>%
    group_by(Model)%>%
    summarise(`Min r`=min(r),
              `Max r`=max(r),
              `Mean r`=mean(r),
              `SD r`=sd(r),
              `Min r2`=min(r2),
              `Max r2`=max(r2),
              `Mean r2`=mean(r2),
              `SD r2`=sd(r2),
              `Min RSE`=min(RSE),
              `Max RSE`=max(RSE),
              `Mean RSE`=mean(RSE),
              `SD RSE`=sd(RSE))
  
  #write out results
  write.csv(cv_results,
            "MLM_FFCV_raw_output.csv",
            row.names = F)
  write.csv(ff_cv_sum,
            "MLM_FFCV_summary.csv",
            row.names = F)
  
}

##########################################
###leave one trial out validation###
##########################################

#check files
file<-list.files()[grep("MLM_LOEO", list.files(), ignore.case = T)]

if(length(file)>0){
  
  #announce
  print("leave one trial out validation has been completed; moving on!")
  remove(file)
  
}else{

  #make empty vector
  loeo_results=c()
  loeo_seeds=c()
  
  #format data
  dat_ml<-dat_ml_blues_envs %>%
    select(environment, fullsamplename, grain_yield, 4:ncol(dat_ml_blues_envs)) %>%
    drop_na()
  
  for(i in unique(dat_ml$environment)){
    
    #announce
    print(paste("Predicting", i))
    
    #select train and test
    test=dat_ml %>% 
      filter(environment==i) %>% 
      select(-environment) 
    train=dat_ml %>% 
      filter(!environment==i) %>% 
      select(-environment) %>% 
      filter(!fullsamplename %in% test$fullsamplename)
    train=train %>% select(-fullsamplename)
    test=test %>% select(-fullsamplename)
    
    #set seed
    seed<-data.frame(seed=seq(1:10000000), rand=rnorm(10000000, mean = 0, sd = 1))
    seed<-seed[order(seed$rand),]
    seed<-seed[1,1]
    set.seed(seed)
    seed<-data.frame(Environment=i,
                     Seed=seed)
    loeo_seeds<-rbind(loeo_seeds, seed)
    remove(seed)
  
    #make folds of the data
    myFolds <- createFolds(train[,1], k = 5)
    
    #tune parameters
    myControl <- trainControl(
      classProb = FALSE,
      verboseIter = FALSE,
      savePredictions = TRUE,
      index = myFolds,
      number = 5,
      repeats = 1000,
      method = "repeatedcv")
    
    #create model statement
    q<-paste(colnames(train)[1], " ~ .", sep = "")
    
    #create k
    k_knn=seq(0, 50, by=5)
    k_knn=c(1, k_knn[2:length(k_knn)])
    k_rf=c(1,5,10,15,20)
    
    #linear regression
    lm_model <- train(formula(q), 
                      data = train,
                      method = "lm",
                      trControl = myControl)
    
    #knn model
    knn_model <- train(formula(q), 
                       data = train,
                       method = "knn",
                       tuneGrid = expand.grid(k = k_knn),
                       trControl = myControl)
    
    
    #rf model
    rf_model <- train(formula(q),
                      train,
                      method = "rf",
                      tuneGrid = expand.grid(mtry = k_rf),
                      trControl = myControl)
    
    #nullify seed
    set.seed(NULL)
    
    #print summary
    print(lm_model)
    print(knn_model)
    print(rf_model)
    
    #make prediction
    pred_lm<-predict(lm_model, test[,-1])
    pred_knn<-predict(knn_model, test[,-1])
    pred_rf<-predict(rf_model, test[,-1])
    
    for(j in c("pred_lm", "pred_knn", "pred_rf")){
    
      if(j=="pred_lm"){
        a<-pred_lm
        name<-"Linear Regression"
      }else if (j=="pred_knn"){
        a<-pred_knn
        name<-"K-Nearest Neighbors"
      }else if (j=="pred_rf"){
        a<-pred_rf
        name<-"Random Forest"
      }else{
        print("Something went wrong...")
      }
      
      #get regression
      fit<-lm(a~test[,1])
      r=summary(fit)$r.squared
      Fstat=summary(fit)$fstatistic
      rse=summary(fit)$sigma
      
      #correlate
      plot(x=a, 
           y=test[,1], 
           xlab = "Predicted Grain Yield (Kg/Ha)", 
           ylab = "Observed Grain Yield (Kg/Ha)",
           main = paste("(Environment = ",
                        i,
                        ") Scatterplot of Predicted vs Observed Grain Yield (Kg/Ha) (", 
                        name,
                        ")",
                        sep = ""))
      
      legend("bottomright" ,legend=c(paste("r =", round(sqrt(r), 2)),
                                     paste("r2 =", round(r, 2)),
                                     paste("F-Statistic =", round(Fstat[1], digits = 1)),
                                     paste("DFn =", round(Fstat[2], digits = 1)),
                                     paste("DFd =", round(Fstat[3], digits = 1))))  
      
      #make dataframe
      a<-data.frame(Environment=i,
                    Model=name,
                    r = sqrt(r),
                    r2 = r,
                    RSE = rse,
                    `F-Statistic`=Fstat[1],
                    DFn=Fstat[2],
                    DFd=Fstat[3],
                    check.names = F)
      
      #bind
      loeo_results<-rbind(loeo_results, a)
      
      #remove
      remove(a, 
             name,
             fit,
             r,
             Fstat,
             rse)
      
    }
    
    #remove
    remove(train, 
           test, 
           myFolds,
           myControl,
           q,
           k_knn,
           k_rf,
           lm_model,
           knn_model,
           rf_model,
           pred_lm,
           pred_knn,
           pred_rf)
    
  }
  
  #format results
  loeo_results<-loeo_results %>%
    mutate(Year=ifelse(Environment %in% unique(loeo_results$Environment[grep("2019", loeo_results$Environment)]), "2019", 
                       ifelse(Environment %in% unique(loeo_results$Environment[grep("2020", loeo_results$Environment)]), "2020",
                              ifelse(Environment %in% unique(loeo_results$Environment[grep("2021", loeo_results$Environment)]), "2021", "ERROR"))),
           Trial=ifelse(Environment %in% unique(loeo_results$Environment[grep("aynd1", loeo_results$Environment)]), "AYN-D1",
                        ifelse(Environment %in% unique(loeo_results$Environment[grep("aynd2", loeo_results$Environment)]), "AYN-D2",
                               ifelse(Environment %in% unique(loeo_results$Environment[grep("csu_elite", loeo_results$Environment)]), "ELITE",
                                      ifelse(Environment %in% unique(loeo_results$Environment[grep("_pyn", loeo_results$Environment)]), "PYN",
                                             ifelse(Environment %in% unique(loeo_results$Environment[grep("ayn", loeo_results$Environment)]), "AYN",
                                                    ifelse(Environment %in% unique(loeo_results$Environment[grep("cl_pyn", loeo_results$Environment)]), "PYN", "ERROR"))))))) %>%
    select(-Environment) %>%
    select(Year, Trial, everything()) %>%
    arrange(Year, Trial)
  
  #write out results
  write.csv(loeo_results,
            "MLM_LOEO_results.csv",
            row.names = F)
  
  #write out results
  write.csv(loeo_seeds,
            "MLM_LOEO_seeds.csv",
            row.names = F)
}


###################################
###leave one year out validation###
###################################

#check files
file<-list.files()[grep("MLM_LOYO", list.files(), ignore.case = T)]

if(length(file)>0){
  
  #announce
  print("Leave one year out validation has been completed; moving on!")
  remove(file)
  
}else{

  #make empty vector
  loyo_results=c()
  loyo_seeds=c()

  #format data
  dat_ml=dat_ml_blues_year
  
  for(i in unique(dat_ml$year)){
    
    #announce
    print(paste("Predicting", i))
    
    #select train and test
    test=dat_ml %>% 
      filter(year==i) %>% 
      select(-year) 
    train=dat_ml %>% 
      filter(!year==i) %>% 
      select(-year) %>% 
      filter(!fullsamplename %in% test$fullsamplename)
    train=train %>% select(-fullsamplename)
    test=test %>% select(-fullsamplename)
    
    #set seed
    seed<-data.frame(seed=seq(1:10000000), rand=rnorm(10000000, mean = 0, sd = 1))
    seed<-seed[order(seed$rand),]
    seed<-seed[1,1]
    set.seed(seed)
    seed<-data.frame(Year=i,
                     Seed=seed)
    loyo_seeds<-rbind(loyo_seeds, seed)
    remove(seed)

    #make folds of the data
    myFolds <- createFolds(train[,1], k = 5)
    
    #tune parameters
    myControl <- trainControl(
      classProb = FALSE,
      verboseIter = FALSE,
      savePredictions = TRUE,
      index = myFolds,
      number = 5,
      repeats = 1000,
      method = "repeatedcv")
    
    #create model statement
    q<-paste(colnames(train)[1], " ~ .", sep = "")
    
    #create k
    k_knn=seq(0, 50, by=5)
    k_knn=c(1, k_knn[2:length(k_knn)])
    k_rf=c(1,5,10,15,20)
    
    #linear regression
    lm_model <- train(formula(q), 
                      data = train,
                      method = "lm",
                      trControl = myControl)
    
    #knn model
    knn_model <- train(formula(q), 
                       data = train,
                       method = "knn",
                       tuneGrid = expand.grid(k = k_knn),
                       trControl = myControl)
    
    
    #rf model
    rf_model <- train(formula(q),
                      train,
                      method = "rf",
                      tuneGrid = expand.grid(mtry = k_rf),
                      trControl = myControl)
    
    #print summary
    print(lm_model)
    print(knn_model)
    print(rf_model)
    
    #unset seed
    set.seed(NULL)
    
    #make prediction
    pred_lm<-predict(lm_model, test[,-1])
    pred_knn<-predict(knn_model, test[,-1])
    pred_rf<-predict(rf_model, test[,-1])
    
    for(j in c("pred_lm", "pred_knn", "pred_rf")){
    
      if(j=="pred_lm"){
        a<-pred_lm
        name<-"Linear Regression"
      }else if (j=="pred_knn"){
        a<-pred_knn
        name<-"K-Nearest Neighbors"
      }else if (j=="pred_rf"){
        a<-pred_rf
        name<-"Random Forest"
      }else{
        print("Something went wrong...")
      }
      
      #get regression
      fit<-lm(a~test[,1])
      r=summary(fit)$r.squared
      Fstat=summary(fit)$fstatistic
      rse=summary(fit)$sigma
      
      #correlate
      plot(x=a, 
           y=test[,1], 
           xlab = "Predicted Grain Yield (Kg/Ha)", 
           ylab = "Observed Grain Yield (Kg/Ha)",
           main = paste("(Year = ",
                        i,
                        ") Scatterplot of Predicted vs Observed Grain Yield (Kg/Ha) (", 
                        name,
                        ")",
                        sep = ""))
      
      legend("bottomright" ,legend=c(paste("r =", round(sqrt(r), 2)),
                                     paste("r2 =", round(r, 2)),
                                     paste("F-Statistic =", round(Fstat[1], digits = 1)),
                                     paste("DFn =", round(Fstat[2], digits = 1)),
                                     paste("DFd =", round(Fstat[3], digits = 1))))  
      
      #make dataframe
      a<-data.frame(Year=i,
                    Model=name,
                    r = sqrt(r),
                    r2 = r,
                    RSE = rse,
                    `F-Statistic`=Fstat[1],
                    DFn=Fstat[2],
                    DFd=Fstat[3],
                    check.names = F)
      
      #bind
      loyo_results<-rbind(loyo_results, a)
      
      #remove
      remove(a, 
             name,
             fit,
             r,
             Fstat,
             rse)
      
    }
    
    #remove
    remove(train, 
           test, 
           myFolds,
           myControl,
           q,
           k_knn,
           k_rf,
           lm_model,
           knn_model,
           rf_model,
           pred_lm,
           pred_knn,
           pred_rf)
    
  }
  
  #format results
  loyo_results<-loyo_results %>%
    select(Year, everything()) %>%
    arrange(Year)
  
  #write out results
  write.csv(loyo_results,
            "MLM_LOYO_results.csv",
            row.names = F)
  
  write.csv(loyo_seeds,
            "MLM_LOYO_seeds.csv",
            row.names = F)
}
```

# GEBV methods

```{r}
#check files
file<-list.files()[grep("GP_FFCV", list.files(), ignore.case = T)]

if(length(file)>0){
  
  print("Genomic prediction, 80:20 split cross validation is complete; Moving on!")
  
}else{
    
  #set number permutation
  perms=30
  
  cv_results<-c()

  for(i in 1:perms){
    
    #print
    print(paste("Permutation", i))
    
    #select train and test
    train=sample(unique(dat_ml_blues$fullsamplename), size = round(length(unique(dat_ml_blues$fullsamplename))*0.8, 0))
    test=dat_ml_blues[!dat_ml_blues$fullsamplename %in% train,"fullsamplename"]
    
    #pull data
    a<-dat_ml_blues[,c("fullsamplename", "grain_yield")]
    colnames(a)[2]="y"
    
    #format columns
    a[,1]<-as.factor(a[,1])
    a[,2]<-as.numeric(a[,2])  
    
    #drop test observations
    a$y<-ifelse(a$fullsamplename %in% test, NA, a$y)
    
    #pull gbs matrix
    b<-gbs_amat[rownames(gbs_amat) %in% a$fullsamplename,
                colnames(gbs_amat) %in% a$fullsamplename]
    b<-as.matrix(b)
  
    #assign to global environment
    assign("a", a, envir = globalenv())
    assign("b", b, envir = globalenv())
    
    #run model
    fit<-asreml::asreml(fixed = y~1,
                        random = ~vm(fullsamplename,b),
                        residual = ~idv(units),
                        data = a,
                        workspace = "4gb")
            
    #predict
    pred<-asreml::predict.asreml(fit, 
                                 classify = "fullsamplename", 
                                 pworkspace="4gb")$pvals[,1:2] 
    pred<-pred[pred$fullsamplename %in% test,]
    tbv<-dat_ml_blues[dat_ml_blues$fullsamplename %in% test,c("fullsamplename", "grain_yield")]
    pred<-merge(pred, tbv, by="fullsamplename")
    
    a<-summary(lm(pred[,2]~pred[,3]))
    a<-data.frame(Permutation=i,
                  Model="gBLUP",
                  r=sqrt(a$r.squared),
                  r2=a$r.squared,
                  RSE=a$sigma,
                  `F-Statistic`=a$fstatistic[1],
                  DFn=a$fstatistic[2],
                  DFd=a$fstatistic[3])
    
    cv_results<-rbind(cv_results, a)

    remove(a,b,fit,pred,tbv)
  }
  
  #summarize
  ff_cv_sum<-cv_results %>%
    group_by(Model)%>%
    summarise(`Min r`=min(r),
              `Max r`=max(r),
              `Mean r`=mean(r),
              `SD r`=sd(r),
              `Min r2`=min(r2),
              `Max r2`=max(r2),
              `Mean r2`=mean(r2),
              `SD r2`=sd(r2),
              `Min RSE`=min(RSE),
              `Max RSE`=max(RSE),
              `Mean RSE`=mean(RSE),
              `SD RSE`=sd(RSE))
  
  #write out results
  write.csv(cv_results,
            "GP_FFCV_raw_output.csv",
            row.names = F)
  write.csv(ff_cv_sum,
            "GP_FFCV_summary.csv",
            row.names = F)
  
}

#leave one trial out
file<-list.files()[grep("GP_LOEO", list.files(), ignore.case = T)]

if(length(file)>0){
  
  print("Genomic prediction, leave one trial out validation is complete; Moving on!")
  
}else{
  
  n_cores = detectCores()-1
  
  #Leave one core to avoid overload your computer
  cluster<-makeCluster(n_cores)
  registerDoParallel(cluster)

  loeo_results<-c()
  for(i in unique(dat_ml_blues_envs$environment)){
    
    #select train and test
    test=dat_ml_blues_envs[dat_ml_blues_envs$environment==i, c("environment", "fullsamplename", "grain_yield")]
    train=dat_ml_blues_envs[!dat_ml_blues_envs$environment==i, c("environment", "fullsamplename", "grain_yield")]
    train=train[!train$fullsamplename %in% test$fullsamplename,]
    train=rbind(train, test)
    train$grain_yield=ifelse(train$environment==i, NA, train$grain_yield)
  
    #format columns
    a<-train
    a[,1:2]<-lapply(a[,1:2], as.factor)
    a[,3]<-as.numeric(a[,3])
    
    #pull gbs matrix
    b<-gbs_amat[rownames(gbs_amat) %in% a$fullsamplename,
                colnames(gbs_amat) %in% a$fullsamplename]
    b<-as.matrix(b)
  
    #assign to global environment
    assign("a", a, envir = globalenv())
    assign("b", b, envir = globalenv())
    
    #run model
    fit<-asreml::asreml(fixed = grain_yield~1,
                        random = ~vm(fullsamplename,b)+idv(environment),
                        residual = ~idv(units),
                        data = a,
                        workspace = "4gb")
                    
    #predict
    pred<-asreml::predict.asreml(fit, 
                                 classify = "fullsamplename", 
                                 pworkspace="4gb")$pvals[,1:2] 
    pred<-pred[pred$fullsamplename %in% test$fullsamplename,]
    tbv<-dat_ml_blues[dat_ml_blues$fullsamplename %in% test$fullsamplename,
                      c("fullsamplename", "grain_yield")]
    pred<-merge(pred, tbv, by="fullsamplename")
    
    a<-summary(lm(pred[,2]~pred[,3]))
    a<-data.frame(Environment=i,
                  Model="gBLUP",
                  r=sqrt(a$r.squared),
                  r2=a$r.squared,
                  RSE=a$sigma,
                  `F-Statistic`=a$fstatistic[1],
                  DFn=a$fstatistic[2],
                  DFd=a$fstatistic[3])
    loeo_results<-rbind(a,loeo_results)

    remove(a,b,fit,pred,tbv)
  }
  
  #format
  loeo_results<-loeo_results %>%
    mutate(Year=ifelse(Environment %in% unique(loeo_results$Environment[grep("2019", loeo_results$Environment)]), "2019", 
                       ifelse(Environment %in% unique(loeo_results$Environment[grep("2020", loeo_results$Environment)]), "2020",
                              ifelse(Environment %in% unique(loeo_results$Environment[grep("2021", loeo_results$Environment)]), "2021", "ERROR"))),
           Trial=ifelse(Environment %in% unique(loeo_results$Environment[grep("aynd1", loeo_results$Environment)]), "AYN-D1",
                        ifelse(Environment %in% unique(loeo_results$Environment[grep("aynd2", loeo_results$Environment)]), "AYN-D2",
                               ifelse(Environment %in% unique(loeo_results$Environment[grep("csu_elite", loeo_results$Environment)]), "ELITE",
                                      ifelse(Environment %in% unique(loeo_results$Environment[grep("_pyn", loeo_results$Environment)]), "PYN",
                                             ifelse(Environment %in% unique(loeo_results$Environment[grep("ayn", loeo_results$Environment)]), "AYN",
                                                    ifelse(Environment %in% unique(loeo_results$Environment[grep("cl_pyn", loeo_results$Environment)]), "PYN", "ERROR"))))))) %>%
    select(-Environment) %>%
    select(Year, Trial, everything())
  
  #write out results
  write.csv(loeo_results,
            "GP_LOEO_results.csv",
            row.names = F)
}

###################################
###leave one year out validation###
###################################

#check files
file<-list.files()[grep("GP_LOYO", list.files(), ignore.case = T)]

if(length(file)>0){
  
  #announce
  print("Leave one year out validation has been completed; moving on!")
  remove(file)
  
}else{

  #make empty vector
  loyo_results=c()

  #format data
  dat_ml=dat_ml_blues_year
  
  for(i in unique(dat_ml$year)){
    
    #announce
    print(paste("Predicting", i))
    
    #select train and test
    test=dat_ml %>% filter(year==i) %>% select(year, fullsamplename, grain_yield)
    train=dat_ml %>% filter(!year==i) %>% select(year, fullsamplename, grain_yield)
    train=train %>% filter(!fullsamplename %in% test$fullsamplename)
    train=rbind(train, test)
    train=train %>% mutate(grain_yield=ifelse(year==i, NA, grain_yield))
    
    #format columns
    train[,1:2]<-lapply(train[,1:2], as.factor)
    train[,3]<-as.numeric(train[,3])  
    
    #pull gbs matrix
    b<-gbs_amat[rownames(gbs_amat) %in% train$fullsamplename,
                colnames(gbs_amat) %in% train$fullsamplename]
    b<-as.matrix(b)
  
    #run model
    fit<-asreml(fixed = grain_yield~1,
                random = ~vm(fullsamplename,b)+idv(year),
                residual = ~idv(units),
                data = train,
                workspace = "12gb")
    
    #predict
    pred<-data.frame(fit$coefficients$random) %>%
      rownames_to_column("id") %>%
      separate(id, into = c("trash", "fullsamplename_1", "fullsamplename_2"), sep = "_") %>%
      filter(trash!="year") %>%
      mutate(fullsamplename=paste(fullsamplename_1, fullsamplename_2, sep = "_")) %>%
      select(fullsamplename, effect)
    u<-as.numeric(fit$coefficients$fixed)
    pred<-pred %>%
      mutate(effect=effect+u)
    colnames(pred)[2]<-"predicted.value"
  
    
    pred<-pred %>%
      filter(fullsamplename %in% test$fullsamplename) %>%
      left_join(test, by="fullsamplename") %>%
      rename(`Grain Yield GEBV`=predicted.value,
             `Grain Yield TBV`=grain_yield) %>%
      select(fullsamplename, 
             `Grain Yield TBV`,
             `Grain Yield GEBV`)
    
    a<-ggplot(pred, aes(x=`Grain Yield TBV`, y=`Grain Yield GEBV`))+
      geom_point()+
      theme_bw()+
      labs(title = paste("(Year = ",
                         i,
                         ") True Breeding Value vs Genomic Estimated Breeding Value"),
           y="GEBV (Kg/Ha)",
           x="TBV (Kg/Ha)")
    print(a)
    
    a<-summary(lm(pred[,2]~pred[,3]))
    a<-data.frame(Year=i,
                  Model="gBLUP",
                  r=sqrt(a$r.squared),
                  r2=a$r.squared,
                  RSE=a$sigma,
                  `F-Statistic`=a$fstatistic[1],
                  DFn=a$fstatistic[2],
                  DFd=a$fstatistic[3])
    
    loyo_results<-rbind(loyo_results, a)
    
    remove(a, b, fit, pred, u, train, test)

  }
  
  #format results
  loyo_results<-loyo_results %>%
    select(Year, everything())
  
  #write out results
  write.csv(loyo_results,
            "GP_LOYO_results.csv",
            row.names = F)
}
```