#Program to analyze the data in Table 8.1 for Example 8.1
library(car)
yd=read.table("C:\\Research\\Book\\Data\\table81.txt")
y  = yd[[2]]
x1 = yd[[3]]
status = yd[[4]]
x2=2-as.numeric(status)
x3 = x1*x2
#cbind(y, x1, x2, x3)

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

#Function to test for the parallel hypothesis x3=0
linearHypothesis(lreg, "x3=0")
#Function to test for the co-incidence hypothesis x2=0,x3=0
linearHypothesis(lreg, c("x2=0","x3=0"))
#Function to test for the common intercept hypothesis x2=0
linearHypothesis(lreg, "x2=0")

#Alternative way to test for the above by using anova function
para = lm(y~x1+x2)
coin = lm(y~x1)
comm = lm(y~x1+x3)

#Function to test for the parallel hypothesis x3=0
anova(para,lreg)
#Function to test for the co-incidence hypothesis x2=0,x3=0
anova(coin,lreg)
#Function to test for the common intercept hypothesis x2=0
anova(comm,lreg)

#For effect coding, we want "N"=1 and "S"=-1. Currently, "N"=1 in x2
#"S"=0 in x2. We need to define a new variable x22 to reflect this.
x22 = 2*x2-1
x33 = x1*x22
cbind(y, x22, x33)
# Above variables can be used as before to fit models and carry out tests

#Regression analysis with y as response and x1, x22, x33 as predictors
md1 = y ~ x1 + x22 + x33
lreg1 = lm(md1)
summary.lm(lreg1)
anova.lm(lreg1)

#Function to test for the parallel hypothesis x33=0
linearHypothesis(lreg1, "x33=0")
#Function to test for the co-incidence hypothesis x22=0,x33=0
linearHypothesis(lreg1, c("x22=0","x33=0"))
#Function to test for the common intercept hypothesis x22=0
linearHypothesis(lreg1, "x22=0")

#Regression analysis with y as response and x1, x22 as predictors
md2 = y ~ x1 + x22
lreg2 = lm(md2)
summary.lm(lreg2)
anova.lm(lreg2)
#Function to test for the common intercept hypothesis x2=0
linearHypothesis(lreg2, "x22=0")

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

#Program to analyze the data in Table 8.2 for Example 8.2
yd=read.table("C:\\Research\\Book\\Data\\table82.txt")
y  = yd[[1]]
x1 = yd[[2]]
state = yd[[3]]
no=length(y)
x2=rep(0,no)
x3=rep(0,no)
for (i in 1:no)
{
if (state[i]=="A") x2[i]=1
if (state[i]=="B") x3[i]=1
}
x4 = x1*x2
x5 = x1*x3

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

#Models to test for the different hypotheses
para = lm(y~x1+x2+x3)
coin = lm(y~x1)
comm = lm(y~x1+x4+x5)

#Function to test for the parallel hypothesis x4=0,x5=0
anova(para,lreg)
#Function to test for the co-incidence hypothesis x2=0,x3=0,x4=0,x5=0
anova(coin,lreg)
#Function to test for the common intercept hypothesis x2=0,x3=0
anova(comm,lreg)

#For effect coding, we use the following transformation
x22=rep(0,no)
x33=rep(0,no)
for (i in 1:no)
{
if (state[i]=="A") x22[i]=1
if (state[i]=="C") x22[i]=-1
if (state[i]=="B") x33[i]=1
if (state[i]=="C") x33[i]=-1
}
x44 = x1*x22
x55 = x1*x33
# Above variables can be used as before to fit models and carry out tests
# We include the full model here and no tests
md1 = y ~ x1 + x22 + x33 + x44 + x55
lreg1 = lm(md1)
lreg1
summary.lm(lreg1)
anova.lm(lreg1)

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

#Program to analyze the data in Table 8.3 for Example 8.3
yd=read.table("C:\\Research\\Book\\Data\\table83.txt")
y  = yd[[2]]
x1 = yd[[3]]
no=length(y)
x2=rep(0,no)
for (i in 1:no)
{
if (x1[i]>=50) x2[i]=1
}
x3 = (x1-50)*x2
#cbind(y,x1,x2,x3)

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