## ----download-----------------------------------------------------------------
library(agro)
ff <- get_data_from_uri("https://data.cipotato.org/dataset.xhtml?persistentId=doi:10.21223/P3/M0HGJ4", ".")
ff


## -----------------------------------------------------------------------------
cluster <- read.delim(ff[1], stringsAsFactors=FALSE)
dict <- read.delim(ff[2], stringsAsFactors=FALSE)
library(readxl)
drought <- read_excel(ff[3])


## ----warn---------------------------------------------------------------------
drought$RYTHA[drought$RYTHA == "*"] <- NA
drought$FYTHA[drought$FYTHA == "*"] <- NA
drought$BIOM[drought$BIOM == "*"] <- NA


## ----convnumb-----------------------------------------------------------------
drought$RYTHA <- as.numeric(drought$RYTHA)
drought$FYTHA <- as.numeric(drought$FYTHA)
drought$BIOM  <- as.numeric(drought$BIOM)


## ----remna--------------------------------------------------------------------
drought <- drought[!is.na(drought$RYTHA),]


## ----cleaning-----------------------------------------------------------------
drought$env <- as.factor(paste(drought$YEAR, drought$TREATMENT))


## ----histogram----------------------------------------------------------------
library(ggplot2)
ggplot(drought, aes(x=RYTHA)) +
  geom_histogram(bins=10) +
  facet_wrap(~drought$env, nrow=1) 



## ----boxplottreatments--------------------------------------------------------
ggplot(drought, aes(y=RYTHA, x=TREATMENT, fill = REP)) + 
  geom_boxplot() + 
  theme_classic()


## ----boxplotcultivars---------------------------------------------------------
#New variable for type of cultivar
cluster$type <- as.numeric(row.names(cluster))
#The first 40 cultivars are landrace, the rest are modern (refer to table 5)
cluster$type <- ifelse(cluster$type <41, "landrace", "modern")

#Remove x-asis title (on first two) and legend (on all) to avoid repetition
p1 <- ggplot(cluster, aes(y=DM, x=type, color=type)) +
  geom_boxplot()+ 
  theme(legend.position="none", axis.title.x=element_blank()) 

p2 <- ggplot(cluster, aes(y=BIOM, x=type, color=type)) +
  geom_boxplot()+ 
  theme(legend.position="none", axis.title.x=element_blank())

p3 <- ggplot(cluster, aes(y=HI, x=type, color=type)) +
  geom_boxplot()+ 
  theme(legend.position="none")

p4 <- ggplot(cluster, aes(y=GM, x=type, color=type)) +
  geom_boxplot()+ 
  theme(legend.position="none")

g_legend <- function(a.gplot) {
  tpm <- ggplot_gtable(ggpl)
}

library(gridExtra)
grid.arrange(p1, p2, p3, p4, nrow=2)


## ----cluster------------------------------------------------------------------
row.names(cluster) <- cluster$Cultivar
j <- hclust(dist(cluster[,c(2,3,5,7,8)]), method="average")
plot(j, xlab="Cultivar", hang= -1, cex = 0.6)


## ----ANOVA--------------------------------------------------------------------
#Storage Root Yield
RYTHA.anova <- aov(RYTHA ~ env + CULTIVAR + CULTIVAR*env, data=drought)
summary(RYTHA.anova)

#Vine Yield
FYTHA.anova <- aov(FYTHA ~ env + CULTIVAR + CULTIVAR*env, data=drought)
summary(FYTHA.anova)

#Biomass
BIOM.anova <- aov(BIOM ~ env + CULTIVAR + CULTIVAR*env, data=drought)
summary(BIOM.anova)


## ----regression---------------------------------------------------------------
#create dummy variables for the different environments
drought$IR6 <- ifelse(drought$env=="2006 Irrigation", 1, 0)
drought$NI6 <- ifelse(drought$env=="2006 Non_Irrigation", 1, 0)
drought$IR8 <- ifelse(drought$env=="2008 Irrigation", 1, 0)
drought$NI8 <- ifelse(drought$env=="2008 Non_Irrigation", 1, 0)

#To see the effect of each environment on root yield, we can use the linear model with root yield as the y axis, and each environment as the x-axis. 
linear <- lm(RYTHA ~IR6 + NI6 + IR8, data=drought)
summary(linear)  


#Because we have so many cultivars, we can make them into dummy variables with the following code: 
cultivar.f = factor(drought$CULTIVAR)
dummies = model.matrix(~cultivar.f)

linear2 <- lm(drought$RYTHA ~ dummies)
#summary(linear2)



## ----ammi---------------------------------------------------------------------
library(agricolae)
#To make for easier viewing, we will use the codes rather than cultivar names, as charactes.
drought$CODE_CULTIVAR <- as.character(drought$CODE_CULTIVAR)
#the function AMMI specifies the environment variable, genotype, and repition, as well as the y-variable of interest.
ammi <- AMMI(ENV=drought$env, GEN=drought$CODE_CULTIVAR, REP=drought$REP, Y = drought$RYTHA, PC=T)
ammi$ANOVA
ammi$analysis


## ----biplot-------------------------------------------------------------------
#The regular plot(ammi) will plot PC1 against PC2
plot(ammi)
#To graph the dependent variable against pc1, the following code is used. 0 is the y-dependent variable (RYTHA), and 1 is PC1. 
plot(ammi, 0, 1, gcol="blue", ecol="green")