rm(list=ls())
options(max.print=2000000)
setwd("/Users/beatriceelohorifie/Documents/Module 10")

#################################################################
##                    Install R packages                       ##
##                      A smart function                       ##
#################################################################

check.packages <- function(pkg){
  new.pkg <- pkg[!(pkg%in%installed.packages()[,"Package"])]
  if (length(new.pkg))
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}
packages <- c('lme4' )
check.packages(packages)

ds1<-read.csv("QG_Mod10_ALA10.3_env_gxe_stdize_ds.csv",header=T)
#cluster environments based on pheno
distance<-dist(ds1,method="euclidean")
group<-hclust(distance,method="ward.D")
plot(group,main="Cluster of envs from ALA 10.1 
with similar Yield responses across lines",xlab="environments")

dslong<-read.csv("QG_Mod10_ALAs_ds.csv",header=T)
summary(dslong)
dslong$line<-factor(dslong$line)
dslong$env<-factor(dslong$env)
dslong$blk<-factor(dslong$blk)
summary(dslong)

dsenvAC<-subset(dslong,env=="A" | env=="C")
dsenvEF<-subset(dslong,env=="E" | env=="F")
dsenvIJ<-subset(dslong,env=='I' | env=="J")
dsenvB<-subset(dslong,env=="B")
dsenvD<-subset(dslong,env=="D")
dsenvG<-subset(dslong,env=="G")
dsenvH<-subset(dslong,env=="H")

gp1<-lmer(data=dsenvAC, Pheno ~ env + (1|line) )
gp1
gpest<-ranef(gp1)
gpest
write.csv(gpest,"line effect estimates for env A and C.csv", row.names = F)

gp2<-lmer(data=dsenvEF, Pheno ~ env + (1|line) )
gpest<-ranef(gp2)
write.csv(gpest,"line effect estimates for env E and F.csv", row.names = F)

gp3<-lmer(data=dsenvIJ, Pheno ~ env + (1|line) )
gpest<-ranef(gp3)
write.csv(gpest,"line effect estimates for env I and J.csv", row.names = F)

gpH<-lmer(data=dsenvH, Pheno ~ blk + (1|line) )
gpest<-ranef(gpH)
write.csv(gpest,"line effect estimates for env H.csv", row.names = F)

all<-lmer(data=dslong, Pheno ~ env + (1|line) )
est<-ranef(all)
write.csv(est,"line effect estimates across all envs.csv", row.names = F)