#Program to analyze the data for Example 14.1
library(car)
library(nortest)
yd=read.table("C:\\Research\\Book\\Data\\example141.txt")
y1 = yd[[1]]
y2 = yd[[2]]
y3 = yd[[3]]
y4 = yd[[4]]
y5 = yd[[5]]
e1= rep(0.0, length(y1))
e2= rep(0.083, length(y1))
e3= rep(0.29, length(y1))
e4= rep(0.50, length(y1))
e5= rep(0.86, length(y1))
y=rbind(cbind(y1),cbind(y2),cbind(y3),cbind(y4),cbind(y5))
fact=as.factor(rbind(cbind(e1),cbind(e2),cbind(e3),cbind(e4),cbind(e5)))
x1=rep(0,length(y))
x2=rep(0,length(y))
x3=rep(0,length(y))
x4=rep(0,length(y))
x1.e=rep(0,length(y))
x2.e=rep(0,length(y))
x3.e=rep(0,length(y))
x4.e=rep(0,length(y))
for (i in 1:length(y))
{
if (fact[i]==0.0) x1[i]=1
if (fact[i]==0.083) x2[i]=1
if (fact[i]==0.29) x3[i]=1
if (fact[i]==0.50) x4[i]=1
if (fact[i]==0.0) x1.e[i]=1
if (fact[i]==0.86) x1.e[i]=-1
if (fact[i]==0.083) x2.e[i]=1
if (fact[i]==0.86) x2.e[i]=-1
if (fact[i]==0.29) x3.e[i]=1
if (fact[i]==0.86) x3.e[i]=-1
if (fact[i]==0.50) x4.e[i]=1
if (fact[i]==0.86) x4.e[i]=-1
}
#Regression with reference cell coding scheme
reg = lm(y ~ x1 + x2 + x3 + x4)
summary(reg)
anova(reg)
linearHypothesis(reg, "x1-x2=0")
linearHypothesis(reg, "x1-x3=0")
linearHypothesis(reg, "x1=0")
linearHypothesis(reg, "x3=0")

#Regression with effect coding scheme
reg.e = lm(y ~ x1.e + x2.e + x3.e + x4.e)
summary(reg.e)
anova(reg.e)
linearHypothesis(reg.e, "x1.e-x2.e=0")
linearHypothesis(reg.e, "x1.e-x3.e=0")
linearHypothesis(reg.e, "2*x1.e+x2.e+x3.e+x4.e=0")
linearHypothesis(reg.e, "x1.e+x2.e+2*x3.e+x4.e=0")

#Analysis of variance with y as response and fact(pressure) as factor
aov1 = aov(y~fact)
aov1
anova(aov1)
resid=residuals(aov1)

#Checking the anova assumptions
lillie.test(y)
cvm.test(y)
ad.test(y)
bartlett.test(y~fact)
TukeyHSD(aov1)

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

#Program to analyze the data in Table 14.2 for Example 14.2
library(car)
yd=read.table("C:\\Research\\Book\\Data\\table142.txt")
pc = yd[,1]
y1 = yd[1,2:6]
y2 = yd[2,2:6]
y3 = yd[3,2:6]
y4 = yd[4,2:6]
y5 = yd[5,2:6]
e1= rep(15, length(y1))
e2= rep(20, length(y1))
e3= rep(25, length(y1))
e4= rep(30, length(y1))
e5= rep(35, length(y1))
y=rbind(t(y1),t(y2),t(y3),t(y4),t(y5))
fact=as.factor(rbind(cbind(e1),cbind(e2),cbind(e3),cbind(e4),cbind(e5)))

#Analysis of variance with y as response and fact(cotton %) as factor
aov1 = aov(y~fact)
aov1
summary(aov1)

linearHypothesis(aov1,c(0,1,0,-1,-2)) #Linear effect
linearHypothesis(aov1,c(0,1,2,1,-2)) #Quadratic effect
linearHypothesis(aov1,c(0,-2,0,2,-1)) #Cubic effect
linearHypothesis(aov1,c(0,4,-6,4,-1)) #Quartic effect

#Linear model by using the powers of the variable 'cotton %'
x1=rbind(cbind(e1),cbind(e2),cbind(e3),cbind(e4),cbind(e5))
x2=x1^2; x3=x1^3; x4=x1^4;
md = y ~ x1+x2+x3+x4
lreg=lm(md)
anova(lreg)
summary(lreg)

linearHypothesis(lreg,"x4=0")
linearHypothesis(lreg,c("x3=0","x4=0"))

md1 = y ~ x1+x2+x3
lreg1=lm(md1)
anova(lreg1)
summary(lreg1)