############## ex6.5 codes

a<-read.table("C:/Documents and Settings/Mandy Jin/My Documents/TA 4315/605.txt")
y<-a$V1
x1<-a$V2
x2<-a$V3
cor.matrix<-cor(a)

f1<-lm(y~x1+x2)     ############## method 1
summary(f1)
resi1<-f1$resi
n<-length(y)
intercept<-rep(1,n)
X<-cbind(intercept,x1,x2)   ######### the design matrix X
b<-solve(t(X)%*%X)%*%t(X)%*%y
resi2<-y-X%*%b
H<-X%*%solve(t(X)%*%X)%*%t(X)
y.hat<-H%*%y
I.n<-diag(n)
resi3<-(I.n-H)%*%y
SSE<-t(y)%*%(I.n-H)%*%y
SSE<-as.numeric(SSE)
MSE<-SSE/13
boxplot(resi3)
par(mfrow=c(2,2))
plot(resi3~y.hat,xlab="fitted value",ylab="residuals",main="plot of residuals against y.hat")
plot(resi3~x1,xlab="x1",ylab="residuals",main="plot of residuals against x1")
plot(resi3~x2,xlab="x2",ylab="residuals",main="plot of residuals against x2")

x1x2<-x1*x2
plot(resi3~x1x2,xlab="x1x2",ylab="residuals",main="plot of residuals against x1x2")
rank.re<-rank(resi3)
normal.resi<-sqrt(MSE)*(qnorm((rank.re-0.375)/(n+0.25)))
plot(resi3~normal.resi,xlab="normal",ylab="residual",main="plot")

e2<-resi1^2
f2<-lm(e2~x1+x2)
a0<-f2$coef[1]
a1<-f2$coef[2]
a2<-f2$coef[3]
SSR.star<-sum((a0+a1*x1+a2*x2-mean(e2))^2)  ####  SSR.star
BP<-SSR.star/2/(SSE/16)^2
p.value<-1-pchisq(BP,2)

qchisq(0.99,2)


aa<-split(a,x2)
bb1<-split(aa[[1]],aa[[1]]$V2)
bb2<-split(aa[[2]],aa[[2]]$V2)

y.mean<-rep(0,8)

for (i in 1:4)
{
y.mean[i]<-mean(bb1[[i]]$V1)
}
for (i in 5:8)
{
y.mean[i]<-mean(bb2[[i-4]]$V1)
}
SSPE<-0
for (i in 1:4)
{

SSPE<- SSPE+sum((bb1[[i]][,1]-y.mean[i])^2)
}
for (i in 5:8)
{

SSPE<-SSPE+sum((bb2[[i-4]][,1]-y.mean[i])^2)

}

SSLF<-SSE-SSPE
 
F<-SSLF/(8-3)/(SSPE/(n-8))