#Chapter 7 question 6
library(ISLR2)
## Warning: package 'ISLR2' was built under R version 4.5.3
library(boot)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
data(Wage)
set.seed(1)
cv.error <- rep(NA, 10)
for(i in 1:10){
fit <- glm(wage ~ poly(age, i), data = Wage)
cv.error[i] <- cv.glm(Wage,
fit,
K = 10)$delta[1]
}
cv.error
## [1] 1676.826 1600.763 1598.399 1595.651 1594.977 1596.061 1594.298 1598.134
## [9] 1593.913 1595.950
cv_df <- data.frame(
Degree = 1:10,
CV_Error = cv.error
)
ggplot(cv_df,
aes(Degree, CV_Error)) +
geom_line(color = "#2E86C1",
linewidth = 1.2) +
geom_point(size = 3,
color = "#1B4F72") +
labs(title = "10-Fold Cross-Validation Errors",
x = "Polynomial Degree",
y = "CV Error") +
theme_minimal(base_size = 14)
best.degree <- which.min(cv.error)
best.degree
## [1] 9
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
#ecause the differences in CV error are very small after degree 4, the 9th-degree model happens to have the lowest prediction error, even though the additional terms are not statistically significant individually.
final.poly <- lm(wage ~ poly(age, best.degree),
data = Wage)
age.grid <- data.frame(
age = seq(min(Wage$age),
max(Wage$age),
length = 200)
)
age.grid$fit <- predict(final.poly,
newdata = age.grid)
ggplot(Wage,
aes(age, wage)) +
geom_point(alpha = 0.30,
color = "gray60") +
geom_line(data = age.grid,
aes(age, fit),
color = "#D35400",
linewidth = 1.5) +
labs(title = paste("Polynomial Regression (Degree",
best.degree,
")"),
x = "Age",
y = "Wage") +
theme_minimal(base_size = 14)
#question 10 of chapter 7
library(ISLR2)
library(leaps)
## Warning: package 'leaps' was built under R version 4.5.3
data(College)
set.seed(1)
train <- sample(
1:nrow(College),
nrow(College) / 2
)
College.train <- College[train, ]
College.test <- College[-train, ]
fit <- regsubsets(
Outstate ~ .,
data = College.train,
nvmax = 17,
method = "forward"
)
fit.summary <- summary(fit)
par(mfrow = c(1, 3))
plot(
fit.summary$cp,
xlab = "Number of Variables",
ylab = "Cp",
type = "l"
)
min.cp <- min(fit.summary$cp)
std.cp <- sd(fit.summary$cp)
abline(
h = min.cp + 0.2 * std.cp,
col = "red",
lty = 2
)
abline(
h = min.cp - 0.2 * std.cp,
col = "red",
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 = "red",
lty = 2
)
abline(
h = min.bic - 0.2 * std.bic,
col = "red",
lty = 2
)
plot(
fit.summary$adjr2,
xlab = "Number of Variables",
ylab = "Adjusted R-squared",
type = "l",
ylim = c(0.4, 0.84)
)
max.adjr2 <- max(fit.summary$adjr2)
std.adjr2 <- sd(fit.summary$adjr2)
abline(
h = max.adjr2 + 0.2 * std.adjr2,
col = "red",
lty = 2
)
abline(
h = max.adjr2 - 0.2 * std.adjr2,
col = "red",
lty = 2
)
par(mfrow = c(1, 1))
coeffs <- coef(fit, id = 6)
coeffs
## (Intercept) PrivateYes Room.Board Terminal perc.alumni
## -4726.8810613 2717.7019276 1.1032433 36.9990286 59.0863753
## Expend Grad.Rate
## 0.1930814 33.8303314
names(coeffs)
## [1] "(Intercept)" "PrivateYes" "Room.Board" "Terminal" "perc.alumni"
## [6] "Expend" "Grad.Rate"
library(gam)
## Warning: package 'gam' was built under R version 4.5.3
## Loading required package: splines
## Loading required package: foreach
## Warning: package 'foreach' was built under R version 4.5.3
## Loaded gam 1.22-7
gam1 <- 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(
gam1,
se = TRUE,
col = "blue"
)
par(mfrow = c(1, 1))
preds <- predict(
gam1,
newdata = College.test
)
# Test mean squared error
test.mse <- mean(
(College.test$Outstate - preds)^2
)
# Test R-squared
tss <- mean(
(College.test$Outstate -
mean(College.test$Outstate))^2
)
test.r2 <- 1 - test.mse / tss
test.mse
## [1] 3349290
test.r2
## [1] 0.7660016
summary(gam1)
##
## 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