BasicASReml

# set working directory to "course material"
setwd("H:\CORP\ITCRD-CORPRND\Agri_Common\GenStat_ASReml_VSNi\ICT DAY2\ICT DAY2 PART1\MET")


library(asreml)
#library(myf)

# read data
sbt <- asreml.read.table("sbt_all.dat",sep="",header=TRUE,na.strings="*")

# sort data frame
sbt <- sbt[order(sbt$Trial,sbt$Row,sbt$Col),]

# get cycle factor for T3
sbt$colcycle <- as.factor(((as.numeric(sbt$Col)-1)%%4))
# get cycle factor for T1
sbt$rowcycle <- as.factor(((as.numeric(sbt$Row)-1)%%2))

# calculate variate of interest
sbt$dmyield <- sbt$yield*sbt$dmpc/100

# get subsets for individual trials
sbt.T1 <- subset(sbt,Trial=="T1")
sbt.T2 <- subset(sbt,Trial=="T2")
sbt.T3 <- subset(sbt,Trial=="T3")
sbt.T4 <- subset(sbt,Trial=="T4")
sbt.T5 <- subset(sbt,Trial=="T5")

# establish models for individual trials 

# T1
T1.asr <- asreml(fixed=dmyield~1+rowcycle, random=~Variety+Rep+Col+units,
                  rcov=~ar1(Row):ar1(Col), na.method.X="include", data=sbt.T1)
summary(T1.asr)$varcomp
metplot(T1.asr)
plot(T1.asr)
anova(T1.asr)

# T2
T2.asr <- asreml(fixed=dmyield~1, random=~Variety+Rep+Col+units,
                 rcov=~ar1(Row):ar1(Col), na.method.X="include", data=sbt.T2)
T2.asr <- update(T2.asr)
summary(T2.asr)$varcomp
plot(T2.asr)
metplot(T2.asr)

# T3
T3.asr <- asreml(fixed=dmyield~1+colcycle, random=~Variety+Rep+Col+units,
                 rcov=~ar1(Row):ar1(Col), na.method.X="include", data=sbt.T3)
summary(T3.asr)$varcomp
plot(T3.asr)
metplot(T3.asr)

# T4
T4.asr <- asreml(fixed=dmyield~1, random=~Variety+Rep,
                 rcov=~ar1(Row):ar1(Col), na.method.X="include", data=sbt.T4)
summary(T4.asr)$varcomp
plot(T4.asr)
metplot(T4.asr)

# T5
T5.asr <- asreml(fixed=dmyield~1+lin(Row), random=~Variety+Rep,
                 rcov=~id(Row):ar1(Col), na.method.X="include", data=sbt.T5)
summary(T5.asr)$varcomp
plot(T5.asr)
#metplot(T5.asr)

# one-stage joint analysis: independence across trials for GxE

diag.asr <- asreml(fixed=dmyield~1+Trial+at(Trial,1):rowcycle+
                         at(Trial,3):colcycle+at(Trial,5):lin(Row),
                   random=~diag(Trial):Variety + at(Trial):Rep
                          + at(Trial,c(1:3)):Col+ at(Trial,c(1:3)):units,
                   rcov=~at(Trial,c(1:4)):ar1(Row):ar1(Col)
                        +at(Trial,5):id(Row):ar1(Col),
                   na.method.X="include", data=sbt)
# use update function to get convergence
diag.asr <- update(diag.asr)
summary(diag.asr)$varcomp