#Program to analyze the data in Table 4.1 for Example 4.1
yd=read.table("C:\\Research\\Book\\Data\\table41.prn")
id = yd[[1]]
y  = yd[[2]]
x1 = yd[[3]]
x2 = yd[[4]]
dma = cbind(y, x1, x2)
cprod = t(dma)%*%dma
cprod  #cross products X'X, X'Y, Y'Y with the intercept
av.y =mean(y)
av.x1=mean(x1)
av.x2=mean(x2)
sd.y =sqrt(var(y))
sd.x1=sqrt(var(x1))
sd.x2=sqrt(var(x2))
ny  = y - av.y
nx1 = x1 - av.x1
nx2 = x2 - av.x2
ndma = cbind(ny, nx1, nx2)
ncprod=t(ndma)%*%ndma
ncprod  #corrected cross products
desc.y =summary(y)
desc.x1=summary(x1)
desc.x2=summary(x2)
rbind(desc.y, desc.x1, desc.x2)  #Descriptives statistics for y, x1 and x2
rbind(sd.y, sd.x1, sd.x2)
#Regression analysis with y as response and x1 and x2 as predictors
linear = y ~ x1 + x2
lreg = lm(linear)
#lreg
anova.lm(lreg)
summary.lm(lreg)
confint(lreg)
#fitted=predict.lm(lreg)
#cbind(dma,fitted)

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

#Program to analyze the data in Table 4.2 for Example 4.2
library(MASS)
yd=read.table("C:\\Research\\Book\\Data\\table42.prn")
id = yd[[1]]
y  = yd[[2]]
x1 = yd[[3]]
x2 = yd[[4]]
x3 = yd[[5]]
xx=cbind(1,x1,x2,x3)
#Regression analysis with y as response and x1, x2 and x3 as predictors
x23 = x2 + x3
x123 = x1 - x2 - x3
lm.1 = y ~ 1
lm.a = y ~ x1 
lm.b = y ~ x1 + x2
lm.c = y ~ x1 + x2 + x3
lm.bc= y ~ x1 + x23
lm.abc= y ~ x123 + offset(-x3)
lreg.1 = lm(lm.1)
lreg.a = lm(lm.a)
lreg.b = lm(lm.b)
anova(lreg.a)
summary.lm(lreg.a)
anova(lreg.b)
summary.lm(lreg.b)
anova(lreg.1, lreg.b)
lreg.c = lm(lm.c)
anova(lreg.c)
summary.lm(lreg.c)
anova(lreg.1, lreg.c)
lreg.bc= lm(lm.bc)
lreg.abc= lm(lm.abc)
anova.lm(lreg.bc, lreg.c)
anova.lm(lreg.abc, lreg.c)
#confint(lreg3)

#Computation of equation (4.32) in the textbook
be=lreg.c[[1]] 
rv=cbind(c(1,0))
bee=cbind(c(be[[1]],be[[2]],be[[3]],be[[4]]))
cv=rbind(c(0,0,1,-1),c(0,1,1,0)) 
inv=ginv(t(xx)%*%xx)
dv=rv - cv%*%bee 
le=inv%*%t(cv) 
ce=ginv(cv%*%le) 
b.est=bee + le%*%(ce%*%dv)
#The coefficient of  estimated restricted least squares model
b.est