rm(list=ls())
options(max.print=2000000)

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

setwd("C:/Users/wdbeavis/Desktop/PBEA/QG/Module 11")

ds<-read.csv("QG_Mod11_ALA11.2.csv",header=T)
summary(ds)
ds$env<-factor(ds$env)
ds$blk<-factor(ds$blk)
ds$hyb<-factor(ds$hyb)
summary(ds)


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

dsenvA<-subset(ds,env=="A")
dsenvB<-subset(ds,env=="B")
dsenvC<-subset(ds,env=="C")
dsenvD<-subset(ds,env=="D") 

##########################
#### Environment A  ######
##########################
hist(dsenvA$pheno.t1,col="beige",
     main="Frequency of Pheno.t1in Env A",
     xlab="Phenotype t1 Values")

hist(dsenvA$pheno.t2,col="beige",
     main="Frequency of Pheno.t2in Env A",
     xlab="Phenotype t2 Values")

summary(aov(data=dsenvA, pheno.t1 ~  hyb))
summary(aov(data=dsenvA, pheno.t2 ~  hyb))

plot(dsenvA$pheno.t2,dsenvA$pheno.t1,pch=20, xlab="t2 in Env A", ylab="t1 in Env A")
cor(dsenvA$pheno.t1,dsenvA$pheno.t2, use="complete.obs")

##########################
#### Environment B  ######
##########################
hist(dsenvB$pheno.t1,col="beige",
     main="Frequency of Pheno.t1 in Env B",
     xlab="Phenotype t1 Values")

hist(dsenvB$pheno.t2,col="beige",
     main="Frequency of Pheno.t2 in Env B",
     xlab="Phenotype t2 Values")

summary(aov(data=dsenvB, pheno.t1 ~  hyb))
summary(aov(data=dsenvB, pheno.t2 ~  hyb))

plot(dsenvB$pheno.t2,dsenvB$pheno.t1,pch=20, xlab="t2 in Env B", ylab="t1 in Env B")
cor(dsenvB$pheno.t1,dsenvB$pheno.t2, use="complete.obs")

##########################
#### Environment C  ######
##########################
hist(dsenvC$pheno.t1,col="beige",
     main="Frequency of Pheno.t1 in Env C",
     xlab="Phenotype t1 Values")

hist(dsenvC$pheno.t2,col="beige",
     main="Frequency of Pheno.t2 in Env C",
     xlab="Phenotype t2 Values")

summary(aov(data=dsenvC, pheno.t1 ~  hyb))
summary(aov(data=dsenvC, pheno.t2 ~  hyb))

plot(dsenvC$pheno.t2,dsenvC$pheno.t1,pch=20, xlab="t2 in Env C", ylab="t1 in Env C")
cor(dsenvC$pheno.t1,dsenvC$pheno.t2, use="complete.obs")

##########################
#### Environment D  ######
##########################
hist(dsenvD$pheno.t1,col="beige",
     main="Frequency of Pheno.t1 in Env D",
     xlab="Phenotypic Values")

hist(dsenvD$pheno.t2,col="beige",
     main="Frequency of Pheno.t2 in Env D",
     xlab="Phenotypic Values")

summary(aov(data=dsenvD, pheno.t1 ~  hyb))
summary(aov(data=dsenvD, pheno.t2 ~  hyb))

plot(dsenvD$pheno.t2,dsenvA$pheno.t1,pch=20, xlab="t2 in Env D", ylab="t1 in Env D")
cor(dsenvD$pheno.t1,dsenvA$pheno.t2, use="complete.obs")

#########################################
#### Combined across Environments  ######
#########################################
boxplot(ds$pheno.t1~ds$env,main="Boxplot of t1 by Environment",
        xlab="Environment",ylab="t1")
boxplot(ds$pheno.t2~ds$env,main="Boxplot of t2 by Environment",
        xlab="Environment",ylab="t2")

plot(ds$pheno.t2,ds$pheno.t1,pch=20, xlab="t2", ylab="t1")
cor(ds$pheno.t1,ds$pheno.t2, use="complete.obs")


######################################################
### Conduct Data Analyses on data combined across ####
### environments to obtain an estimate of V(GxE)  ####
######################################################

anova(lm(data=ds,pheno.t1~env+ blk/env + hyb + hyb*env))

anova(lm(data=ds,pheno.t2~env+ blk/env + hyb + hyb*env))


attach(ds)
outt1<-lsmeans(lm(pheno.t1 ~ env + blk + hyb), ~ as.factor(hyb))
outt2<-lsmeans(lm(pheno.t2 ~ env + blk + hyb), ~ as.factor(hyb))
write.csv(outt1, file= "QG_Mod11_lsmeanst1.csv")
write.csv(outt2, file= "QG_Mod11_lsmeanst2.csv")
detach(ds)

dslsmt1<-read.csv("QG_Mod11_lsmeanst1.csv",header=T)
dslsmt1$lsmeant1<-dslsmt1$lsmean
summary(dslsmt1)
dslsmt1<-dslsmt1[-c(1,3:7)]
summary(dslsmt1)

dslsmt2<-read.csv("QG_Mod11_lsmeanst2.csv",header=T)
dslsmt2$lsmeant2<-dslsmt2$lsmean
summary(dslsmt2)
dslsmt2<-dslsmt2[-c(1,3:7)]
summary(dslsmt2)

combine<-merge(dslsmt1,dslsmt2,by= c("hyb"),all=TRUE)

plot(combine$lsmeant2,combine$lsmeant1,pch=20,xlab="lsmean for t2", ylab="lsmean for t1")
cor(combine$lsmeant2,combine$lsmeant1)

write.csv(combine, file= "QG_Mod11_lsmeanst1andt2.csv")