# Chapter 7 page 297: 6,10
# 6
library(ISLR)
## Warning: package 'ISLR' was built under R version 4.1.3
fit1=lm(wage~age,data=Wage)
fit2=lm(wage~poly(age,2),data=Wage)
fit3=lm(wage~poly(age,3),data=Wage)
fit4=lm(wage~poly(age,4),data=Wage)
fit5=lm(wage~poly(age,5),data=Wage)
anova(fit1,fit2,fit3,fit4,fit5)
## Analysis of Variance Table
##
## Model 1: wage ~ age
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.5931 < 2.2e-16 ***
## 3 2996 4777674 1 15756 9.8888 0.001679 **
## 4 2995 4771604 1 6070 3.8098 0.051046 .
## 5 2994 4770322 1 1283 0.8050 0.369682
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(boot)
set.seed(1)
deltas=rep(NA,10)
for(i in 1:10)
{
fit=glm(wage~poly(age,i),data=Wage)
deltas[i]=cv.glm(Wage,fit,K=10)$delta[1]
}
plot(1:10,deltas,xlab="Degree",ylab="Test MSE",type="l")
d.min=which.min(deltas)
points(which.min(deltas),deltas[which.min(deltas)],col="pink",cex=2,pch=20)
cvs=rep(NA,10)
for(i in 2:10)
{
Wage$age.cut=cut(Wage$age,i)
fit=glm(wage~age.cut,data=Wage)
cvs[i]=cv.glm(Wage,fit,K=10)$delta[1]
}
plot(2:10,cvs[-1],xlab="Cuts",ylab="Test MSE",type="l")
d.min=which.min(cvs)
points(which.min(cvs),cvs[which.min(cvs)],col="pink",cex=2,pch=20)
plot(wage~age,data=Wage,col="red")
agelims<-range(Wage$age)
age.grid<-seq(from=agelims[1],to=agelims[2])
fit<-glm(wage~cut(age,8),data=Wage)
preds<-predict(fit,data.frame(age=age.grid))
lines(age.grid,preds,col="pink",lwd=2)
# 10
library(leaps)
## Warning: package 'leaps' was built under R version 4.1.3
set.seed(1)
attach(College)
train=sample(length(Outstate),length(Outstate)/2)
test<--train
College.train=College[train,]
College.test=College[test,]
fit=regsubsets(Outstate~.,data=College.train,nvmax=17,method="forward")
fit.summary=summary(fit)
par(mfrow=c(1,3))
plot(fit.summary$cp,xlaba="Number of variables",ylab="Cp",type="l")
## Warning in plot.window(...): "xlaba" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "xlaba" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "xlaba" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "xlaba" is not a
## graphical parameter
## Warning in box(...): "xlaba" is not a graphical parameter
## Warning in title(...): "xlaba" is not a graphical parameter
min.cp=min(fit.summary$cp)
std.cp=sd(fit.summary$cp)
abline(h=min.cp+0.2*std.cp,col="pink",lty=2)
abline(h=min.cp-0.2*std.cp,col="pink",lty=2)
plot(fit.summary$bic,xlab="Number of variables",ylab="BIC",type="l")
min.bic=min(fit.summary$bic)
std.bic=sd(fit.summary$bic)
abline(h=min.bic+0.2*std.bic,col="pink",lty=2)
abline(h=min.bic-0.2*std.bic,col="pink",lty=2)
plot(fit.summary$adjr2,xlab="Number of variables",ylab="Adjusted R2",type="l",ylim=c(0.4,0.84))
max.adjr2=max(fit.summary$adjr2)
std.adjr2=max(fit.summary$adjr2)
abline(h=max.adjr2+0.2*std.adjr2,col="pink",lty=2)
abline(h=max.adjr2-0.2*std.adjr2,col="pink",lty=2)
fit=regsubsets(Outstate~.,data=College,method="forward")
coeffs=coef(fit,id=6)
names(coeffs)
## [1] "(Intercept)" "PrivateYes" "Room.Board" "PhD" "perc.alumni"
## [6] "Expend" "Grad.Rate"
library(gam)
## Warning: package 'gam' was built under R version 4.1.3
## Loading required package: splines
## Loading required package: foreach
## Warning: package 'foreach' was built under R version 4.1.3
## Loaded gam 1.20.1
fit=gam(Outstate~Private+s(Room.Board,df=2)+s(PhD,df=2)+s(perc.alumni,df=2)+s(Expend,df=5)+s(Grad.Rate,df=2),data=College.train)
par(mfrow=c(2,3))
plot(fit,se=T,col="green")
preds=predict(fit,College.test)
err=mean((College.test$Outstate-preds)^2)
err
## [1] 3349290
summary(fit)
##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 2) + s(PhD,
## df = 2) + s(perc.alumni, df = 2) + s(Expend, df = 5) + s(Grad.Rate,
## df = 2), data = College.train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7402.89 -1114.45 -12.67 1282.69 7470.60
##
## (Dispersion Parameter for gaussian family taken to be 3711182)
##
## Null Deviance: 6989966760 on 387 degrees of freedom
## Residual Deviance: 1384271126 on 373 degrees of freedom
## AIC: 6987.021
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 1778718277 1778718277 479.286 < 2.2e-16 ***
## s(Room.Board, df = 2) 1 1577115244 1577115244 424.963 < 2.2e-16 ***
## s(PhD, df = 2) 1 322431195 322431195 86.881 < 2.2e-16 ***
## s(perc.alumni, df = 2) 1 336869281 336869281 90.771 < 2.2e-16 ***
## s(Expend, df = 5) 1 530538753 530538753 142.957 < 2.2e-16 ***
## s(Grad.Rate, df = 2) 1 86504998 86504998 23.309 2.016e-06 ***
## Residuals 373 1384271126 3711182
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## Private
## s(Room.Board, df = 2) 1 1.9157 0.1672
## s(PhD, df = 2) 1 0.9699 0.3253
## s(perc.alumni, df = 2) 1 0.1859 0.6666
## s(Expend, df = 5) 4 20.5075 2.665e-15 ***
## s(Grad.Rate, df = 2) 1 0.5702 0.4506
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1