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

#make sure the working directory is correct
setwd("/Users/aamahama/Documents/PBEA-SIL_Instructor_Guide-Project/PBEA-SIL Project/Instructor Guides/Quantitative Genetics - Instructor Resources/Module 12 Multi Environment Trials Linear Mixed Models/ALAs")

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)

ds<-read.csv("MET ds5.csv",header=T)
summary(ds)
ds$line<-factor(ds$line)
ds$env<-factor(ds$env)
ds$rep<-factor(ds$rep)
ds$mscore<-as.numeric(ds$mscore)
summary(ds)
############################################################
### in addition to the usual exploratory data analyses #####
### previously conducted in other assignments          #####
### be sure to look at the linear relationship between #####
### the predictor variable mscore and the response variable#
############################################################
dsenv1<-subset(ds,env=="1")
dsenv2<-subset(ds,env=="2")
dsenv3<-subset(ds,env=="3")
dsenv4<-subset(ds,env=="4")
dsenv5<-subset(ds,env=="5")
dsenv6<-subset(ds,env=="6")
dsenv7<-subset(ds,env=="7")
dsenv8<-subset(ds,env=="8")
dsenv9<-subset(ds,env=="9")
dsenv10<-subset(ds,env=="10")

############################################################
# consider analysis of variance of a model in which lines
# are considered categorical variables, that is to say
# factors in a factorial treatement design
# linear model: y=mean+env+line +res: What about GxE?
# linear model: y=mean+env+line + env:line + res
# since the data are balanced the MoM can be used 
############################################################
attach(ds)
M1yld.allenv<-lm(yld~env + mscore + env:mscore)
anova(M1yld.allenv)
M1yld.allenv
detach(ds)