STA6543 — Assignment 6 (Chapter 7: Moving Beyond Linearity)
Yifan Gao
June 2026
1 Exercise 6 — Polynomial
regression and step functions on Wage
Wage <- na.omit(Wage)
dim(Wage)## [1] 3000 111.1 (a) Polynomial
regression of wage on age
1.1.1 Choosing the degree by 10-fold cross-validation
set.seed(1)
max.d <- 10
cv.err <- rep(NA, max.d)
for (d in 1:max.d) {
fit <- glm(wage ~ poly(age, d), data = Wage)
cv.err[d] <- cv.glm(Wage, fit, K = 10)$delta[1]
}
raw.min <- which.min(cv.err)
# the CV curve is essentially flat beyond degree 3-4: adopt the smallest degree whose
# CV error is within 0.2% of the minimum (a parsimonious "within-noise" choice)
poly.d <- min(which(cv.err <= 1.002 * min(cv.err)))
plot(1:max.d, cv.err, type = "b", xlab = "Degree", ylab = "10-fold CV error",
main = "Polynomial degree selection")
points(raw.min, cv.err[raw.min], col = "grey40", pch = 1)
points(poly.d, cv.err[poly.d], col = "red", pch = 19)
c(raw_argmin = raw.min, chosen_degree = poly.d, cv_error = round(cv.err[poly.d], 2))## raw_argmin chosen_degree cv_error
## 9.00 4.00 1595.65The CV error drops sharply from degree 1 to 2–3 and is then essentially flat: every degree from 3 to 10 lies within about \(0.3\%\) of the overall minimum, so the nominal argmin at degree 9 is not meaningfully better than a low-order fit (the difference is CV noise). Applying a “within-noise” rule — the smallest degree whose CV error is within \(0.2\%\) of the minimum — gives degree 4, a parsimonious choice consistent with the elbow of the curve.
1.1.2 Comparison with ANOVA hypothesis testing
fits <- lapply(1:5, function(d) lm(wage ~ poly(age, d), data = Wage))
do.call(anova, fits)## Analysis of Variance Table
##
## Model 1: wage ~ poly(age, d)
## Model 2: wage ~ poly(age, d)
## Model 3: wage ~ poly(age, d)
## Model 4: wage ~ poly(age, d)
## Model 5: wage ~ poly(age, d)
## 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 ' ' 1The ANOVA of nested models shows that going from degree 1 to 2 and from 2 to 3 is highly significant (\(p < 0.001\)), degree 3 to 4 is marginal (around the \(5\%\) level), and degree 4 to 5 is clearly not significant. ANOVA therefore points to a cubic (degree 3, or at most 4) fit — consistent with the cross-validation result: both approaches agree that a degree-3/4 polynomial captures the age–wage relationship and higher-order terms add nothing.
1.1.3 Plot of the chosen polynomial fit
age.grid <- seq(min(Wage$age), max(Wage$age), length = 100)
fit <- lm(wage ~ poly(age, poly.d), data = Wage)
pred <- predict(fit, newdata = list(age = age.grid), se = TRUE)
se.bands <- cbind(pred$fit + 2 * pred$se.fit, pred$fit - 2 * pred$se.fit)
plot(Wage$age, Wage$wage, col = "darkgrey", cex = 0.5,
xlab = "Age", ylab = "Wage", main = paste0("Degree-", poly.d, " polynomial fit"))
lines(age.grid, pred$fit, col = "blue", lwd = 2)
matlines(age.grid, se.bands, col = "blue", lty = 2)
The fit rises steeply through the twenties and thirties, peaks around age 40–50, and declines afterward.
1.2 (b) Step function of
wage on age
1.2.1 Choosing the number of cuts by 10-fold cross-validation
set.seed(1)
max.cuts <- 10
cv.err <- rep(NA, max.cuts)
for (k in 2:max.cuts) {
Wage$age.cut <- cut(Wage$age, k)
fit <- glm(wage ~ age.cut, data = Wage)
cv.err[k] <- cv.glm(Wage, fit, K = 10)$delta[1]
}
Wage$age.cut <- NULL
best.k <- which.min(cv.err)
plot(2:max.cuts, cv.err[2:max.cuts], type = "b", xlab = "Number of intervals",
ylab = "10-fold CV error", main = "Step-function cut selection")
points(best.k, cv.err[best.k], col = "red", pch = 19)
c(best_cuts = best.k, cv_error = round(cv.err[best.k], 2))## best_cuts cv_error
## 8.00 1601.32Cross-validation selects 8 intervals (i.e. 7 internal cut points). The error falls quickly up to about 7–8 bins and then flattens.
1.2.2 Plot of the chosen step-function fit
fit <- lm(wage ~ cut(age, best.k), data = Wage)
pred <- predict(fit, newdata = list(age = age.grid), se = TRUE)
se.bands <- cbind(pred$fit + 2 * pred$se.fit, pred$fit - 2 * pred$se.fit)
plot(Wage$age, Wage$wage, col = "darkgrey", cex = 0.5,
xlab = "Age", ylab = "Wage", main = paste0("Step function with ", best.k, " bins"))
lines(age.grid, pred$fit, col = "red", lwd = 2)
matlines(age.grid, se.bands, col = "red", lty = 2)
The step function reproduces the same overall shape as the polynomial — low wages for the young, a plateau in middle age, and a decline near retirement — but as a piecewise-constant approximation.
2 Exercise 10 — GAM for
out-of-state tuition in College
College <- na.omit(College)
dim(College)## [1] 777 182.1 (a) Forward stepwise selection on the training set
set.seed(1)
n <- nrow(College)
train <- sample(n, n / 2)
test <- setdiff(seq_len(n), train)
fwd <- regsubsets(Outstate ~ ., data = College[train, ],
nvmax = 17, method = "forward")
fwd.sum <- summary(fwd)
par(mfrow = c(1, 3))
plot(fwd.sum$cp, type = "b", xlab = "# predictors", ylab = "Cp", main = "Cp")
plot(fwd.sum$bic, type = "b", xlab = "# predictors", ylab = "BIC", main = "BIC")
plot(fwd.sum$adjr2,type = "b", xlab = "# predictors", ylab = "Adj R2", main = "Adj R2")
par(mfrow = c(1, 1))
c(min_cp = which.min(fwd.sum$cp),
min_bic = which.min(fwd.sum$bic),
max_adjr2 = which.max(fwd.sum$adjr2))## min_cp min_bic max_adjr2
## 14 6 14Cp and adjusted \(R^2\)
keep improving out to large models, but all three criteria show a sharp
“elbow”: the fit barely improves beyond about six
predictors. For a parsimonious yet satisfactory model we take the best
6-variable model:
coef(fwd, 6)## (Intercept) PrivateYes Room.Board Terminal perc.alumni
## -4726.8810613 2717.7019276 1.1032433 36.9990286 59.0863753
## Expend Grad.Rate
## 0.1930814 33.8303314The six selected predictors are Private,
Room.Board, Terminal,
perc.alumni, Expend, and
Grad.Rate.
2.2 (b) Fit a GAM on the selected predictors and plot
We fit each continuous predictor with a smoothing spline (4 df) and
keep the qualitative Private as a linear term.
gam.fit <- gam(Outstate ~ Private + s(Room.Board, 4) + s(Terminal, 4) +
s(perc.alumni, 4) + s(Expend, 4) + s(Grad.Rate, 4),
data = College[train, ])
par(mfrow = c(2, 3))
plot(gam.fit, se = TRUE, col = "blue")
par(mfrow = c(1, 1))Findings. Out-of-state tuition rises with
Room.Board, Terminal,
perc.alumni, and Grad.Rate, and private
colleges charge substantially more than public ones. The most striking
pattern is Expend: tuition increases steeply with
instructional spending up to about $15{,}000–20{,}000 per student and
then levels off (even bends down) — a clearly
non-linear, saturating relationship.
2.3 (c) Evaluate on the test set
pred <- predict(gam.fit, newdata = College[test, ])
test.mse <- mean((College$Outstate[test] - pred)^2)
tss <- mean((College$Outstate[test] - mean(College$Outstate[test]))^2)
test.r2 <- 1 - test.mse / tss
c(test_MSE = round(test.mse), test_RMSE = round(sqrt(test.mse)), test_R2 = round(test.r2, 3))## test_MSE test_RMSE test_R2
## 3353802.000 1831.000 0.766On the held-out test set the GAM explains about 77% of the variation in out-of-state tuition (test \(R^2 \approx 0.766\)), with a typical prediction error (RMSE) of roughly $1,831. This is a good fit and improves on a purely linear model of the same predictors, reflecting the benefit of allowing non-linear terms.
lm.fit <- lm(Outstate ~ Private + Room.Board + Terminal + perc.alumni + Expend + Grad.Rate,
data = College[train, ])
lm.mse <- mean((College$Outstate[test] - predict(lm.fit, College[test, ]))^2)
c(GAM_test_MSE = round(test.mse), linear_test_MSE = round(lm.mse))## GAM_test_MSE linear_test_MSE
## 3353802 38448572.4 (d) Which variables show evidence of a non-linear relationship?
summary(gam.fit)##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, 4) + s(Terminal,
## 4) + s(perc.alumni, 4) + s(Expend, 4) + s(Grad.Rate, 4),
## data = College[train, ])
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7128.62 -1133.86 -74.25 1231.50 7369.50
##
## (Dispersion Parameter for gaussian family taken to be 3724586)
##
## Null Deviance: 6989966760 on 387 degrees of freedom
## Residual Deviance: 1363197370 on 365.9997 degrees of freedom
## AIC: 6995.069
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 1764398916 1764398916 473.717 < 2.2e-16 ***
## s(Room.Board, 4) 1 1616561254 1616561254 434.024 < 2.2e-16 ***
## s(Terminal, 4) 1 287918343 287918343 77.302 < 2.2e-16 ***
## s(perc.alumni, 4) 1 354690429 354690429 95.230 < 2.2e-16 ***
## s(Expend, 4) 1 601731164 601731164 161.556 < 2.2e-16 ***
## s(Grad.Rate, 4) 1 90312393 90312393 24.248 1.284e-06 ***
## Residuals 366 1363197370 3724586
## ---
## 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, 4) 3 1.9107 0.1274
## s(Terminal, 4) 3 1.4636 0.2241
## s(perc.alumni, 4) 3 0.3498 0.7893
## s(Expend, 4) 3 26.1184 2.442e-15 ***
## s(Grad.Rate, 4) 3 0.9075 0.4375
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The “Anova for Nonparametric Effects” table tests,
for each smooth term, whether a non-linear component is needed beyond
the linear term. Only Expend is
significant (\(p \ll 0.001\)), giving
strong evidence of a non-linear relationship with
tuition — exactly the saturating curve seen in part (b). All the other
smooth terms have large \(p\)-values
(Room.Board \(\approx
0.13\), Terminal \(\approx
0.22\), perc.alumni \(\approx 0.79\), Grad.Rate
\(\approx 0.44\)), so there is
no evidence of non-linearity for them — their effects
on out-of-state tuition are adequately described as linear. In short,
Expend is the one clearly non-linear predictor, and the
rest are effectively linear.