HW6
NastaranTaabodi -dew050
2026-07-20
#Question 6:
library(ISLR2)
data(Wage)
set.seed(1)
library(boot)
cv.errors <- rep(NA, 10)
for (d in 1:10) {
glm.fit <- glm(wage ~ poly(age, d), data = Wage)
cv.errors[d] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}
cv.errors## [1] 1676.826 1600.763 1598.399 1595.651 1594.977 1596.061 1594.298 1598.134
## [9] 1593.913 1595.950which.min(cv.errors)## [1] 9fit.1 <- lm(wage ~ poly(age, 1), data = Wage)
fit.2 <- lm(wage ~ poly(age, 2), data = Wage)
fit.3 <- lm(wage ~ poly(age, 3), data = Wage)
fit.4 <- lm(wage ~ poly(age, 4), data = Wage)
fit.5 <- lm(wage ~ poly(age, 5), data = Wage)
anova(fit.1, fit.2, fit.3, fit.4, fit.5)## Analysis of Variance Table
##
## Model 1: wage ~ poly(age, 1)
## 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 ' ' 1agelims <- range(Wage$age)
age.grid <- seq(from = agelims[1], to = agelims[2])
fit <- lm(wage ~ poly(age, 4), data = Wage)
preds <- predict(fit, newdata = list(age = age.grid), se = TRUE)
se.bands <- cbind(preds$fit + 2 * preds$se.fit, preds$fit - 2 * preds$se.fit)
plot(Wage$age, Wage$wage, xlim = agelims, cex = 0.5, col = "darkgrey",
xlab = "Age", ylab = "Wage", main = "Degree-4 Polynomial Fit")
lines(age.grid, preds$fit, lwd = 2, col = "blue")
matlines(age.grid, se.bands, lwd = 1, col = "blue", lty = 3)
## CV selected degree 9 as the numerical minimum, but the CV errors from
degree 3 onward were essentially flat (differences of ~0.3–3 units out
of ~1594–1598, well within noise). This is consistent with ANOVA’s
conclusion: the real improvement happens up through degree 3 (or 4), and
everything beyond that is just fitting noise CV’s “optimal” degree 9
isn’t meaningfully better than degree 3, it’s just where random CV-fold
variability happened to bottom out. In practice, degree 3 or 4 is the
more defensible, parsimonious choice, and both methods agree on this
once you look past CV’s raw minimum to its shape.
#set.seed(1)
cv.errors <- rep(NA, 9)
for (c in 2:10) {
Wage$age.cut <- cut(Wage$age, c)
glm.fit <- glm(wage ~ age.cut, data = Wage)
cv.errors[c - 1] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}
cv.errors## [1] 1735.201 1682.307 1636.967 1631.602 1625.438 1613.727 1600.229 1612.512
## [9] 1607.027which.min(cv.errors) + 1## [1] 8fit <- glm(wage ~ cut(age, 8), data = Wage)
preds <- predict(fit, newdata = list(age = age.grid), se = TRUE)
se.bands <- cbind(preds$fit + 2 * preds$se.fit, preds$fit - 2 * preds$se.fit)
plot(Wage$age, Wage$wage, xlim = agelims, cex = 0.5, col = "darkgrey",
xlab = "Age", ylab = "Wage", main = "Step Function Fit (8 cuts)")
lines(age.grid, preds$fit, lwd = 2, col = "red")
matlines(age.grid, se.bands, lwd = 1, col = "red", lty = 3)
#Cross-validation selected degree 9 as the polynomial with the lowest CV error, but the CV errors were nearly flat from degree 3 onward, suggesting this “minimum” is mostly noise rather than a real improvement. ANOVA confirmed this: the jump from degree 1→2 and 2→3 was highly significant, but degree 3→4 was only borderline (p≈0.05) and 4→5 was not significant at all. Together, both methods suggest a degree 3 or 4 polynomial is sufficient to capture the relationship between age and wage, with higher-order terms adding little real value. The degree-4 fit shows wage rising through the 20s–30s, peaking in middle age, and declining toward age 80.
(b): Cross-validation selected 8 cuts as optimal for the step function (CV error ≈ 1600.2), though 6–7 cuts performed nearly as well, again suggesting the CV curve is fairly flat once enough cuts are used. The resulting step function captures the same broad age-wage pattern as the polynomial fit rising through the 20s–30s, peaking/plateauing in middle age, then declining toward age 80 but as a piecewise-constant approximation with visible jumps at each bin boundary rather than a smooth curve.
#QUESTION 10 :
library(ISLR2)
library(leaps)
data(College)
set.seed(1)
train <- sample(1:nrow(College), nrow(College) * 0.7)
test <- -train
College.train <- College[train, ]
College.test <- College[test, ]
fwd.fit <- regsubsets(Outstate ~ ., data = College.train, nvmax = 17, method = "forward")
fwd.summary <- summary(fwd.fit)
fwd.summary## Subset selection object
## Call: regsubsets.formula(Outstate ~ ., data = College.train, 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 ) "*" "*" "*"par(mfrow = c(1, 3))
plot(fwd.summary$cp, xlab = "Number of Variables", ylab = "Cp", type = "b")
points(which.min(fwd.summary$cp), fwd.summary$cp[which.min(fwd.summary$cp)], col = "red", pch = 20, cex = 2)
plot(fwd.summary$bic, xlab = "Number of Variables", ylab = "BIC", type = "b")
points(which.min(fwd.summary$bic), fwd.summary$bic[which.min(fwd.summary$bic)], col = "red", pch = 20, cex = 2)
plot(fwd.summary$adjr2, xlab = "Number of Variables", ylab = "Adjusted R2", type = "b")
points(which.max(fwd.summary$adjr2), fwd.summary$adjr2[which.max(fwd.summary$adjr2)], col = "red", pch = 20, cex = 2)
which.min(fwd.summary$cp)## [1] 13which.min(fwd.summary$bic)## [1] 13which.max(fwd.summary$adjr2)## [1] 13coef(fwd.fit, 13)## (Intercept) PrivateYes Apps Accept Top10perc
## -1739.5725417 2276.7996721 -0.3358567 0.7814587 28.9687655
## F.Undergrad Room.Board Personal PhD Terminal
## -0.1559550 0.9134134 -0.3484815 11.8113175 24.9233138
## S.F.Ratio perc.alumni Expend Grad.Rate
## -55.0649149 48.6046652 0.1744677 20.9498491#Forward stepwise selection was performed on the training set, and Cp, BIC, and adjusted R² all agreed unanimously that the optimal model size is 13 predictors. The selected variables are: Private, Apps, Accept, Top10perc, F.Undergrad, Room.Board, Personal, PhD, Terminal, S.F.Ratio, perc.alumni, Expend, and Grad.Rate. This gives us a satisfactory, more parsimonious subset (13 of 17 available predictors) to carry forward into the GAM.
library(gam)## Loading required package: splines## Loading required package: foreach## Loaded gam 1.22-7gam.fit <- gam(Outstate ~ Private + s(Apps, df = 4) + s(Accept, df = 4) +
s(Top10perc, df = 4) + s(F.Undergrad, df = 4) +
s(Room.Board, df = 4) + s(Personal, df = 4) + s(PhD, df = 4) +
s(Terminal, df = 4) + s(S.F.Ratio, df = 4) +
s(perc.alumni, df = 4) + s(Expend, df = 4) + s(Grad.Rate, df = 4),
data = College.train)
graphics.off()
par(
mfrow = c(3, 3),
mar = c(2.5, 2.5, 2, 1),
mgp = c(1.6, 0.5, 0)
)
plot(gam.fit, se = TRUE, col = "blue")#The GAM reveals that several predictors most notably Expend have a clearly non-linear relationship with out-of-state tuition, while others like Room.Board, perc.alumni, and Accept appear to relate roughly linearly. Predictors like Apps, Personal, PhD, and Terminal show weak or highly uncertain relationships (wide confidence bands), suggesting they may contribute relatively little unique explanatory power once the other variables are accounted for. Private status also has a strong, simple effect private schools charge substantially more out-of-state tuition than public schools.
gam.pred <- predict(gam.fit, newdata = College.test)
gam.err <- mean((College.test$Outstate - gam.pred)^2)
gam.err## [1] 3146155gam.tss <- mean((College.test$Outstate - mean(College.test$Outstate))^2)
gam.rss <- mean((College.test$Outstate - gam.pred)^2)
test.rsq <- 1 - (gam.rss / gam.tss)
test.rsq## [1] 0.7720737#The GAM performs quite well, explaining roughly 77% of the variance in out-of-state tuition using just the 13 predictors selected via forward stepwise selection. This is a solid result for a real-world dataset with substantial institution-to-institution variability, and it demonstrates that allowing for non-linear relationships (via the smoothing splines) captures meaningful structure beyond what a simple linear model would particularly given the strong non-linear pattern we saw for Expend. The GAM’s flexibility appears well-justified rather than overfit, since a reasonably high R² held up on the test set (not just training data)
summary(gam.fit)##
## Call: gam(formula = Outstate ~ Private + s(Apps, df = 4) + s(Accept,
## df = 4) + s(Top10perc, df = 4) + s(F.Undergrad, df = 4) +
## s(Room.Board, df = 4) + s(Personal, df = 4) + s(PhD, df = 4) +
## s(Terminal, df = 4) + s(S.F.Ratio, df = 4) + s(perc.alumni,
## df = 4) + s(Expend, df = 4) + s(Grad.Rate, df = 4), data = College.train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5910.22 -1088.99 73.01 1113.20 7156.97
##
## (Dispersion Parameter for gaussian family taken to be 3195196)
##
## Null Deviance: 9260683704 on 542 degrees of freedom
## Residual Deviance: 1575229760 on 492.9994 degrees of freedom
## AIC: 9723.111
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 2531384624 2531384624 792.2471 < 2.2e-16 ***
## s(Apps, df = 4) 1 920829741 920829741 288.1920 < 2.2e-16 ***
## s(Accept, df = 4) 1 111647663 111647663 34.9424 6.346e-09 ***
## s(Top10perc, df = 4) 1 1058267443 1058267443 331.2058 < 2.2e-16 ***
## s(F.Undergrad, df = 4) 1 273900647 273900647 85.7226 < 2.2e-16 ***
## s(Room.Board, df = 4) 1 513712100 513712100 160.7764 < 2.2e-16 ***
## s(Personal, df = 4) 1 54852808 54852808 17.1673 4.028e-05 ***
## s(PhD, df = 4) 1 61155825 61155825 19.1399 1.483e-05 ***
## s(Terminal, df = 4) 1 22024889 22024889 6.8931 0.0089216 **
## s(S.F.Ratio, df = 4) 1 89787883 89787883 28.1009 1.741e-07 ***
## s(perc.alumni, df = 4) 1 139078026 139078026 43.5272 1.081e-10 ***
## s(Expend, df = 4) 1 468470487 468470487 146.6171 < 2.2e-16 ***
## s(Grad.Rate, df = 4) 1 37708950 37708950 11.8018 0.0006417 ***
## Residuals 493 1575229760 3195196
## ---
## 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(Apps, df = 4) 3 2.1126 0.097725 .
## s(Accept, df = 4) 3 7.7040 4.868e-05 ***
## s(Top10perc, df = 4) 3 1.9302 0.123747
## s(F.Undergrad, df = 4) 3 2.0602 0.104615
## s(Room.Board, df = 4) 3 1.8472 0.137676
## s(Personal, df = 4) 3 2.7447 0.042544 *
## s(PhD, df = 4) 3 1.9813 0.115856
## s(Terminal, df = 4) 3 1.8216 0.142254
## s(S.F.Ratio, df = 4) 3 4.8727 0.002386 **
## s(perc.alumni, df = 4) 3 1.3279 0.264519
## s(Expend, df = 4) 3 25.8833 1.443e-15 ***
## s(Grad.Rate, df = 4) 3 1.3980 0.242634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#There is clear statistical evidence of a non-linear relationship with out-of-state tuition for Expend, Accept, S.F.Ratio, and Personal with Expend showing by far the strongest and most visually obvious non-linear (diminishing-returns) pattern. The remaining predictors (Apps, Top10perc, F.Undergrad, Room.Board, PhD, Terminal, perc.alumni, Grad.Rate) do not show statistically significant non-linearity, meaning a simple linear term would likely be adequate for them despite some visual “wiggliness” in a few of the plots (e.g., PhD, Terminal) those wiggles are within the range of what could occur by chance given the wide confidence bands.