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

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

ds$blk<-factor(ds$blk)
ds$fam<-factor(ds$fam)
ds$dh<-factor(ds$dh)
summary(ds)

attach(ds)
fam1<-subset(ds,ds$fam=="F1")
fam2<-subset(ds,ds$fam=="F2")
fam3<-subset(ds,ds$fam=="F3")
detach(ds)

########################################################################################
# next consider conducting EDA 
# do this so that you have estimates of experimental 
# variability that can be used to make decisions about 
# the quality of the data 

# the model for ANOVA is y = mean + coop + coop(field) + g ==sample + residual
# in this case, the model represents an RCBD
########################################################################################

attach(fam1)
hist(ds$PV,col="beige",main="Frequency of Phenotypic Value in Family 1",
     xlab="Phenotypic Values")
anova(lm(PV~ blk + dh))
detach(fam1)

attach(fam2)
hist(ds$PV,col="beige",main="Frequency of Phenotypic Value in Family 2",
     xlab="Phenotypic Values")
anova(lm(PV~ blk + dh))
detach(fam2)

attach(fam3)
hist(ds$PV,col="beige",main="Frequency of Phenotypic Value in Family 3",
     xlab="Phenotypic Values")
anova(lm(PV~ blk + dh))
detach(fam3)


attach(ds)
hist(ds$PV,col="beige",main="Frequency of Phenotypic Value ",
     xlab="Phenotypic Values")
boxplot(ds$PV~ds$fam,main="Boxplot by family",
      xlab="Sample",ylab="Phenotypic Value")
anova(lm(PV~ blk + fam + fam:dh))
anova(lm(PV~blk + dh))
detach(ds)