This document demonstrates Model Selection techniques using Best Subset Selection, Forward Selection, and Backward Selection methods. We will evaluate models using criteria such as Cp, BIC, and Adjusted R².
First, we need to load the necessary package and set up the dataset.
# Setup CRAN Mirror
options(repos = c(CRAN = "https://cran.rstudio.com/"))
# Load required packages
install.packages("leaps")
The downloaded binary packages are in
/var/folders/m3/k788kw6103zdvc0bwpc1lzd00000gn/T//RtmpLhaThL/downloaded_packagesrequire(leaps)We will generate a dataset where Y is a cubic polynomial of X with added noise (eps).
set.seed(1)
X = rnorm(100)
eps = rnorm(100)
Y = 4 + 9 * X - 2 * X^2 + X^3 + eps
# Plot X and Y
plot(X, Y, pch = 20, col = "blue", main = "Scatter plot of X and Y")
We perform Best Subset Selection to find the optimal polynomial degree for the model.
# Best Subset Selection
regfit.full = regsubsets(Y ~ poly(X, 10, raw = TRUE), data = data.frame(Y, X), nvmax = 10)
reg.summary = summary(regfit.full)
# Display summary
reg.summarySubset selection object
Call: regsubsets.formula(Y ~ poly(X, 10, raw = TRUE), data = data.frame(Y,
X), nvmax = 10)
10 Variables (and intercept)
Forced in Forced out
poly(X, 10, raw = TRUE)1 FALSE FALSE
poly(X, 10, raw = TRUE)2 FALSE FALSE
poly(X, 10, raw = TRUE)3 FALSE FALSE
poly(X, 10, raw = TRUE)4 FALSE FALSE
poly(X, 10, raw = TRUE)5 FALSE FALSE
poly(X, 10, raw = TRUE)6 FALSE FALSE
poly(X, 10, raw = TRUE)7 FALSE FALSE
poly(X, 10, raw = TRUE)8 FALSE FALSE
poly(X, 10, raw = TRUE)9 FALSE FALSE
poly(X, 10, raw = TRUE)10 FALSE FALSE
1 subsets of each size up to 10
Selection Algorithm: exhaustive
poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
1 ( 1 ) "*" " "
2 ( 1 ) "*" "*"
3 ( 1 ) "*" "*"
4 ( 1 ) "*" "*"
5 ( 1 ) "*" "*"
6 ( 1 ) "*" "*"
7 ( 1 ) "*" "*"
8 ( 1 ) "*" "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)3 poly(X, 10, raw = TRUE)4
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) "*" " "
4 ( 1 ) "*" " "
5 ( 1 ) "*" " "
6 ( 1 ) "*" " "
7 ( 1 ) "*" " "
8 ( 1 ) "*" "*"
9 ( 1 ) " " "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)5 poly(X, 10, raw = TRUE)6
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) " " " "
4 ( 1 ) "*" " "
5 ( 1 ) "*" "*"
6 ( 1 ) " " " "
7 ( 1 ) "*" "*"
8 ( 1 ) " " "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)7 poly(X, 10, raw = TRUE)8
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) " " " "
4 ( 1 ) " " " "
5 ( 1 ) " " " "
6 ( 1 ) "*" "*"
7 ( 1 ) " " "*"
8 ( 1 ) " " "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)9 poly(X, 10, raw = TRUE)10
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) " " " "
4 ( 1 ) " " " "
5 ( 1 ) " " " "
6 ( 1 ) "*" " "
7 ( 1 ) " " "*"
8 ( 1 ) "*" "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*" We will now plot and identify the best models based on Cp, BIC, and Adjusted R².
# Set up plotting area
par(mfrow = c(1, 3))
# Plot Cp
min.cp = which.min(reg.summary$cp)
plot(reg.summary$cp, xlab = "Number of Poly(X)", ylab = "Best Subset Cp", type = "l", main = "Cp Plot")
points(min.cp, reg.summary$cp[min.cp], col = "red", pch = 4, lwd = 5)
# Plot BIC
min.bic = which.min(reg.summary$bic)
plot(reg.summary$bic, xlab = "Number of Poly(X)", ylab = "Best Subset BIC", type = "l", main = "BIC Plot")
points(min.bic, reg.summary$bic[min.bic], col = "red", pch = 4, lwd = 5)
# Plot Adjusted R²
max.adjr2 = which.max(reg.summary$adjr2)
plot(reg.summary$adjr2, xlab = "Number of Poly(X)", ylab = "Best Subset Adjusted R²", type = "l", main = "Adjusted R² Plot")
points(max.adjr2, reg.summary$adjr2[max.adjr2], col = "red", pch = 4, lwd = 5)
We can extract the coefficients of the best models based on Cp, BIC, and Adjusted R².
# Coefficients for best models
coef(regfit.full, min.cp) # Best model by Cp (Intercept) poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
4.07200775 9.38745596 -2.15424359
poly(X, 10, raw = TRUE)3 poly(X, 10, raw = TRUE)5
0.55797426 0.08072292 coef(regfit.full, min.bic) # Best model by BIC (Intercept) poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
4.061507 8.975280 -2.123791
poly(X, 10, raw = TRUE)3
1.017639 coef(regfit.full, max.adjr2) # Best model by Adjusted R² (Intercept) poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
4.07200775 9.38745596 -2.15424359
poly(X, 10, raw = TRUE)3 poly(X, 10, raw = TRUE)5
0.55797426 0.08072292 We will now perform forward and backward selection to compare with the best subset selection.
# Forward Selection
regfit.fwd = regsubsets(Y ~ poly(X, 10, raw = TRUE), data = data.frame(Y, X), nvmax = 10, method = "forward")
fwd.summary = summary(regfit.fwd)
# Display summary
fwd.summarySubset selection object
Call: regsubsets.formula(Y ~ poly(X, 10, raw = TRUE), data = data.frame(Y,
X), nvmax = 10, method = "forward")
10 Variables (and intercept)
Forced in Forced out
poly(X, 10, raw = TRUE)1 FALSE FALSE
poly(X, 10, raw = TRUE)2 FALSE FALSE
poly(X, 10, raw = TRUE)3 FALSE FALSE
poly(X, 10, raw = TRUE)4 FALSE FALSE
poly(X, 10, raw = TRUE)5 FALSE FALSE
poly(X, 10, raw = TRUE)6 FALSE FALSE
poly(X, 10, raw = TRUE)7 FALSE FALSE
poly(X, 10, raw = TRUE)8 FALSE FALSE
poly(X, 10, raw = TRUE)9 FALSE FALSE
poly(X, 10, raw = TRUE)10 FALSE FALSE
1 subsets of each size up to 10
Selection Algorithm: forward
poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
1 ( 1 ) "*" " "
2 ( 1 ) "*" "*"
3 ( 1 ) "*" "*"
4 ( 1 ) "*" "*"
5 ( 1 ) "*" "*"
6 ( 1 ) "*" "*"
7 ( 1 ) "*" "*"
8 ( 1 ) "*" "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)3 poly(X, 10, raw = TRUE)4
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) "*" " "
4 ( 1 ) "*" " "
5 ( 1 ) "*" " "
6 ( 1 ) "*" " "
7 ( 1 ) "*" " "
8 ( 1 ) "*" " "
9 ( 1 ) "*" " "
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)5 poly(X, 10, raw = TRUE)6
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) " " " "
4 ( 1 ) "*" " "
5 ( 1 ) "*" "*"
6 ( 1 ) "*" "*"
7 ( 1 ) "*" "*"
8 ( 1 ) "*" "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)7 poly(X, 10, raw = TRUE)8
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) " " " "
4 ( 1 ) " " " "
5 ( 1 ) " " " "
6 ( 1 ) " " " "
7 ( 1 ) "*" " "
8 ( 1 ) "*" "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)9 poly(X, 10, raw = TRUE)10
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) " " " "
4 ( 1 ) " " " "
5 ( 1 ) " " " "
6 ( 1 ) "*" " "
7 ( 1 ) "*" " "
8 ( 1 ) "*" " "
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*" # Backward Selection
regfit.bwd = regsubsets(Y ~ poly(X, 10, raw = TRUE), data = data.frame(Y, X), nvmax = 10, method = "backward")
bwd.summary = summary(regfit.bwd)
# Display summary
bwd.summarySubset selection object
Call: regsubsets.formula(Y ~ poly(X, 10, raw = TRUE), data = data.frame(Y,
X), nvmax = 10, method = "backward")
10 Variables (and intercept)
Forced in Forced out
poly(X, 10, raw = TRUE)1 FALSE FALSE
poly(X, 10, raw = TRUE)2 FALSE FALSE
poly(X, 10, raw = TRUE)3 FALSE FALSE
poly(X, 10, raw = TRUE)4 FALSE FALSE
poly(X, 10, raw = TRUE)5 FALSE FALSE
poly(X, 10, raw = TRUE)6 FALSE FALSE
poly(X, 10, raw = TRUE)7 FALSE FALSE
poly(X, 10, raw = TRUE)8 FALSE FALSE
poly(X, 10, raw = TRUE)9 FALSE FALSE
poly(X, 10, raw = TRUE)10 FALSE FALSE
1 subsets of each size up to 10
Selection Algorithm: backward
poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
1 ( 1 ) "*" " "
2 ( 1 ) "*" "*"
3 ( 1 ) "*" "*"
4 ( 1 ) "*" "*"
5 ( 1 ) "*" "*"
6 ( 1 ) "*" "*"
7 ( 1 ) "*" "*"
8 ( 1 ) "*" "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)3 poly(X, 10, raw = TRUE)4
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) " " " "
4 ( 1 ) " " " "
5 ( 1 ) " " " "
6 ( 1 ) " " " "
7 ( 1 ) " " " "
8 ( 1 ) " " " "
9 ( 1 ) " " "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)5 poly(X, 10, raw = TRUE)6
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) "*" " "
4 ( 1 ) "*" " "
5 ( 1 ) "*" " "
6 ( 1 ) "*" " "
7 ( 1 ) "*" " "
8 ( 1 ) "*" "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)7 poly(X, 10, raw = TRUE)8
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) " " " "
4 ( 1 ) "*" " "
5 ( 1 ) "*" "*"
6 ( 1 ) "*" "*"
7 ( 1 ) "*" "*"
8 ( 1 ) "*" "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*"
poly(X, 10, raw = TRUE)9 poly(X, 10, raw = TRUE)10
1 ( 1 ) " " " "
2 ( 1 ) " " " "
3 ( 1 ) " " " "
4 ( 1 ) " " " "
5 ( 1 ) " " " "
6 ( 1 ) "*" " "
7 ( 1 ) "*" "*"
8 ( 1 ) "*" "*"
9 ( 1 ) "*" "*"
10 ( 1 ) "*" "*" We will now plot the results for forward and backward selection using Cp, BIC, and Adjusted R².
# Set up plotting area
par(mfrow = c(3, 2))
# Forward Selection Cp
min.cp = which.min(fwd.summary$cp)
plot(fwd.summary$cp, xlab = "Number of Poly(X)", ylab = "Forward Selection Cp", type = "l", main = "Forward Cp Plot")
points(min.cp, fwd.summary$cp[min.cp], col = "red", pch = 4, lwd = 5)
# Backward Selection Cp
min.cp = which.min(bwd.summary$cp)
plot(bwd.summary$cp, xlab = "Number of Poly(X)", ylab = "Backward Selection Cp", type = "l", main = "Backward Cp Plot")
points(min.cp, bwd.summary$cp[min.cp], col = "red", pch = 4, lwd = 5)
# Forward Selection BIC
min.bic = which.min(fwd.summary$bic)
plot(fwd.summary$bic, xlab = "Number of Poly(X)", ylab = "Forward Selection BIC", type = "l", main = "Forward BIC Plot")
points(min.bic, fwd.summary$bic[min.bic], col = "red", pch = 4, lwd = 5)
# Backward Selection BIC
min.bic = which.min(bwd.summary$bic)
plot(bwd.summary$bic, xlab = "Number of Poly(X)", ylab = "Backward Selection BIC", type = "l", main = "Backward BIC Plot")
points(min.bic, bwd.summary$bic[min.bic], col = "red", pch = 4, lwd = 5)
# Forward Selection Adjusted R²
min.adjr2 = which.max(fwd.summary$adjr2)
plot(fwd.summary$adjr2, xlab = "Number of Poly(X)", ylab = "Forward Selection Adjusted R²", type = "l", main = "Forward Adjusted R² Plot")
points(min.adjr2, fwd.summary$adjr2[min.adjr2], col = "red", pch = 4, lwd = 5)
# Backward Selection Adjusted R²
min.adjr2 = which.max(bwd.summary$adjr2)
plot(bwd.summary$adjr2, xlab = "Number of Poly(X)", ylab = "Backward Selection Adjusted R²", type = "l", main = "Backward Adjusted R² Plot")
points(min.adjr2, bwd.summary$adjr2[min.adjr2], col = "red", pch = 4, lwd = 5)
Finally, we extract the coefficients for the best models from forward and backward selection.
# Coefficients for best models from forward selection
coef(regfit.fwd, which.min(fwd.summary$cp)) # Best model by Cp (Intercept) poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
4.07200775 9.38745596 -2.15424359
poly(X, 10, raw = TRUE)3 poly(X, 10, raw = TRUE)5
0.55797426 0.08072292 coef(regfit.fwd, which.min(fwd.summary$bic)) # Best model by BIC (Intercept) poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
4.061507 8.975280 -2.123791
poly(X, 10, raw = TRUE)3
1.017639 coef(regfit.fwd, which.max(fwd.summary$adjr2)) # Best model by Adjusted R² (Intercept) poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
4.07200775 9.38745596 -2.15424359
poly(X, 10, raw = TRUE)3 poly(X, 10, raw = TRUE)5
0.55797426 0.08072292 # Coefficients for best models from backward selection
coef(regfit.bwd, which.min(bwd.summary$cp)) # Best model by Cp (Intercept) poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
4.06013965 9.70201369 -2.13582615
poly(X, 10, raw = TRUE)5 poly(X, 10, raw = TRUE)7
0.30330191 -0.02419389 coef(regfit.bwd, which.min(bwd.summary$bic)) # Best model by BIC (Intercept) poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
4.0738073 9.9427063 -2.1784855
poly(X, 10, raw = TRUE)5
0.1721833 coef(regfit.bwd, which.max(bwd.summary$adjr2)) # Best model by Adjusted R² (Intercept) poly(X, 10, raw = TRUE)1 poly(X, 10, raw = TRUE)2
4.06013965 9.70201369 -2.13582615
poly(X, 10, raw = TRUE)5 poly(X, 10, raw = TRUE)7
0.30330191 -0.02419389 In this lab, we have explored Model Selection techniques, including Best Subset Selection, Forward Selection, and Backward Selection. We evaluated models using Cp, BIC, and Adjusted R² criteria and extracted the best model coefficients for each approach.