#Program to analyze the data in Table 5.1 for Example 5.1
yd=read.table("C:\\Research\\Book\\Data\\table51.txt")
id  = yd[[1]]
x = yd[[2]]
y = yd[[3]]
#Plot of the scatter diagram
plot(x, y)

#Regression analysis with y as response and x as predictor
simple = y ~ x
lreg = lm(simple)
lreg
summary.lm(lreg)
anova.lm(lreg)
predicted=predict.lm(lreg)
residuals=residuals(lreg)
resid.s0=rep(0,length(y))
#Plot of residuals versus the predicted
plot(predicted, residuals)
lines(predicted, resid.s0)

x2=x*x
lm2 = y ~ x + x2
lreg2=lm(lm2)
lreg2
summary.lm(lreg2)
anova.lm(lreg2)
predicted2=predict.lm(lreg2)
residuals2=residuals(lreg2)
#Plot of residuals versus the predicted in the quadratic model
plot(predicted2, residuals2)
lines(predicted2, resid.s0)

linear.hypothesis(lreg2, "x2=0")

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

#Program to analyze data in Table 5.2 for Example 5.2
yd=read.table("C:\\Research\\Book\\Data\\table52.txt")
id= yd[[1]]
x = yd[[2]]
y = yd[[3]]

av.x=mean(x)
av.y=mean(y)
sd.x=sqrt(var(x))
sd.y=sqrt(var(y))
desc.x=cbind(rbind(summary(x)),sd.x)
desc.y=cbind(rbind(summary(y)),sd.y)
rbind(desc.x, desc.y)  #Descriptives statistics for x and y

#Regression analysis with y as response and x as predictor
md1 = y ~ x
lreg = lm(md1)
lreg
summary.lm(lreg)
anova.lm(lreg)
predicted=predict.lm(lreg)
residuals=residuals(lreg)
resid.s0=rep(0,length(y))
#Plot of residuals versus the predicted
plot(predicted, residuals)
lines(predicted, resid.s0)

shapiro.test(residuals)
qqnorm(residuals)
qqline(residuals)

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

#Program to analyze the data in Table 5.4 for Example 5.3
yd=read.table("C:\\Research\\Book\\Data\\table54.txt")
age = yd[[1]]
sbp = yd[[2]]

#Regression analysis with sbp as response and age as predictor
md1 = sbp ~ age
lreg = lm(md1)
lreg
summary.lm(lreg)
anova.lm(lreg)
predicted=predict.lm(lreg)
residuals=residuals(lreg)
ri=rstandard(lreg)
influ=influence(lreg)
leverage=influ[[1]]
sigma=summary(lreg)$sigma
di=residuals/sigma
pre=residuals/(1-leverage)
result1=cbind(age, sbp, predicted, residuals, di, ri, pre, leverage)
result1=round(result1, digits=4)
result1

resid.s2n=rep(-2,length(sbp))
resid.s0=rep(0,length(sbp))
resid.s2p=rep(2,length(sbp))
#Plot of residuals versus the predicted
plot(predicted, ri, ylim=c(-5, 2.5))
lines(predicted, resid.s2n)
lines(predicted, resid.s0)
lines(predicted, resid.s2p)

RStudent=rstudent(lreg)
dffits=dffits(lreg)
dfbetas=dfbetas(lreg)
covrat=covratio(lreg)
result2=cbind(residuals, RStudent, leverage, covrat, dffits, dfbetas)
result2=round(result2, digits=4)
result2

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

#Program to analyze the data on page 76 for Example 5.4
yd=read.table("C:\\Research\\Book\\Data\\example54.txt")
y = yd[[1]]
x1= yd[[2]]
x2= yd[[3]]
x3= yd[[4]]

#Regression analysis with y as response and x1, x2, x3 as predictors
md3 = y ~ x1+x2+x3
lreg3 = lm(md3)
lreg3
summary.lm(lreg3)
anova.lm(lreg3)
md2 = y ~ x1+x2
lreg = lm(md2)
lreg
summary.lm(lreg)
anova.lm(lreg)
pred=predict.lm(lreg)
resid=residuals(lreg)
stud=rstandard(lreg)
influ=influence(lreg)
rstud=rstudent(lreg)
lev=influ[[1]]
cooks=cooks.distance(lreg)
dffits=dffits(lreg)
#sigma=summary(lreg)$sigma
#di=residuals/sigma
result=cbind(y, pred, resid, stud, rstud, lev, cooks, dffits)
result=round(result, digits=4)
result

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

#Program to analyze the data in Table 5.6 for Example 5.5
library(Hmisc)
library(car)
yd=read.table("C:\\Research\\Book\\Data\\table56.txt")
id= yd[[1]]
x1= yd[[2]]
x2= yd[[3]]
x3= yd[[4]]
x4= yd[[5]]
y = yd[[6]]
xv=cbind(x1,x2,x3,x4)
cor(xv)
rcorr(xv)

#Regression analysis with y as response and x1, x2, x3, x4 as predictors
md = y ~ x1+x2+x3+x4
lreg = lm(md)
lreg
summary.lm(lreg)
anova.lm(lreg)
vif=vif(lreg)
tol=1/vif
result=cbind(tol,vif)
result

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

#Program to analyze data in Table 5.2 for Example 5.2 [page 84]
library(MASS)
yd=read.table("C:\\Research\\Book\\Data\\table52.txt")
id= yd[[1]]
x = yd[[2]]
y = yd[[3]]

md1 = y ~ x
#Create the indicator variable for the variable x
group=rep(0,length(x))
wt=rep(0,length(x))
for (i in 1:length(x))
{
if (20<=x[i] & x[i]<30) 
{group[i]=1 
wt[i]=1/(25^2)}
if (30<=x[i] & x[i]<40) 
{group[i]=2 
wt[i]=1/(35^2)}
if (40<=x[i] & x[i]<50) 
{group[i]=3 
wt[i]=1/(45^2)}
if (50<=x[i] & x[i]<60) 
{group[i]=4 
wt[i]=1/(55^2)}
}
group=factor(group)
nyd=cbind(yd,group)
x1 = nyd[[2]]
y1 = nyd[[3]]
new=function(ud)
{
x1 = ud[[2]]
y1 = ud[[3]]
lreg1=lm(y1~x1)
resid1=residuals(lreg1)
size=length(x1)
aver=mean(x1)
vari=var(resid1)
col5=vari/aver
col6=vari/(aver^2)
col7=vari/sqrt(aver)
res=cbind(size, aver, vari, col5, col6, col7)
}
u1=split(nyd,group)$'1'
res1=new(u1)
u2=split(nyd,group)$'2'
res2=new(u2)
u3=split(nyd,group)$'3'
res3=new(u3)
u4=split(nyd,group)$'4'
res4=new(u4)
gr=seq(1:4)
ress = cbind(gr,rbind(res1, res2, res3, res4))
ress

lreg.wt = lm(md1, weights=wt)
summary.lm(lreg.wt)
anova.lm(lreg.wt)
predict.wt=predict.lm(lreg.wt)
resid.wt=residuals(lreg.wt)
resid.wt2 = resid.wt*sqrt(wt)
shapiro.test(resid.wt2)

#Plot of residuals versus the predicted in the weighted model
resid.s0=rep(0,length(y))
plot(predict.wt, resid.wt2)
lines(predict.wt, resid.s0)

wt1=1/(x^2)
lreg.wt1=lm(md1, weights=wt1)
summary.lm(lreg.wt1)
anova.lm(lreg.wt1)

#Box-Cox transformation on page 88
cxy=boxcox(y~x, lambda=seq(-3,3,0.05))
cx=cxy[[1]]; cy=cxy[[2]]
ld=cbind(cx,cy)
m.cy=max(cy)
m.cy
for (i in 1:length(cy))
{
if (cy[i]==max(cy)) max.lda=cx[i]
}
max.lda
new.y = 1/y
new.m = new.y~x
new.reg = lm(new.m)
summary.lm(new.reg)
anova.lm(new.reg)
predict.new=predict.lm(new.reg)
resid.new=residuals(new.reg)
#Plot of residuals versus the predicted in the weighted model
plot(x, resid.new)
lines(x, resid.s0)

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

#Program to analyze the data in Table 5.9 for Example 5.6
library(MASS)
yd=read.table("C:\\Research\\Book\\Data\\table59.txt")
car=yd[[1]]
y = yd[[2]]
x = yd[[3]]

#Scatterplot of y against x
plot(x, y)

#Regression analysis with y as response and x as predictor
md1 = y ~ x
lreg = lm(md1)
lreg
summary.lm(lreg)
anova.lm(lreg)
predicted=predict.lm(lreg)
residuals=residuals(lreg)
resid.s0=rep(0,length(y))
#Plot of residuals versus the predicted
plot(predicted, residuals)
lines(predicted, resid.s0)

#Box-Cox transformation on page 93
cxy=boxcox(y~x, lambda=seq(-2,2,0.1))
cx=cxy[[1]]; cy=cxy[[2]]
m.cy=max(cy)
m.cy
for (i in 1:length(cy))
{
if (cy[i]==max(cy)) max.lda=cx[i]
}
max.lda

new1.y = 1/sqrt(y); new2.y=log(y); new3.y=y^(0.3); new4.y=sqrt(y)
n.1=new1.y~x; n.2=new2.y~x; n.3=new3.y~x; n.4=new4.y~x;
n1.reg = lm(n.1); anova.lm(n1.reg)
n2.reg = lm(n.2); anova.lm(n2.reg)
n3.reg = lm(n.3); anova.lm(n3.reg)
n4.reg = lm(n.4); anova.lm(n4.reg)

predict.n2=predict.lm(n2.reg)
resid.n2=residuals(n2.reg)
#Plot of residuals versus the predicted in the model with lambda=0
plot(predict.n2, resid.n2)
lines(predict.n2, resid.s0)

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

#Program to analyze loss data on page 101 for Robust regression
library(robustbase)
yd=read.table("C:\\Research\\Book\\Data\\loss.txt")
y=yd[[1]]
x1=yd[[2]]
x2=yd[[3]]
x3=yd[[4]]
md=y ~ x1+x2+x3
lreg = lm(md)
summary.lm(lreg)
anova.lm(lreg)
resid=residuals(lreg)
resid; svar=abs(resid)
my=cbind(y,x1,x2,x3,svar)

#We remove observations 1, 3, 4, 21 by using the value of residual
new.yd=my[which(my[,5]<3.2), ]
new.yd
y=new.yd[,1]
nx1=new.yd[,2]
nx2=new.yd[,3]
nx3=new.yd[,4]
md1=y ~ nx1+nx2+nx3
lreg = lm(md1)
summary.lm(lreg)
anova.lm(lreg)

#Robust regression using "LTS" method.
y=yd[[1]]
x1=yd[[2]]
x2=yd[[3]]
x3=yd[[4]]
md=y ~ x1+x2+x3
lts=ltsReg(md)
summary(lts)
fits=fitted(lts); resid=residuals(lts)
result1=cbind(y,fits,resid)
result1

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

#Program to analyze the data in Table 5.10 for Example 5.7
library(quantreg)
yd=read.table("C:\\Research\\Book\\Data\\table510.txt")
y = yd[[1]]
x = yd[[2]]
x2=x^2

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

qreg50 = rq(md, tau=0.5)
qreg50
#summary(qreg50)
pred50=predict(qreg50)
resid50=residuals(qreg50)
#cbind(y, x, x2, pred50, resid50)

qreg05 = rq(md, tau=0.05)
pred05=predict(qreg05)

qreg95 = rq(md, tau=0.95)
pred95=predict(qreg95)

cbind(x, y, pred05, pred50, pred95)

xt=(x-1780)/10
plot(xt,pred05,type="l", ylab="Population", xlab="X-Transformed Years")
lines(xt,pred50)
lines(xt,pred95)

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

#Program to analyze the data in Table 5.11 for Example 5.8
library(quantreg)
yd=read.table("C:\\Research\\Book\\Data\\table511.txt")
y = yd[[1]]
x = yd[[2]]

#Regression analysis with y as response and x as predictor
md = y ~ x 
lreg = lm(md)
lreg
summary.lm(lreg)
anova.lm(lreg)

qreg50 = rq(md, tau=0.5)
summary(qreg50)
pred50=predict(qreg50)
resid50=residuals(qreg50)
cbind(y, pred50, resid50);

resid.s0=rep(0,length(x))
#Plot of residuals versus the predicted
plot(x, resid50)
lines(x, resid.s0)

###############################################################
#Program to analyze data in Table 5.12 on page 113 [Autocorrelation]
library(car)
library(bstats)
yd=read.table("C:\\Research\\Book\\Data\\table512.txt")
year = yd[[1]]
y = yd[[2]]
t = year-1973

#Regression analysis with y as response and t as predictor
md = y ~ t
lreg = lm(md)
summary.lm(lreg)
anova.lm(lreg)
pred=predict.lm(lreg)
resid=residuals(lreg)
resid.s0=rep(0,length(y))
#Plot of residuals versus the predicted
plot(t, resid)
lines(t,resid)
lines(t, resid.s0)

durbinWatsonTest(lreg, max.lag=1, alternative="two.sided")
dw.test(md, alternative="two.sided")

#ty=as.ts(y)
arima(y, order=c(1,0,0), method="ML", xreg=t)
resid=arima(y, order=c(1,0,0), method="ML", xreg=t)$residuals
resid

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

#Program to analyze the data in Table 5.6 for Example 5.9 [Ridge Regression]
library(MASS)
yd=read.table("C:\\Research\\Lawal_Book\\Data\\table56.txt")
id = yd[[1]]
x1 = yd[[2]]
x2 = yd[[3]]
x3 = yd[[4]]
x4 = yd[[5]]
y =  yd[[6]]
lda=seq(0, .4, .01)
lda
#Regression analysis with y as response and x1-x4 as predictors
md = y ~ x1 + x2 + x3 + x4
lreg = lm(md)
summary.lm(lreg)
anova.lm(lreg)

#Ridge regression for the above linear model
lm.ridge(md, lambda=lda)