#Program to analyze the data in Table 16.1 for Example 16.1
library(car)
library(nortest)
yd=read.table("C:\\Research\\Book\\Data\\table161.txt")
f0 = yd[,1]
y  = yd[,2]
xx = yd[,3]
ave= mean(xx)
x1 = xx-ave
fact=as.factor(f0)
z1=rep(0,length(y))
z2=rep(0,length(y))
z1.e=rep(0,length(y))
z2.e=rep(0,length(y))
for (i in 1:length(y))
{
if (fact[i]=='brus') z1[i]=1
if (fact[i]=='spry') z2[i]=1
if (fact[i]=='brus') z1.e[i]=1
if (fact[i]=='dipp') z1.e[i]=-1
if (fact[i]=='spry') z2.e[i]=1
if (fact[i]=='dipp') z2.e[i]=-1
}
#cbind(y,x1,z1,z2,z1.e,z2.e)
xz1 = x1*z1
xz2 = x1*z2
xz1.e = x1*z1.e
xz2.e = x1*z2.e

#Regression with reference cell coding scheme
reg = lm(y ~ z1 + z2 + x1 + xz1 + xz2)
summary(reg)
anova(reg)
linearHypothesis(reg, c("xz1=0","xz2=0"))
linearHypothesis(reg, "z1-z2=0")
linearHypothesis(reg, "z1=0")
linearHypothesis(reg, "z2=0")

reg1 = lm(y ~ z1 + z2 + x1)
summary(reg1)
anova(reg1)
vcov(reg1)
linearHypothesis(reg1, c("z1=0","z2=0"))
linearHypothesis(reg1, "z1-z2=0")
linearHypothesis(reg1, "z1=0")
linearHypothesis(reg1, "z2=0")

#Regression with effect coding scheme
reg.e = lm(y ~ z1.e + z2.e + x1 + xz1.e + xz2.e)
summary(reg.e)
anova(reg.e)
linearHypothesis(reg.e, c("xz1.e=0","xz2.e=0"))

reg.e1 = lm(y ~ z1.e + z2.e + x1)
summary(reg.e1)
anova(reg.e1)
vcov(reg.e1)
linearHypothesis(reg.e1, c("z1.e=0","z2.e=0"))
linearHypothesis(reg.e1, "z1.e-z2.e=0")
linearHypothesis(reg.e1, "2*z1.e+z2.e=0")
linearHypothesis(reg.e1, "z1.e+2*z2.e=0")

#Analysis of covariance with y=response, x=covariate, fact=factor
aov1 = aov(y~fact+x1)
aov1
anova(aov1)
Anova(aov1,type="III")
TukeyHSD(aov1,"fact")
resid=residuals(aov1)

#Checking the anova assumptions
shapiro.test(resid)
lillie.test(resid)
cvm.test(resid)
ad.test(resid)

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

#Program to analyze data in Table 16.2 for Example 16.2
library(car)
library(nortest)
yd=read.table("C:\\Research\\Book\\Data\\table162.txt")
u1 = yd[,1]
y1 = yd[,2]
u2 = yd[,3]
y2 = yd[,4]
u3 = yd[,5]
y3 = yd[,6]
y  = rbind(cbind(y1),cbind(y2),cbind(y3))
xx = rbind(cbind(u1),cbind(u2),cbind(u3))
bb = rbind(cbind(rep(1,10)),cbind(rep(2,10)),cbind(rep(3,10)))
a1 = rbind(cbind(rep(1,5)),cbind(rep(2,5)))
aa = rbind(a1,a1,a1)
ave= mean(xx)
x1 = xx-ave
fact.a=as.factor(aa)
fact.b=as.factor(bb)
fact.a=recode(fact.a,"1='A1'; 2='A2'")
fact.b=recode(fact.b,"1='B1'; 2='B2';3='B3'")

tapply(y, fact.a, mean)
tapply(y, fact.b, mean)
#Analysis of variance with y=response, fact.a and fact.b as factors
aov1 = aov(y~fact.a*fact.b)
aov1
anova(aov1)
TukeyHSD(aov1)

#Analysis of covariance with y=response, x=covariate, fact=factor
aov2 = aov(y~x1 + fact.a*fact.b)
aov2
anova(aov2)
resid=residuals(aov2)
Anova(aov2, type="III")
coef(aov2)

#Checking the anova assumptions
shapiro.test(resid)
lillie.test(resid)
cvm.test(resid)
ad.test(resid)

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

#Program to analyze data in Table 16.3 on unequal cell counts
library(car)
yd=read.table("C:\\Research\\Book\\Data\\table163.txt")
meth = yd[,1]
inst = yd[,2]
xv = yd[,3]
y  = yd[,4]
ave= mean(xv)
x1 = xv-ave
meth=as.factor(meth)
inst=as.factor(inst)
#Analysis of variance with y=response, meth and inst as factors
aov1 = aov(y~meth*inst)
aov1
anova(aov1)
TukeyHSD(aov1,"inst")

#Analysis of covariance with y=response, x1=covariate, meth and inst are factors
#Some of the results are not exactly the same as those obtained from SAS
aov2 = aov(y~x1+meth*inst)
aov2
anova(aov2)
resid=residuals(aov2)
Anova(aov2, type="III")
coef(aov2)