STA 6543 Predictive Modeling
Homework 6
Haniya Farhat
July 16, 2026
Problem 6: Polynomial Regression and Step Functions on the Wage Data Set
## Dimensions: 3000 rows x 11 columns(a) Polynomial regression: choosing the degree by cross-validation
set.seed(1)
cv_errs <- rep(NA, 10)
for (d in 1:10) {
fit_d <- glm(wage ~ poly(age, d), data = Wage)
cv_errs[d] <- cv.glm(Wage, fit_d, K = 10)$delta[1]
}
round(cv_errs, 2)## [1] 1676.83 1600.76 1598.40 1595.65 1594.98 1596.06 1594.30 1598.13 1593.91
## [10] 1595.95plot(1:10, cv_errs, type = "b", pch = 19,
xlab = "Degree", ylab = "10-fold CV MSE",
main = "Polynomial Regression: CV Error by Degree")
## Degree that minimizes CV error: 9The raw CV minimum falls at degree 9, but the curve is essentially flat from degree 3 onward — every degree from 3 to 10 lands within about half a unit of MSE of each other, a difference that isn’t meaningful. That flatness lines up with what a sequential ANOVA F-test says:
fit_1 <- lm(wage ~ age, 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 ~ 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 ' ' 1Going from linear to quadratic and quadratic to cubic are both highly significant. Cubic to quartic is borderline (p ≈ 0.051), and quartic to quintic is clearly not significant (p ≈ 0.37). So both the CV curve’s flatness and the ANOVA test point to the same conclusion even though the raw CV minimum technically sits at a higher degree: degree 3 or 4 captures essentially all of the real signal, and any improvement past that is noise. I’ll go with degree 4 as the final choice, matching the borderline ANOVA cutoff and giving a small amount of extra flexibility without meaningfully overfitting.
best_d <- 4
age_grid <- seq(min(Wage$age), max(Wage$age), length.out = 100)
fit_poly <- lm(wage ~ poly(age, best_d), data = Wage)
pred_poly <- predict(fit_poly, newdata = list(age = age_grid), se = TRUE)
band <- cbind(pred_poly$fit - 2 * pred_poly$se.fit,
pred_poly$fit + 2 * pred_poly$se.fit)
plot(wage ~ age, data = Wage, col = "darkgrey", cex = 0.5,
xlab = "Age", ylab = "Wage",
main = paste0("Degree-", best_d, " Polynomial Fit"))
lines(age_grid, pred_poly$fit, col = "blue", lwd = 2)
matlines(age_grid, band, col = "blue", lty = 2)
(b) Step function: choosing the number of cuts by cross-validation
set.seed(1)
cv_errs_step <- rep(NA, 10)
for (cts in 2:10) {
Wage$age.cut <- cut(Wage$age, cts)
fit_cts <- glm(wage ~ age.cut, data = Wage)
cv_errs_step[cts] <- cv.glm(Wage, fit_cts, K = 10)$delta[1]
}
round(cv_errs_step[2:10], 2)## [1] 1734.49 1684.27 1635.55 1632.08 1623.42 1615.00 1601.32 1613.95 1606.33plot(2:10, cv_errs_step[2:10], type = "b", pch = 19,
xlab = "Number of Cuts", ylab = "10-fold CV MSE",
main = "Step Function: CV Error by Number of Cuts")
best_cuts <- which.min(cv_errs_step)
cat("Number of cuts that minimizes CV error:", best_cuts, "\n")## Number of cuts that minimizes CV error: 8Here the CV curve has a clean minimum at 8 cuts, so I’ll use that directly.
fit_step <- glm(wage ~ cut(age, best_cuts), data = Wage)
pred_step <- predict(fit_step, newdata = list(age = age_grid), se.fit = TRUE)
band_step <- cbind(pred_step$fit - 2 * pred_step$se.fit,
pred_step$fit + 2 * pred_step$se.fit)
plot(wage ~ age, data = Wage, col = "darkgrey", cex = 0.5,
xlab = "Age", ylab = "Wage",
main = paste0("Step Function Fit (", best_cuts, " Cuts)"))
lines(age_grid, pred_step$fit, col = "red", lwd = 2)
matlines(age_grid, band_step, col = "red", lty = 2)
The step function fit tells the same basic story as the polynomial: wage rises from the 20s through the 40s-50s, then flattens and drifts down toward retirement age. With 8 cuts the fit has enough resolution to pick up that rise-then-plateau shape without carving the age range into so many pieces that individual bins become noisy.
Problem 10: Predicting Out-of-State Tuition in the College Data Set
data(College, package = "ISLR2")
cat("Dimensions:", nrow(College), "rows x", ncol(College), "columns\n")## Dimensions: 777 rows x 18 columns(a) Forward stepwise selection
set.seed(1)
n_col <- nrow(College)
train10 <- sample(n_col, n_col * 0.7)
col_tr <- College[train10, ]
col_te <- College[-train10, ]
cat("Training:", nrow(col_tr), " Test:", nrow(col_te), "\n")## Training: 543 Test: 234fwd_fit <- regsubsets(Outstate ~ ., data = col_tr, nvmax = 17, method = "forward")
fwd_sum <- summary(fwd_fit)
par(mfrow = c(1, 3))
plot(fwd_sum$cp, xlab = "Number of Variables", ylab = "Cp", type = "b", pch = 19)
points(which.min(fwd_sum$cp), min(fwd_sum$cp), col = "red", pch = 19, cex = 1.5)
plot(fwd_sum$bic, xlab = "Number of Variables", ylab = "BIC", type = "b", pch = 19)
points(which.min(fwd_sum$bic), min(fwd_sum$bic), col = "red", pch = 19, cex = 1.5)
plot(fwd_sum$adjr2, xlab = "Number of Variables", ylab = "Adjusted R2", type = "b", pch = 19)
points(which.max(fwd_sum$adjr2), max(fwd_sum$adjr2), col = "red", pch = 19, cex = 1.5)
The global optimum on all three criteria technically falls around 13 variables, but the curves make clear that most of the gain happens in the first 6 variables — Cp and BIC both flatten out sharply after that point, and adjusted \(R^2\) barely moves from 6 variables onward. Since the goal here is “a satisfactory model using just a subset of predictors,” I’ll use the elbow of the curves rather than the strict optimum and go with the 6-variable model:
## (Intercept) PrivateYes Room.Board PhD perc.alumni
## -3782.6544026 2808.0522815 0.9722009 38.0615210 59.1918276
## Expend Grad.Rate
## 0.2032837 28.6598253That gives Private, Room.Board,
PhD, perc.alumni, Expend, and
Grad.Rate as the predictors carried forward into the
GAM.
(b) Fitting a GAM
gam_fit <- gam(Outstate ~ Private + s(Room.Board, df = 5) + s(PhD, df = 5) +
s(perc.alumni, df = 5) + s(Expend, df = 5) + s(Grad.Rate, df = 5),
data = col_tr)
summary(gam_fit)##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 5) + s(PhD,
## df = 5) + s(perc.alumni, df = 5) + s(Expend, df = 5) + s(Grad.Rate,
## df = 5), data = col_tr)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7589.84 -1158.20 30.67 1238.92 7465.43
##
## (Dispersion Parameter for gaussian family taken to be 3560371)
##
## Null Deviance: 9260683704 on 542 degrees of freedom
## Residual Deviance: 1837148606 on 515.9992 degrees of freedom
## AIC: 9760.632
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 2424267855 2424267855 680.903 < 2.2e-16 ***
## s(Room.Board, df = 5) 1 1724227100 1724227100 484.283 < 2.2e-16 ***
## s(PhD, df = 5) 1 603100356 603100356 169.393 < 2.2e-16 ***
## s(perc.alumni, df = 5) 1 452642481 452642481 127.133 < 2.2e-16 ***
## s(Expend, df = 5) 1 760908252 760908252 213.716 < 2.2e-16 ***
## s(Grad.Rate, df = 5) 1 100712369 100712369 28.287 1.561e-07 ***
## Residuals 516 1837148606 3560371
## ---
## 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 = 5) 4 2.7346 0.02836 *
## s(PhD, df = 5) 4 1.1911 0.31373
## s(perc.alumni, df = 5) 4 1.1891 0.31462
## s(Expend, df = 5) 4 26.3798 < 2e-16 ***
## s(Grad.Rate, df = 5) 4 1.7908 0.12929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Private and Room.Board behave close to
linearly — more expensive room and board is associated with steadily
higher out-of-state tuition, and private schools charge more than public
ones across the board. PhD and perc.alumni are
mostly increasing but with a bit of wiggle at the low end where there’s
less data. Grad.Rate looks close to linear with a slight
flattening at the high end. Expend stands out as clearly
non-linear: tuition rises steeply with instructional expenditure at the
low-to-middle range and then levels off — schools that already spend a
lot per student don’t charge proportionally more out-of-state tuition on
top of that.
(c) Evaluating the model on the test set
pred_gam <- predict(gam_fit, newdata = col_te)
tss10 <- sum((col_te$Outstate - mean(col_te$Outstate))^2)
mse_gam <- mean((col_te$Outstate - pred_gam)^2)
r2_gam <- 1 - sum((col_te$Outstate - pred_gam)^2) / tss10
cat("GAM test MSE:", round(mse_gam, 0), " Test R2:", round(r2_gam, 4), "\n")## GAM test MSE: 3187953 Test R2: 0.769lm_fit10 <- lm(Outstate ~ Private + Room.Board + PhD + perc.alumni + Expend + Grad.Rate,
data = col_tr)
pred_lm10 <- predict(lm_fit10, newdata = col_te)
mse_lm10 <- mean((col_te$Outstate - pred_lm10)^2)
r2_lm10 <- 1 - sum((col_te$Outstate - pred_lm10)^2) / tss10
cat("Linear model test MSE:", round(mse_lm10, 0), " Test R2:", round(r2_lm10, 4), "\n")## Linear model test MSE: 3635355 Test R2: 0.7366The GAM explains about 76.9% of the variance in
Outstate on the test set, versus 73.7% for
a plain linear model using the same six predictors — a solid improvement
in test MSE (4.47402^{5} lower, about a 12.3% reduction). That gap is
exactly what the term plots suggested: letting Expend in
particular bend instead of forcing a straight line captures real
curvature in how spending relates to tuition, and the GAM picks that up
without needing any manual transformation.
(d) Evidence of non-linear relationships
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## Private
## s(Room.Board, df = 5) 4 2.7346 0.02836 *
## s(PhD, df = 5) 4 1.1911 0.31373
## s(perc.alumni, df = 5) 4 1.1891 0.31462
## s(Expend, df = 5) 4 26.3798 < 2e-16 ***
## s(Grad.Rate, df = 5) 4 1.7908 0.12929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The “Anova for Nonparametric Effects” table tests each smooth term
against a linear fit. Expend shows overwhelming evidence of
non-linearity (p < 2e-16) — consistent with the sharp bend seen in
its term plot. Room.Board also shows a statistically
significant non-linear component (p ≈ 0.028), though its curve in the
plot is fairly mild. PhD, perc.alumni, and
Grad.Rate all have non-parametric p-values well above 0.05,
so there’s no strong evidence their relationships with
Outstate depart from linear — a straight-line term would
likely fit about as well for those three. Overall, Expend
is the variable where allowing non-linearity clearly pays off, and it’s
the main driver of the GAM’s improvement over the linear model in part
(c).