rm(list=ls())
options(max.print=2000000)
setwd("C:/Users/wdbeavis/Desktop/PBEA/QG/Module 9")

#################################################################
##                    Installing 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)

#load in the data set from a spreadsheet (excel.csv)
ds<-read.csv("QG_Mod9_ALA9.3_ds.csv",header=T)
summary(ds)

###################################################

ds$dh<-as.character(ds$dh)
ds$env<-factor(ds$env)
ds$block<-factor(ds$block)
ds$fam<-as.character(ds$fam)

summary(ds)
###################################################

dsfamA<-subset(ds,fam=="A")
dsfamB<-subset(ds,fam=="B")

dsenv1<-subset(ds,env==1)
dsenv2<-subset(ds,env==2)

dsfamAenv1<-subset(dsfamA,env==1)
dsfamAenv2<-subset(dsfamA,env==2)

dsfamBenv1<-subset(dsfamB,env==1)
dsfamBenv2<-subset(dsfamB,env==2)
###################################################


attach(dsfamBenv2)
summary(dsfamBenv2)
hist(dsfamBenv2$Pheno,col="beige",
     main="Frequency of Phenotypic Value in Family B Environment 2",
     xlab="Phenotypic Values")
#boxplot()
anova(lm(Pheno ~ block + dh))
# REML code  
detach(dsfamBenv2)
####################################################
##repeat the above code for all combinations
##of family and environment ##
####################################################


#### analyses by environment ####
attach(dsenv1)
#histogram
#boxplot
anova(lm(Pheno ~ block + fam + fam:dh))
anova(lm(Pheno ~ fam + fam:dh))
# REML code
detach(dsenv1)

attach(dsenv2)
#histogram
#boxplot
anova(lm(Pheno ~ fam + fam:dh))
# REML code
detach(dsenv2)


#### analyses by family ####
attach(dsfamA)
#histogram
#boxplot
anova(lm(Pheno ~ env + env:block + dh))
# REML code
detach(dsfamA)

attach(dsfamB)
#histogram
#boxplot
anova(lm(Pheno ~ env + env:block + dh))
# REML code
detach(dsfamB)


attach(ds)
#histogram
#boxplot
anova(lm(Pheno ~ env + env:block + fam + fam:dh + env*fam + env*fam:dh))
anova(lm(Pheno ~ env + fam + fam:dh + env*fam ))

model <- lmer(data=ds, Pheno ~ env + (1|fam) + (1|fam:dh) + (1|env:fam) )
                                                               
anova(model)
model

detach(ds)