#Program to analyze the data in Table 9.1 for Example 9.1
library(car)
yd=read.table("C:\\Research\\Book\\Data\\table91.txt")
x1 = yd[[1]]
y  = yd[[2]]
x2 = x1^2
x3 = x1^3
x4 = x1^4
x5 = x1^5
#cbind(y, x1, x2, x3, x4, x5)

#Regression analysis with y as response and x1-x5 as predictors
md = y ~ x1 + x2 + x3 + x4 + x5
lreg = lm(md)
lreg
summary.lm(lreg)
anova.lm(lreg)
vif(lreg)

# Test for the hypothesis x5=0
x5m = lm(y~x1+x2+x3+x4)
summary.lm(x5m)
anova.lm(x5m)
vif(x5m)
anova(x5m,lreg)

# Test for the hypothesis x4=0
x4m = lm(y~x1+x2+x3)
anova(x4m,x5m)

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

#Program to analyze the data in Table 9.2 for lack of fit test
yd=read.table("C:\\Research\\Book\\Data\\table92.txt")
x1 = yd[[1]]
y  = yd[[2]]
x2 = x1^2
lackfit = as.factor(x1)
lackfit

#Analysis of variance with y as response and x1 and lackfit as predictors
md = y ~ x1 + lackfit
linear = aov(md)
linear
summary(linear)
anova(linear)

md2 = y ~ x1 + x2 + lackfit
quadratic = aov(md2)
summary(quadratic)
#anova(quadratic)

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

#Program to analyze the data in Table 9.6 for Example 9.2
library(lattice)
yd=read.table("C:\\Research\\Book\\Data\\table96.txt")
x1 = yd[[2]]
x2 = yd[[3]]
y  = yd[[4]]
x11=x1^2
x12=x1*x2
x22=x2^2

#Regression analysis with y as response and x1-x22 as predictors
md = y ~ x1 + x2 + x12 + x11 + x22
lreg = lm(md)
lreg
summary(lreg)
anova(lreg)

u1=seq(80,100,0.1)
u2=seq(50,60,0.1)
u1.n=length(u1)
u2.n=length(u2)
uu=matrix(0, nrow=u1.n, ncol=u2.n)
for (i in 1:u1.n)
{
for (j in 1:u2.n)
{
cc = -5127.899+31.096*u1[i]+139.747*u2[j]-0.146*u1[j]*u2[j]
uu[i,j]=cc-0.133*(u1[i]^2)-1.144*(u2[j]^2)
}
}
uu.t = t(uu)
wireframe(uu.t)

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

#Program to analyze the data in Table 9.2 for Example 9.4
yd=read.table("C:\\Research\\Book\\Data\\table92.txt")
x1 = yd[[1]]
y  = yd[[2]]

#LOESS regression analysis with y as response and x1 as predictor
loes=function(av)
{
md = y ~ x1
lreg = loess(md, span=av, degree=1)
lreg
summary(lreg)
pred.y = predict(lreg)
resi.y = y-pred.y
cbind(x1,y,pred.y,resi.y)
}
loes4=loes(0.4); loes5=loes(0.5); loes6=loes(0.6); loes7=loes(0.7);
all=rbind(loes4,loes5,loes6,loes7)
all