# October 22, 2024

### Example 1
### Page 157

library(alr4)

physics
?physics

summary(lm(y~x, data=physics))

physics.lm<-lm(y~x, weights=1/SD^2, data=physics)

summary(physics.lm)
anova(physics.lm)

y<-physics$y
x<-physics$x
w<-physics$SD

n<-length(y)
p<-1

X<-matrix(c(rep(1,n),x),ncol=2)
W<-diag(1/w^2)
W
bhat<-solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%y
bhat

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

sqrt(sigma2hat)

# Watch out for W
rss<-t(ehat)%*%W%*%(ehat)
sigma2hat<-rss/(n-p-1)

sqrt(sigma2hat)

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

se<-sqrt(diag(varhatbhat))
se

tstat<-bhat/se
tstat

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

syy<-sum((y-mean(y))^2)
ssreg<-syy-rss
fstat<-(ssreg/p)/sigma2hat
fstat

# weighted mean
wmean.y<-sum(W%*%y)/sum(diag(W))
syy<-t(y-wmean.y)%*%W%*%(y-wmean.y)  
# or sum(W%*%(y-wmean.y)^2)
ssreg<-syy-rss
fstat<-(ssreg/p)/sigma2hat
fstat

r2<-1-rss/syy
r2
adjr2<-1-(1-r2)*(n-1)/(n-p-1)
adjr2


### Example 2

Mitchell
?Mitchell

plot(Temp~Month, data=Mitchell)
plot(Temp~Month, data=Mitchell, type="l")

mitchell.lm<-lm(Temp~Month, data=Mitchell)
summary(mitchell.lm)

mitchell.wlm<-lm(Temp~Month, weight=1/(Month+1), data=Mitchell)
summary(mitchell.wlm)

plot(Temp~Month, data=Mitchell, type="l")
abline(mitchell.lm, col="blue")
abline(mitchell.wlm, col="red")

fit.lm<-fitted(mitchell.lm)
resid.lm<-resid(mitchell.lm)

plot(fit.lm, resid.lm)
abline(h=0, lty=2)

fit.wlm<-fitted(mitchell.wlm)
resid.wlm<-resid(mitchell.wlm)

plot(fit.wlm, resid.wlm)
abline(h=0, lty=2)

par(mfrow=c(1,2))
plot(fit.lm, resid.lm)
abline(h=0, lty=2)
plot(fit.wlm, resid.wlm)
abline(h=0, lty=2)

# reset plot
par(mfrow=c(1,1))


### Example 3

brains
?brains

summary(brains)
plot(brains)

plot(BrainWt~BodyWt, data=brains)
text(brains[brains$BrainWt>1000,2],brains[brains$BrainWt>1000,1],row.names(brains[brains$BrainWt>1000,]))

plot(log(BrainWt)~log(BodyWt), data=brains)

brainslog.lm<-lm(log(BrainWt)~log(BodyWt), data=brains)
summary(brainslog.lm)
abline(brainslog.lm)

fit.brainslog<-fitted(brainslog.lm)
resid.brainslog<-resid(brainslog.lm)
plot(fit.brainslog,resid.brainslog)
abline(h=0, lty=2)

brains.lm<-lm(BrainWt~BodyWt, data=brains)
summary(brains.lm)

fit.brains<-fitted(brains.lm)
resid.brains<-resid(brains.lm)
plot(fit.brains,resid.brains)
abline(h=0, lty=2)


### Example 4

airquality
?airquality

summary(airquality)

plot(airquality[,1:4])

fit.aq <- lm(Ozone ~ Solar.R + Wind + Temp, data=airquality)
summary(fit.aq)

plot(fitted(fit.aq),resid(fit.aq))
abline(h=0,lty=2)

fit.aq.log<-lm(log(Ozone) ~ Solar.R + Wind + Temp, data=airquality)
summary(fit.aq.log)

plot(fitted(fit.aq.log),resid(fit.aq.log))
abline(h=0,lty=2)

fit.aq.sqrt<-lm(sqrt(Ozone) ~ Solar.R + Wind + Temp, data=airquality)
summary(fit.aq.sqrt)

plot(fitted(fit.aq.sqrt),resid(fit.aq.sqrt))
abline(h=0,lty=2)

fit.aq.inv<-lm(I(1/Ozone) ~ Solar.R + Wind + Temp, data=airquality)
summary(fit.aq.inv)

plot(fitted(fit.aq.inv),resid(fit.aq.inv))
abline(h=0,lty=2)

b<-1/3
fit.aq.p<-lm(I(Ozone^b) ~ Solar.R + Wind + Temp, data=airquality)
summary(fit.aq.p)

plot(fitted(fit.aq.p),resid(fit.aq.p))
abline(h=0,lty=2)

par(mfrow=c(2,2))
plot(fitted(fit.aq.log),resid(fit.aq.log))
abline(h=0,lty=2)
plot(fitted(fit.aq.sqrt),resid(fit.aq.sqrt))
abline(h=0,lty=2)
plot(fitted(fit.aq.inv),resid(fit.aq.inv))
abline(h=0,lty=2)
plot(fitted(fit.aq.p),resid(fit.aq.p))
abline(h=0,lty=2)

par(mfrow=c(1,1))

### Example 5
### Problem 8.5, page 202

BigMac2003
?BigMac2003

summary(BigMac2003)
plot(BigMac2003)

fit.bm.f<-lm(BigMac~FoodIndex, data=BigMac2003)
summary(fit.bm.f)

plot(fitted(fit.bm.f),resid(fit.bm.f))
abline(h=0,lty=2)

library(MASS)
?boxcox

boxcox(fit.bm.f)
# boxcox(BigMac~FoodIndex, data=BigMac2003)

bc.bm.f<-boxcox(fit.bm.f)

bc.bm.f$x[which.max(bc.bm.f$y)]
bc.bm.f.lambda<-bc.bm.f$x[which.max(bc.bm.f$y)]

# "Box-Cox"
# Note the Box-Cox transform usually refers to 
# the "scale power" on page 189

fit.bm.f.bc<-lm(((BigMac^bc.bm.f.lambda-1)/bc.bm.f.lambda)~FoodIndex, data=BigMac2003)
summary(fit.bm.f.bc)

plot(BigMac2003$FoodIndex, ((BigMac2003$BigMac^bc.bm.f.lambda-1)/bc.bm.f.lambda))
plot(BigMac2003$FoodIndex, BigMac2003$BigMac)

plot(fitted(fit.bm.f.bc),resid(fit.bm.f.bc))
abline(h=0,lty=2)

par(mfrow=c(1,2))
plot(fitted(fit.bm.f),resid(fit.bm.f))
abline(h=0,lty=2)
plot(fitted(fit.bm.f.bc),resid(fit.bm.f.bc))
abline(h=0,lty=2)

# You may continue with further parts if you wish...

par(mfrow=c(1,1))

#