# November 12, 2024

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)

ex1.lm<-lm(y~x1+x2)
summary(ex1.lm)

X<-matrix(c(rep(1,5),x1,x2),ncol=3)
n<-nrow(X)
p<-ncol(X)-1
bhat<-solve(t(X)%*%X)%*%t(X)%*%y
yhat<-X%*%bhat
ehat<-y-yhat
rss<-t(ehat)%*%(ehat)
sigma2hat<-rss/(n-p-1)

aic<-n*log(rss/n)+2*p
bic<-n*log(rss/n)+p*log(n)
cp<-rss/sigma2hat+2*p-n

aic
bic
cp

summary(lm(y~x1))

X1<-matrix(c(rep(1,5),x1),ncol=2)
p1<-ncol(X1)-1
bhat1<-solve(t(X1)%*%X1)%*%t(X1)%*%y
yhat1<-X1%*%bhat1
ehat1<-y-yhat1
rss1<-t(ehat1)%*%(ehat1)

aic1<-n*log(rss1/n)+2*p1
bic1<-n*log(rss1/n)+p1*log(n)
cp1<-rss1/sigma2hat+2*p1-n

aic1
aic

bic1
bic

cp1
cp

ex1.1.lm<-lm(y~x1)
summary(ex1.1.lm)

AIC(ex1.lm)
BIC(ex1.lm)

AIC(ex1.1.lm)
BIC(ex1.1.lm)

AIC(ex1.lm, ex1.1.lm)
BIC(ex1.lm, ex1.1.lm)


### Example 2

swiss
?swiss
?step

dim(swiss)
plot(swiss)
cor(swiss)
summary(swiss)

summary(lm1 <- lm(Fertility ~ ., data = swiss))
slm1 <- step(lm1)
summary(slm1)
slm1$anova

slm1be <- step(lm1, direction="backward")
summary(slm1be)
slm1be$anova

summary(lm0 <- lm(Fertility ~ 1, data = swiss))

slm0fs <- step(lm0, scope=~Agriculture+Examination+Education+Catholic+Infant.Mortality, direction="forward")
summary(slm0fs)
slm0fs$anova

summary(slm1)
summary(slm1be)
summary(slm0fs)

AIC(lm1)
AIC(slm1)
BIC(lm1)
BIC(slm1)

AIC(lm1, slm1)
BIC(lm1, slm1)


### Example 3

library(boot)

?cv.glm

data(mammals, package="MASS")
mammals
?mammals
summary(lm(log(brain) ~ log(body), data = mammals))
summary(glm(log(brain) ~ log(body), data = mammals))
mammals.glm <- glm(log(brain) ~ log(body), data = mammals)

(cv.err <- cv.glm(mammals, mammals.glm)$delta)
(cv.err.6 <- cv.glm(mammals, mammals.glm, K = 6)$delta)

set.seed(123)
(cv.err.6 <- cv.glm(mammals, mammals.glm, K = 6)$delta)

muhat <- fitted(mammals.glm)
mammals.diag <- glm.diag(mammals.glm)
?glm.diag
mean((mammals.glm$y - muhat)^2/(1 - mammals.diag$h)^2)
cv.err


### Example 4

cv.glm(swiss,lm1)
cv.glm(swiss,lm1)$delta

glm1<-glm(Fertility ~ ., data = swiss)
summary(glm1)
summary(lm1)

cv.glm(swiss,glm1)
cv.glm(swiss,glm1)$delta

sglm1 <- step(glm1)
summary(sglm1)

cv.glm(swiss,sglm1)$delta

cv.glm(swiss,glm1)$delta[1]
cv.glm(swiss,sglm1)$delta[1]


### Example 5

ais
?ais

plot(ais)

ais.glm<-glm(Ferr~Ht+Wt+LBM+BMI+SSF+Bfat, data=ais)
summary(ais.glm)

ais.glm.step <- step(ais.glm)
summary(ais.glm.step)
ais.glm.step$anova

ais.glm.be <- step(ais.glm, direction="backward")
summary(ais.glm.be)
ais.glm.be$anova

ais.glm0<-glm(Ferr~1, data=ais)
summary(ais.glm0)
ais.glm.fs <- step(ais.glm0, scope=~Ht+Wt+LBM+BMI+SSF+Bfat, direction="forward")
summary(ais.glm.fs)
ais.glm.fs$anova

plot(ais$BMI, ais$Bfat)
cor(ais$BMI, ais$Bfat)
cor.test(ais$BMI, ais$Bfat)

plot(fitted(ais.glm.step), resid(ais.glm.step))
abline(h=0, lty=2)

qqnorm(resid(ais.glm.step))
qqline(resid(ais.glm.step))

cv.glm(ais,ais.glm)$delta[1]
cv.glm(ais,ais.glm.step)$delta[1]