###################################################
### 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
pairs(EX3)


###################################################
### chunk number 3: 
###################################################
#plot the 3d data 
#install.packages("scatterplot3d")
#library(scatterplot3d)
#scatterplot3d(CAS, DIS, TIM, highlight.3d=TRUE )

#install.packages("rgl")
#library(rgl)
#rgl.plot3d(CAS, DIS, TIM)


###################################################
### chunk number 4: 
###################################################
fm1=lm(TIM~CAS+DIS, data=EX3)
anova(fm1)  
summary(fm1)

fitted(fm1)	#Y.hat
resid(fm1)	#e :residual
df.residual(fm1)  #n-p
coef(fm1)  	  #beta.hat's
SSres=sum(resid(fm1)^2);SSres	#sum of squares due to residuals
MSres=SSres/df.residual(fm1); MSres  #Residual Mean Square 
sqrt(MSres)  


X=model.matrix(fm1) #assign the model matrix of the linear model to X
head(X)
t(X) %*% X	#X'X
t(X) %*% TIM 	#X'Y


solve(t(X) %*% X) %*% t(X) %*% TIM   #inv(X'X)X'Y
solve(t(X) %*% X)     #inv(X'X)

C=solve(t(X) %*% X)   #assign C to be the inverse of the X'X matrix
sqrt(MSres*diag(C))   #Standard error of the estimators of beta's
		      #diag(C) diagonal elemenets of matrix C

vcov(fm1)  #variance and covariance matrix for beta.hat's
sqrt(diag(vcov(fm1))) # standard errors of beta.hat's

confint(fm1)   #confident interval for the estimators of beta's



###################################################
### chunk number 5: 
###################################################
#install.packages("ellipse")  #install the package "ellipse" 
library(ellipse)  #call-in the library
plot(ellipse(fm1, which=c(2,3)), type="l", xlim=c(1,2.2))  
		  #confidence region for CAS and DIS
abline(v=confint(fm1)[2,], lty=2) #confidence interval for CAS
abline(h=confint(fm1)[3,], lty=2) #confidence interval for DIS

cor(CAS, DIS)


###################################################
### chunk number 6: 
###################################################
par(mfrow=c(2,2))
plot(fm1, ask=F)


###################################################
### chunk number 7: 
###################################################
fm2=lm(TIM~CAS)
anova(fm2)
anova(fm2, fm1)