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

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','lsmeans')
check.packages(packages)

#Next load in the data set from a spreadsheet (excel.csv)
ds1 <-read.csv("QG_Mod13_ALA13.1_ds.csv",header=T)
summary(ds1)
ds1$Sampleid<-factor(ds1$Sampleid)
ds1$Env<-factor(ds1$Env)
ds1$Rep<-factor(ds1$Rep)
summary(ds1)


######################################################
### Conduct Exploratory Data Analyses by location ####
### Some useful tools are histograms and boxplots ####
######################################################

hist(ds1$yld,col="beige",main="Frequency of yld",
     xlab="Phenotypic value")

boxplot(ds1$yld~ds1$Env,main="Boxplot of yield in four environments",
        xlab="Environment",ylab="Yield")

mean(yld)
var(yld)
sd(yld)

#create new data sets for each environment
env1<-subset(ds1,ds1$Env=='E1')
env2<-subset(ds1,ds1$Env=='E2')
env3<-subset(ds1,ds1$Env=='E3')
env4<-subset(ds1,ds1$Env=='E4')
detach(ds1)

attach(env1)
hist(env1$yld,col="beige",
     main="Frequency of yield in Env E1",
     xlab="Yield Values")
detach(env1)

attach(env2)
hist(env2$yld,col="beige",
     main="Frequency of yield in Env E2",
     xlab="Yield Values")
detach(env2)

attach(env3)
hist(env3$yld,col="beige",
     main="Frequency of yield in Env E3",
     xlab="Yield Values")
detach(env3)

attach(env4)
hist(env4$yld,col="beige",
     main="Frequency of yield in Env E4",
     xlab="Yield Values")
detach(env4)

##########################################
#### Estimating the correlations      #### 
#### between all 6 pairs of locations ####
##########################################

plot(env1$yld,env2$yld,pch=20, 
     main="Correlation of yield in E1 and E2",
     xlab="Env 1", ylab="Env 2")
cor(env1$yld,env2$yld, use="complete.obs")

plot(env1$yld,env3$yld,pch=20,  
     main="Correlation of yield in E1 and E3",
     xlab="Env 1", ylab="Env 3")
cor(env1$yld,env3$yld, use="complete.obs")

plot(env1$yld,env4$yld,pch=20,
     main="Correlation of yield in E1 and E4",
     xlab="Env 1", ylab="Env 4")
cor(env1$yld,env4$yld, use="complete.obs")

plot(env2$yld,env3$yld,pch=20,
     main="Correlation of yield in  E2 and E3",
     xlab="Env 2", ylab="Env 3")
cor(env2$yld,env3$yld, use="complete.obs")

plot(env2$yld,env4$yld,pch=20,
     main="Correlation of yield in E2 and E4",
     xlab="Env 2", ylab="Env 4")
cor(env2$yld,env4$yld, use="complete.obs")

plot(env3$yld,env4$yld,pch=20,
     main="Correlation of yield in E3 and E4",
     xlab="Env 3", ylab="Env 4")
cor(env3$yld,env4$yld, use="complete.obs")


#obtain estimates of residual variance for each environment
################################
#### Environment E1  ###########
################################
attach(env1)
env1.M1<-lm(yld~Sampleid)
anova(env1.M1)
env1.M1

env1.M2 = lmer(yld ~ (1|Sampleid))
ranef(env1.M2)
anova(env1.M2)
env1.M2
detach(env1)
write.csv(ranef(env1.M2), "Genotype effect estimates for env 1.csv", row.names = F)
detach(env1)

################################
#### Environment E2  ###########
################################
attach(env2)
env2.M1<-lm(yld~Sampleid)
anova(env2.M1)
env2.M1

env2.M2 = lmer(yld ~ (1|Sampleid))
ranef(env2.M2)
anova(env2.M2)
env2.M2
write.csv(ranef(env2.M2), "Genotype effect estimates for env 2.csv", row.names = F)
detach(env2)

################################
#### Environment E3  ###########
################################
attach(env3)
env3.M1<-lm(yld~Sampleid)
anova(env3.M1)
env3.M1

env3.M2 = lmer(yld ~ (1|Sampleid))
ranef(env3.M2)
anova(env3.M2)
env3.M2
write.csv(ranef(env3.M2), "Genotype effect estimates for env 3.csv", row.names = F)
detach(env3)

################################
#### Environment E4  ###########
################################
attach(env4)
env4.M1<-lm(yld~Sampleid)
anova(env4.M1)
env4.M1

env4.M2 = lmer(yld ~ (1|Sampleid))
ranef(env4.M2)
anova(env4.M2)
env4.M2
write.csv(ranef(env4.M2), "Genotype effect estimates for env 4.csv", row.names = F)
detach(env4)



################################
#### Estimate GxE variance#####
################################
attach(ds1)
yld.M1<-lm(yld~Env+Sampleid+Sampleid*Env)
anova(yld.M1)
yld.M1

yld.M2 = lmer(yld ~ Env+(1|Sampleid)+(1|Sampleid:Env))
ranef(yld.M2)
anova(yld.M2)
yld.M2
detach(ds17)