###################################################
### chunk number 1: 
###################################################
options(show.signif.stars=F)
options(width=75)


###################################################
### chunk number 2: 
###################################################
#read in the data from a csv file
dat1=read.table("Data Sets/Chapter 3/Examples/ex-3-1.csv", header=T, sep=",")
str(dat1)  #checking data format
TIM=dat1[,2] #assign the second column in dat1 to be variable TIM 
CAS=dat1[,3] #assign the third column in dat1 to be variable CAS 
DIS=dat1[,4] #assign the forth column to be variable DIS
EX3=data.frame(TIM, CAS, DIS) 
head(EX3)  #printout the first few row in dat1
cor(EX3)


###################################################
### chunk number 3: 
###################################################
pairs(EX3)


###################################################
### chunk number 4: 
###################################################
library(xtable)
g=lm(TIM~CAS+DIS, data=EX3)
xtable(anova(g))
xtable(summary(g))
yhat=fitted(g)  #fitted values
ei=resid(g)     #residuals
ri=rstandard(g) #studentized residuals: ei/(s*sqrt(1-hii))
ti=rstudent(g)  #external studentized residuals: ei/(s(i)*sqrt(1-hii))
hii=hatvalues(g)    #hii: diagnoal of H
out=cbind(TIM, yhat, ei, ri, ti, hii)
xtable(out)  #two decimal places for printout


###################################################
### chunk number 5: 
###################################################
par(mfrow=c(2,2),mar=c(3,4,2,2))
plot(yhat, ei, main="Residuals vs fitted")
plot(yhat, ri, main="Studentized Res vs fitted")
plot(yhat, ti, main="Exernal Studentized Res vs fitted")
plot(yhat, TIM, main="Yobs ei vs fitted")


###################################################
### chunk number 6: 
###################################################
par(mfrow=c(1,1),mar=c(5,4,2,2))
qqnorm(ri);qqline(ri,col=2)


###################################################
### chunk number 7: 
###################################################
my.ppnorm=function(r) {
  n=length(r)
  order.res=sort(r)
  i=1:n
  pi=pnorm((i-1/2)/n)
  plot(order.res, pi, xlab = "Sample Residuals", ylab = "Theoretical Probability", pch=16)
}
my.ppnorm(ei)


###################################################
### chunk number 8: 
###################################################
par(mfrow=c(1,1))
curve(dnorm(x), -4,4)
curve(dcauchy(x), -4,4, add=TRUE, col=2) 
curve(dt(x,2), -4,4, add=TRUE ,col=3)
curve(dt(x,4), -4,4, add=TRUE ,col=3, lty=3)
curve(dexp(x),0,4, add=TRUE,col=4)
legend("topright", legend=c("Normal","Cauchy", "t","Exp"), lty=1, col=c(1:4))


###################################################
### chunk number 9: 
###################################################
par(mfrow=c(3,3), mar=c(2,2,1,1))
for (i in 1:9) { e=3+rnorm(50); qqnorm(e, main="");qqline(e,col=2) }


###################################################
### chunk number 10: 
###################################################
par(mfrow=c(3,3), mar=c(2,2,1,1))
for (i in 1:9) { e=exp(rnorm(50)); qqnorm(e, main="");qqline(e,col=2) }


###################################################
### chunk number 11: 
###################################################
par(mfrow=c(3,3), mar=c(2,2,1,1))
for (i in 1:9) { e=3+rcauchy(50); qqnorm(e, main="");qqline(e,col=2) }


###################################################
### chunk number 12: 
###################################################
par(mfrow=c(3,3), mar=c(2,2,1,1))
for (i in 1:9) { e=runif(50,-1,1); qqnorm(e, main="");qqline(e,col=2) }


###################################################
### chunk number 13: 
###################################################
e=resid(lm(TIM~DIS))
d=resid(lm(CAS~DIS))
plot(d, e, xlab="CAS", ylab="TIM Residuals")
abline(0, coef(g)['CAS'])


###################################################
### chunk number 14: 
###################################################
coef(lm(e~d))
coef(g)


###################################################
### chunk number 15: 
###################################################
library(car)    #include car library 
av.plot(g, CAS, identify.points=F)


###################################################
### chunk number 16: 
###################################################
b=coef(g)['CAS']
plot(CAS, resid(g)+b*CAS, xlab="CAS", ylab="RES+COMPONENT")


###################################################
### chunk number 17: 
###################################################
cr.plot(g, CAS)


###################################################
### chunk number 18: 
###################################################
x1=runif(20,-1,1)
x2=runif(20,-1,1)
x3=runif(20,-1,1)
y=2+x2+x3+I(x1^2)+rnorm(20, sd=0.1)
cor(cbind(y, x1, x2, x3))


###################################################
### chunk number 19: 
###################################################
g=lm(y~x1+x2+x3)
summary(g)
anova(g)


###################################################
### chunk number 20: 
###################################################
par(mfrow=c(1,3),mar=c(3,4,1,2))
av.plot(g, x1, identify.points=FALSE)  
av.plot(g, x2, identify.points=FALSE)
av.plot(g, x3, identify.points=FALSE)


###################################################
### chunk number 21: 
###################################################
par(mfrow=c(1,3),mar=c(3,4,1,2))
cr.plot(g, x1) #partial regression plot
cr.plot(g, x2) 
cr.plot(g, x3) 


###################################################
### chunk number 22: 
###################################################
sum(ei^2)   #Residual Sum of squares: SSres
press.ei=ei/(1-hii)
sum(press.ei^2) #PRESS
n=length(TIM)  #number of obs
(SSTO=(n-1)*var(TIM))  #SSTO=Syy 
1-sum(press.ei^2)/SSTO  #prediction R2


###################################################
### chunk number 23: 
###################################################
par(mfrow=c(1,1))
x=c(1,1,2, 3.3, 3.3, 4,4,4,4.7,5,5.6,5.6,5.6,6,6,6.5,6.9)
y=c(10.84,9.3,16.35, 22.88, 24.35, 24.56, 25.86, 29.16,24.59,22.25,25.9, 27.2,
25.61,25.45,26.56, 21.03, 21.46)
print(rbind(x,y))
plot(x,y, pch=16)
g=lm(y~x)
abline(coef(g), col=2)


###################################################
### chunk number 24: 
###################################################
xtable(summary(g)) #test for the significance of regression?
xtable(anova(g))


###################################################
### chunk number 25: 
###################################################
X=factor(x)  #re-assigne x as a factor X 
levels(X)  #level for the factor X
length(levels(X))
f=lm(y~X)
xtable(anova(f))


###################################################
### chunk number 26: 
###################################################
anova(g,f)