Chapter 7 Hmwk. STA-6543-91T-202630
Mark Hertzfeld (OPT121)
2026-07-03
6. In this exercise, you will further analyze the Wage data set considered throughout this chapter.
library("ISLR2")
library("leaps")
library("gam")
library("boot")
data("Wage")
data("College")(a) Perform polynomial regression to predict wage using age. Use cross-validation to select the optimal degree d for the polynomial. What degree was chosen, and how does this compare to the results of hypothesis testing using ANOVA? Make a plot of the resulting polynomial fit to the data.
Polynomial Regression:
set.seed(70)
cv.errorWage = rep(NA, 12)
for(i in 1:12){
Wage.fitplot = glm(wage ~ poly(age, i), data = Wage)
cv.errorWage[i] = cv.glm(Wage, Wage.fitplot, K = 10)$delta[1]
}
cv.errorWage## [1] 1676.377 1600.109 1596.298 1595.526 1595.232 1594.232 1596.368 1595.009
## [9] 1593.771 1596.706 1598.234 1596.669which.min(cv.errorWage)## [1] 9The lowest degree of freedom for the polynomial regression from part A is to have a DF of 9.
However, it looks like the most significant DF values are 2 or 3 because of the decrease that occurs in coefficient values.
I would say that the optimal value for the degrees of freedom is actually 3, just because that is the last value that appears to have a decrease before the values begin plateauing.
Poly Regress D = 9, Personally think is best D = 3
ANOVA:
Wagefit.1 = lm(wage~age, data = Wage)
Wagefit.2 = lm(wage~poly(age,2), data = Wage)
Wagefit.3 = lm(wage~poly(age,3), data = Wage)
Wagefit.4 = lm(wage~poly(age,4), data = Wage)
Wagefit.5 = lm(wage~poly(age,5), data = Wage)
# Nested Models
anova(Wagefit.1, Wagefit.2, Wagefit.3, Wagefit.4, Wagefit.5)## 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 ' ' 1From the ANOVA test, it appears that the last DF value considered statistically significant is 3.
The values for the optimized value of the Degrees of Freedom appears to be 3 for both.
Wage.fitplot2 = lm(wage~poly(age,3), data = Wage)
agelims=range(Wage$age)
#agelims
age.grid = seq(from=agelims[1], to = agelims[2])
preds.plot=predict(Wage.fitplot2,newdata=list(age=age.grid),se=TRUE)
se.bands=cbind(preds.plot$fit+2*preds.plot$se.fit,
preds.plot$fit-2*preds.plot$se.fit)
plot(Wage$age,Wage$wage,xlim=agelims, cex=.5, col ="darkgrey")
title('Degree-3 Polynomial')
lines(age.grid, preds.plot$fit, lwd = 2, col ='darkblue')
matlines(age.grid, se.bands, lwd = 1, col = "lightblue", lty = 3) 
(b) Fit a step function to predict wage using age, and perform cross validation to choose the optimal number of cuts. Make a plot of the fit obtained.
cv.errorWageB <- rep(NA, 10)
for(i in 2:10){
breaks <- seq(min(Wage$age), max(Wage$age),length.out = i + 1)
Wage$age.cut= cut(Wage$age, breaks = breaks,include.lowest = TRUE)
fitWageB = glm(wage~age.cut, data = Wage)
cv.errorWageB[i-1] = cv.glm(Wage, fitWageB)$delta[1]
}
cv.errorWageB## [1] 1734.064 1682.763 1635.894 1631.450 1623.291 1611.604 1601.006 1610.213
## [9] 1604.804 NAwhich.min(cv.errorWageB) + 1## [1] 8step.fitWage <- glm(wage ~ cut(age, 8), data = Wage)
#summary(step.fitWage)
age.grid2 = seq(from=agelims[1], to = agelims[2])
preds.plot2=predict(step.fitWage,newdata=list(age=age.grid),se=TRUE)
se.bands2=cbind(preds.plot2$fit+2*preds.plot2$se.fit,
preds.plot2$fit-2*preds.plot2$se.fit)
plot(Wage$age,Wage$wage,xlim=agelims, cex=.5, col ="darkgrey")
title('Model with 8 Cuts')
lines(age.grid2, preds.plot2$fit, lwd = 2, col ='darkblue')
matlines(age.grid2, se.bands2, lwd = 1, col = "lightblue", lty = 3) 
After using cross validation on the model, it appears that the optimal value for cuts is 8.
I then fitted the above model with 8 cuts.
10. This question relates to the College data set. (a) Split the data into a training set and a test set. Using out-of-state tuition as the response and the other variables as the predictors, perform forward stepwise selection on the training set in order to identify a satisfactory model that uses just a subset of the predictors.
set.seed(70)
College=na.omit(College)
trainCol=sample(1:nrow(College), nrow(College)/2)
testCol=setdiff(1:nrow(College), trainCol)
test.ColMatx=model.matrix(Outstate~.,data = College, subset = testCol)fit.CollegeFwd = regsubsets(Outstate~.,data = College,
subset = trainCol, nvmax=17,method="forward")
summary(fit.CollegeFwd)## Subset selection object
## Call: regsubsets.formula(Outstate ~ ., data = College, subset = trainCol,
## nvmax = 17, method = "forward")
## 17 Variables (and intercept)
## Forced in Forced out
## PrivateYes FALSE FALSE
## Apps FALSE FALSE
## Accept FALSE FALSE
## Enroll FALSE FALSE
## Top10perc FALSE FALSE
## Top25perc FALSE FALSE
## F.Undergrad FALSE FALSE
## P.Undergrad FALSE FALSE
## Room.Board FALSE FALSE
## Books FALSE FALSE
## Personal FALSE FALSE
## PhD FALSE FALSE
## Terminal FALSE FALSE
## S.F.Ratio FALSE FALSE
## perc.alumni FALSE FALSE
## Expend FALSE FALSE
## Grad.Rate FALSE FALSE
## 1 subsets of each size up to 17
## Selection Algorithm: forward
## PrivateYes Apps Accept Enroll Top10perc Top25perc F.Undergrad
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " " " " "
## 4 ( 1 ) "*" " " " " " " " " " " " "
## 5 ( 1 ) "*" " " " " " " " " " " " "
## 6 ( 1 ) "*" " " " " " " " " " " " "
## 7 ( 1 ) "*" " " " " " " " " " " " "
## 8 ( 1 ) "*" " " " " " " " " " " " "
## 9 ( 1 ) "*" " " "*" " " " " " " " "
## 10 ( 1 ) "*" " " "*" "*" " " " " " "
## 11 ( 1 ) "*" "*" "*" "*" " " " " " "
## 12 ( 1 ) "*" "*" "*" "*" " " "*" " "
## 13 ( 1 ) "*" "*" "*" "*" " " "*" "*"
## 14 ( 1 ) "*" "*" "*" "*" " " "*" "*"
## 15 ( 1 ) "*" "*" "*" "*" " " "*" "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## P.Undergrad Room.Board Books Personal PhD Terminal S.F.Ratio
## 1 ( 1 ) " " "*" " " " " " " " " " "
## 2 ( 1 ) " " "*" " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " " " " "
## 4 ( 1 ) " " "*" " " " " " " " " " "
## 5 ( 1 ) " " "*" " " " " " " " " " "
## 6 ( 1 ) " " "*" " " " " "*" " " " "
## 7 ( 1 ) " " "*" " " "*" "*" " " " "
## 8 ( 1 ) " " "*" " " "*" "*" "*" " "
## 9 ( 1 ) " " "*" " " "*" "*" "*" " "
## 10 ( 1 ) " " "*" " " "*" "*" "*" " "
## 11 ( 1 ) " " "*" " " "*" "*" "*" " "
## 12 ( 1 ) " " "*" " " "*" "*" "*" " "
## 13 ( 1 ) " " "*" " " "*" "*" "*" " "
## 14 ( 1 ) "*" "*" " " "*" "*" "*" " "
## 15 ( 1 ) "*" "*" " " "*" "*" "*" "*"
## 16 ( 1 ) "*" "*" " " "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## perc.alumni Expend Grad.Rate
## 1 ( 1 ) " " " " " "
## 2 ( 1 ) "*" " " " "
## 3 ( 1 ) "*" "*" " "
## 4 ( 1 ) "*" "*" " "
## 5 ( 1 ) "*" "*" "*"
## 6 ( 1 ) "*" "*" "*"
## 7 ( 1 ) "*" "*" "*"
## 8 ( 1 ) "*" "*" "*"
## 9 ( 1 ) "*" "*" "*"
## 10 ( 1 ) "*" "*" "*"
## 11 ( 1 ) "*" "*" "*"
## 12 ( 1 ) "*" "*" "*"
## 13 ( 1 ) "*" "*" "*"
## 14 ( 1 ) "*" "*" "*"
## 15 ( 1 ) "*" "*" "*"
## 16 ( 1 ) "*" "*" "*"
## 17 ( 1 ) "*" "*" "*"val.errorsCollege=rep(NA,17) # 17 variable models
for(i in 1:17){
coefi=coef(fit.CollegeFwd,id = i) # Take coefficients from best fit model
predCollege=test.ColMatx[,names(coefi)]%*%coefi #all rows and names of coefficients we just got
val.errorsCollege[i]=mean((College$Outstate[testCol]-predCollege)^2)
}
val.errorsCollege## [1] 22083056 24851771 26199505 26929794 27689832 28140450 28087969 28124835
## [9] 28318099 29304111 28984067 29089193 29025317 29019568 29054024 29064324
## [17] 29062875which.min(val.errorsCollege)## [1] 1coef(fit.CollegeFwd, id = 4)## (Intercept) PrivateYes Room.Board perc.alumni Expend
## -747.5697120 2283.2696057 1.2684884 81.4175640 0.2250401The lowest value returned is 1, but a more flexible model with more predictors appear to be with 4 predictors.
So, I chose to go with 4 over 1.
(b) Fit a GAM on the training data, using out-of-state tuition as the response and the features selected in the previous step as the predictors. Plot the results, and explain your findings.
gam.fitCol = gam(Outstate~s(Room.Board), data = College, subset= trainCol)
plot(gam.fitCol, se=T, col = "blue")
gam.fitCol1A = gam(Outstate~s(Room.Board)+s(perc.alumni)+s(Expend)+Private,
data = College, subset= trainCol)
par(mfrow=c(1,4))
plot(gam.fitCol1A, se=T, col = "blue")
The results appear to show:
As Room and Board prices increases, so to does out of state tuition, which makes sense when you think about it.
As the percentage of alumni who donate to the school increases, so to does the cost for out of state tuition. This is an interesting one because you would think that out of state tuition would decrease due to alumni giving, but then there’s the age old question “how loyal/much are you willing to pay to get something you want”. Got to love greed.
As the educational expenditures per student increase, so to does out of state tuition. This makes sense as if there’s an increase price to house students, this will be reflected in tuition across the board.
If the school is private, out of tuition is extremely high. This makes sense due to the lack of funding from the state that public schools are able to experience.
In general, out of state tuition is more expensive then in state tuition.
(c) Evaluate the model obtained on the test set, and explain the results obtained.
gam.fitCol1B = gam(Outstate~s(Room.Board)+s(perc.alumni)+s(Expend)+Private,
data = College, subset= testCol)
par(mfrow=c(1,4))
plot(gam.fitCol1B, se=T, col = "red")
The results appear to be pretty similar to what was seen when fitted on the training model.
The only clear differences are:
The results for the plot of Room and Board vs. Out of State tuition are lower/have less of a shepherd’s crook
The results for Expenditure vs Out of State tuition are missing the hump that was seen in the training model.
(d) For which variables, if any, is there evidence of a non-linear relationship with the response? There appears to be evidence that the ‘perc.alumni’, ‘Room.Board’, and ‘Expend’ variables have a non-linear response with ‘Outstate’.
Perc.alumni appears to be the most non-linear out of the 3 overall.
Expend appears to start linear, but suddenly deviate the larger it gets.