#Program to analyze the data for Example 15.3
library(car)
yd=read.table("C:\\Research\\Book\\Data\\example153.txt")
y1 = yd[,1]; y2=yd[,2]; y3=yd[,3]; y4=yd[,4]
y = rbind(cbind(y1),cbind(y2),cbind(y3),cbind(y4))
x1 = rep(1:4,each=15)
x2 = rep(rep(1:3,each=5), 4)

#cbind(y,x1,x2)
loca=factor(x1)
spec=factor(x2)
spec=recode(spec,"1='A';2='B';3='C'")

md=y~spec*loca
anov = aov(md)
summary(anov)
summary.lm(anov)

ma1=cbind(y,loca); ma2=cbind(y,spec); ma3=cbind(y,loca,spec)
by(ma1, INDICES=loca, FUN=mean)
by(ma2, INDICES=spec, FUN=mean)
by(ma3, INDICES=spec:loca, FUN=mean)

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

#Program to analyze the data in Table 15.6 for Example 15.4
library(car)
yd=read.table("C:\\Research\\Book\\Data\\table156.txt")
y1 = yd[,1]; y2=yd[,2]; y3=yd[,3]; y4=yd[,4]
y = rbind(cbind(y1),cbind(y2),cbind(y3),cbind(y4))
x1 = rep(1:4,each=6)
x2 = rep(seq(1:3), 8)
x3 = rep(rep(1:2,each=3), 4)
#cbind(y,x1,x2,x3)

fac.a=factor(x1)
fac.b=factor(x2)
repl =factor(x3)
fac.a=recode(fac.a,"1='A1';2='A2';3='A3';4='A4'")
fac.b=recode(fac.b,"1='B1';2='B2';3='B3'")

md=y~repl+fac.a+fac.b+fac.a:fac.b
anov = aov(md)
summary(anov)
summary.lm(anov)

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

#Program to analyze the data in Table 15.12 for Example 15.5
library(car)
yd=read.table("C:\\Research\\Book\\Data\\table1512.txt")
y1 = yd[,1]
y  = cbind(y1)
x1 = rep(1:3,each=16)
x2 = c(rep(1:4,each=4),rep(1:4,each=4),rep(1:4,each=4))
x3 = rep(seq(1:4),12)
no=length(y)
#cbind(y,x1,x2,x3)
time=factor(x1)
variety=factor(x2)
block=factor(x3)
time=recode(time,"1='Feb';2='Mar';3='Apr'")
variety=recode(variety,"1='V';2='R';3='F';4='G'")
block=recode(block,"1='I';2='II';3='III';4='IV'")

#Analysis of variance with y as response
aov1 = aov(y~block + time + variety + variety:time)
aov1
anova(aov1)
TukeyHSD(aov1,"time")
TukeyHSD(aov1,"variety")
coef(aov1)
Anova(aov1, type="II")

linearHypothesis(aov1,c(0,0,0,0,-1,-1,0,0,0,0,0,0,0,-1,-1)) #contrast one
linearHypothesis(aov1,c(0,0,0,0,-1,-1,0,0,0,0,0,-1,-1,0,0)) #contrast two

ma1=cbind(y,time); ma2=cbind(y,variety); ma3=cbind(y,time,variety)
by(ma1, INDICES=time, FUN=mean)
by(ma2, INDICES=variety, FUN=mean)
by(ma3, INDICES=time:variety, FUN=mean)