dat1=read.table("B3.csv", header=T, sep=",")
head(dat1) #see the first six rows of dat1
y=dat1[,1] #first column of dat1
x1=dat1[,2] #second colum of dat1
x6=dat1[,7]
##########
#attach the data.frame to use the header directly
#########
# attach(dat1)
# y
# x1

f=lm(y~x1)
summary(f)
coef(f)  #coefficients for the linear regression model
anova(f)
confint(f)  #confidence intervals for beta0, beta1
##CI for mean response and PI for a new obs ##
predict(f, data.frame(x1=275), interval="confidence")
predict(f, data.frame(x1=275), interval="prediction")

################
##done by hand##
################
new=matrix(c(1, 275));new
X=model.matrix(f)

s2=sum(resid(f)^2)/df.residual(f);s2  ##MSres=SSres/df.res
s=sqrt(s2);s
t.lim=qt(0.975,30); t.lim
y.new= t(new) %*% coef(f); y.new
hii.new=t(new) %*% solve( t(X) %*% X ) %*% new; hii.new     #x0'(X'X)^{-1}x0

se.Ey=s*sqrt(hii.new)
c(y.new - t.lim * se.Ey, y.new + t.lim * se.Ey)

se.y=s*sqrt(1+hii.new)
c(y.new - t.lim * se.y, y.new + t.lim * se.y)



f3=lm(y~x1+x6)
summary(f3)
anova(f3)
coef(f3)
confint(f3)
predict(f3, data.frame(x1=275,x6=2), interval="confidence")
predict(f3, data.frame(x1=275,x6=2), interval="prediction")


plot(x1, y)  #x-axis vs y-axis
#identify(x1, y) #Esc for stop identify
X=model.matrix(f)
H=X %*% solve(t(X) %*% X) %*% t(X)
diag(H) #diagnoal elements of H
sum(diag(H))
x22=matrix(c(1, 360), nrow=2) #2*1 vector
x17=matrix(c(1, 500), nrow=2) #2*1 vector
x16=matrix(c(1, 302), nrow=2) #2*1 vector
x12=matrix(c(1, 85), nrow=2) #2*1 vector
xi=x22
hii=t(xi) %*% solve(t(X)%*%X) %*% xi ; hii
xi=x17
hii=t(xi) %*% solve(t(X)%*%X) %*% xi ; hii
xi=x12; xj=x16;
hij=t(xi) %*% solve(t(X)%*%X) %*% xj ; hij



dat1=read.table("ex-3-1.csv", header=T, sep=",")
head(dat1) #see the first six rows of dat1
y=dat1[,2] 
can=dat1[,3]
dist=dat1[,4]

f=lm(y~can+dist)

yhat=fitted(f)
ei=resid(f);ei
di=rstandard(f);di
ti=rstudent(f);ti
hii=hatvalues(f); hii
out=cbind("y"=y, "yhat"=yhat, "ei"=ei, "rstudent"=di, "deletedr"=ti, "hii"=hii)
out
par(mfrow=c(2,2)) #
plot( yhat, ei)
plot( yhat, di)
plot( yhat, ti)
qqnorm(di);qqline(di) #normal probability plot