# Sep 24, 2024

### Example 1
### Simple regression without matrix
### See HW 1 and 2

x<-c(1,2,3,4,5)
y<-c(1.9, 2.2, 2.9, 3.2, 5.2)

xbar<-mean(x)
ybar<-mean(y)

x-xbar
(x-xbar)^2

sxx<-sum((x-xbar)^2)
syy<-sum((y-ybar)^2)
sxy<-sum((x-xbar)*(y-ybar))

beta1hat<-sxy/sxx
beta0hat<-ybar-beta1hat*xbar

n<-length(y)

rss<-syy-sxy^2/sxx
sigma2hat<-rss/(n-2)

se1<-sqrt(sigma2hat/sxx)
se0<-sqrt(sigma2hat*(1/n+mean(x)^2/sxx))

t0<-beta0hat/se0
t1<-beta1hat/se1

pt(abs(t0),n-2,lower.tail=F)*2
pt(abs(t1),n-2,lower.tail=F)*2

r2<-1-rss/syy
adjr2<-1-(rss/(n-2))/(syy/(n-1))


### Example 2
### Simple regression with matrix

X<-matrix(c(rep(1,5),x),ncol=2)
X

bhat<-solve(t(X)%*%X)%*%t(X)%*%y
bhat

beta0hat<-bhat[1,1]
beta1hat<-bhat[2,1]

yhat<-X%*%bhat
ehat<-y-yhat
t(ehat)%*%(ehat)
t(y)%*%y - t(yhat)%*%yhat
rss
rss<-t(ehat)%*%(ehat)

sigma2hat
sigma2hat<-rss/(5-1-1)

# residual se
sqrt(sigma2hat)

sigma2hat*solve(t(X)%*%X)
sigma2hat[1,1]*solve(t(X)%*%X)
varhatbhat<-sigma2hat[1,1]*solve(t(X)%*%X)

se0
se1
sqrt(varhatbhat[1,1])
sqrt(varhatbhat[2,2])
se0<-sqrt(varhatbhat[1,1])
se1<-sqrt(varhatbhat[2,2])
cov01<-varhatbhat[1,2]
cov01
-sigma2hat[1,1]*mean(x)/sxx

# same as before

t0<-beta0hat/se0
t1<-beta1hat/se1
t0
t1

pt(abs(t0),5-1-1,lower.tail=F)*2
pt(abs(t1),5-1-1,lower.tail=F)*2

syy<-sum((y-mean(y))^2)
ssreg<-syy-rss

r2<-1-rss/syy
adjr2<-1-(rss/(5-1-1))/(syy/(5-1))
r2
adjr2
1-(1-r2)*(5-1)/(5-1-1)

# 95% prediction at xstar=3.6

ystar<-beta0hat+beta1hat*3.6
sepred<-sqrt(sigma2hat*(1+1/5+(3.6-mean(x))^2/sxx))
qt(0.025,5-2,lower.tail=F)
ystar-qt(0.025,5-2,lower.tail=F)*sepred
ystar+qt(0.025,5-2,lower.tail=F)*sepred

xstar<-c(1, 3.6)
t(xstar)%*%bhat
ystar
sqrt(sigma2hat*(1+t(xstar)%*%solve(t(X)%*%X)%*%xstar))
sepred



### Example 3

library(alr4)

Forbes

x<-Forbes$bp
y<-Forbes$lpres

xbar<-mean(x)
ybar<-mean(y)

sxx<-sum((x-xbar)^2)
syy<-sum((y-ybar)^2)
sxy<-sum((x-xbar)*(y-ybar))

beta1hat<-sxy/sxx
beta0hat<-ybar-beta1hat*xbar

n<-length(y)

rss<-syy-sxy^2/sxx
sigma2hat<-rss/(n-2)

se1<-sqrt(sigma2hat/sxx)
se0<-sqrt(sigma2hat*(1/n+mean(x)^2/sxx))

t0<-beta0hat/se0
t1<-beta1hat/se1

pt(abs(t0),n-2,lower.tail=F)*2
pt(abs(t1),n-2,lower.tail=F)*2

# 90% CI for beta0

qt(0.05,n-2,lower.tail=F)
beta0hat-qt(0.05,n-2,lower.tail=F)*se0
beta0hat+qt(0.05,n-2,lower.tail=F)*se0

# 95% CI for beta1

qt(0.025,n-2,lower.tail=F)
beta1hat-qt(0.025,n-2,lower.tail=F)*se1
beta1hat+qt(0.025,n-2,lower.tail=F)*se1

# 99% Prediction Interval, at Xstar=200

ystar200<-beta0hat+beta1hat*200
sepred200<-sqrt(sigma2hat*(1+1/n+(200-mean(x))^2/sxx))

qt(0.005,n-2,lower.tail=F)
ystar200-qt(0.005,n-2,lower.tail=F)*sepred200
ystar200+qt(0.005,n-2,lower.tail=F)*sepred200

ll<-ystar200-qt(0.005,n-2,lower.tail=F)*sepred200
lu<-ystar200+qt(0.005,n-2,lower.tail=F)*sepred200

10^(ll/100)
10^(lu/100)

r2<-1-rss/syy
adjr2<-1-(rss/(n-2))/(syy/(n-1))

yhat<-beta0hat+beta1hat*x
y-yhat
summary(y-yhat)

sum((y-yhat)^2)