# November 5, 2024

setwd("~/Class/MATH5910/Fall2024/R")

library(alr4)

### Example 1

y<-c(1.9, 2.2, 2.9, 3.2, 5.2)
x1<-c(1:5)
x2<-c(5.1, 2.3, 4.2, 3.3, 4.8)

X<-matrix(c(rep(1,5),x1,x2),ncol=3)
X

H<-X%*%solve(t(X)%*%X)%*%t(X)
H

t(H)
H^2
H%*%H

diag(H)
sum(diag(H))
apply(H,1,sum)
apply(H,2,sum)

Id<-diag(5)
Id

Id-H
t(Id-H)
(Id-H)%*%(Id-H)

diag(Id-H)
sum(diag(Id-H))
apply(Id-H,1,sum)
apply(Id-H,2,sum)

H%*%(Id-H)
round(H%*%(Id-H),2)

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



### Example 2

fuel2001

Tax<-fuel2001$Tax
Income<-fuel2001$Income/1000
Fuel<-1000*fuel2001$FuelC/fuel2001$Pop
Dlic<-1000*fuel2001$Drivers/fuel2001$Pop
LogMiles<-log(fuel2001$Miles)

fuel.lm<-lm(Fuel~Tax+Dlic+Income+LogMiles)
summary(fuel.lm)

par(mar=c(4,4,2,2))
par(mfrow=c(3,2))
plot(Tax, resid(fuel.lm))
abline(h=0,lty=2)
plot(Dlic, resid(fuel.lm))
abline(h=0,lty=2)
plot(Income, resid(fuel.lm))
abline(h=0,lty=2)
plot(LogMiles, resid(fuel.lm))
abline(h=0,lty=2)
plot(fitted(fuel.lm), resid(fuel.lm))
abline(h=0,lty=2)



### Example 3
### see page 217

Forbes

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

par(mfrow=c(1,1))
plot(x,y)

n<-nrow(Forbes)
p<-1

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

H<-X%*%solve(t(X)%*%X)%*%t(X)

hii<-diag(H)
hii
sum(diag(H))
sum(diag(diag(n)-H))

bhat<-solve(t(X)%*%X)%*%t(X)%*%y
yhat<-X%*%bhat
ehat<-y-yhat
rss<-t(ehat)%*%(ehat)
sigma2hat<-rss/(n-p-1)
sigmahat<-sqrt(sigma2hat)
sigmahat

ri<-ehat/(sigmahat[1,1]*sqrt(1-diag(H)))
ri

ti<-ri*sqrt((n-p-2)/(n-p-1-ri^2))
ti

pt(abs(ti),n-p-2,lower.tail=F)*2
n*pt(abs(ti),n-p-2,lower.tail=F)*2
round(n*pt(abs(ti),n-p-2,lower.tail=F)*2, 2)
pmin(1,round(n*pt(abs(ti),n-p-2,lower.tail=F)*2, 2))

Di<-ri^2/(p+1)*(hii/(1-hii))
Di
round(Di,2)

forbes.lm<-lm(lpres~bp, data=Forbes)
summary(forbes.lm)

library(MASS)
studres(forbes.lm)
ti

studres(forbes.lm)-ti

influence.measures(forbes.lm)
Di
hii

### Example 4

snowgeese<-read.table("snowgeese.txt", header=T)
snowgeese

m1 <- lm(photo ~ obs1, snowgeese)
summary(m1)
par(mfrow=c(1,1))
plot(fitted(m1),resid(m1))
abline(h=0,lty=2)

plot(photo ~ obs1, snowgeese)

influence.measures(m1) 

plot(photo ~ obs1, snowgeese)
points(snowgeese$obs1[c(28,29,41)],snowgeese$photo[c(28,29,41)], pch=4)

inf.m1<-influence.measures(m1) 

plot(inf.m1[,5])

inf.m1$infmat
inf.m1$infmat[,5]
inf.m1$infmat[,6]

studres(m1)

# see page 223

par(mfrow=c(3,1))
plot(studres(m1), type="h", xlab=" ", ylab=" ", main="Studentized Residuals")
points(studres(m1), cex=0.7)
abline(h=0, lty=2)
plot(inf.m1$infmat[,6], type="h", xlab=" ", ylab=" ", main="Hat")
points(inf.m1$infmat[,6], cex=0.7)
plot(inf.m1$infmat[,5], type="h", xlab=" ", ylab=" ", main="Cook's Distance")
points(inf.m1$infmat[,5], cex=0.7)

m2<-lm(photo ~ obs1, snowgeese[-29,])
summary(m2)
par(mfrow=c(1,1))
plot(fitted(m2),resid(m2))

m2log<-lm(log(photo) ~ obs1, snowgeese[-29,])
summary(m2log)
plot(fitted(m2log),resid(m2log))

m2sqrt<-lm(sqrt(photo) ~ obs1, snowgeese[-29,])
summary(m2sqrt)
plot(fitted(m2sqrt),resid(m2sqrt))

m3<-lm(sqrt(photo) ~ sqrt(obs1), snowgeese)
summary(m3)
plot(fitted(m3),resid(m3))

m3sub<-lm(sqrt(photo) ~ sqrt(obs1), snowgeese[-29,])
summary(m3sub)
plot(fitted(m3sub),resid(m3sub))
abline(h=0, lty=2)



### Example 5

rat

rat.lm<-lm(y~., data=rat)
summary(rat.lm)
summary(lm(y~BodyWt+Dose, data=rat))
plot(rat)

studres(rat.lm)
rat.inf<-influence.measures(rat.lm) 
rat.inf
rat.inf$infmat
rat.inf$infmat[,7]
rat.inf$infmat[,8]

par(mfrow=c(3,1))
plot(studres(rat.lm), type="h", xlab=" ", ylab=" ", main="Studentized Residuals")
points(studres(rat.lm), cex=0.7)
abline(h=0, lty=2)
plot(rat.inf$infmat[,8], type="h", xlab=" ", ylab=" ", main="Hat")
points(rat.inf$infmat[,8], cex=0.7)
plot(rat.inf$infmat[,7], type="h", xlab=" ", ylab=" ", main="Cook's Distance")
points(rat.inf$infmat[,7], cex=0.7)

rat.lm3<-lm(y~., data=rat[-3,])
summary(rat.lm3)


### Example 6

par(mfrow=c(1,1))

qqnorm(resid(forbes.lm))
qqline(resid(forbes.lm))

qqnorm(resid(m2sqrt))
qqline(resid(m2sqrt))

qqnorm(resid(m3sub))
qqline(resid(m3sub))

qqnorm(resid(rat.lm3))
qqline(resid(rat.lm3))