Group Members

Registration Number Name
223014431 Niyobyose Brancal
223010775 Iradukunda Claudine
223021861 Ishimwe Aime Louange
223012050 Ntihinyurwa Ezechiel
223020453 Iradukunda Boni Magnifique 
library(leaps)
library(glmnet)
library(pls)

Question 1

(a) and (b) Generate X, noise, and Y

We choose the true model: Y = 2 + 3X - 2X² + X³ + noise.

set.seed(123)
n <- 100
X <- rnorm(n)
noise <- rnorm(n)
Y <- 2 + 3 * X - 2 * X^2 + X^3 + noise

# Make a data frame with X, X^2, ..., X^10
sim_data <- data.frame(Y = Y)
for (i in 1:10) {
  sim_data[[paste0("X", i)]] <- X^i
}

head(sim_data)
##           Y          X1          X2            X3           X4            X5
## 1 -1.196163 -0.56047565 0.314132950 -0.1760638685 9.867951e-02 -5.530746e-02
## 2  1.448193 -0.23017749 0.052981677 -0.0121951893 2.807058e-03 -6.461216e-04
## 3  5.357283  1.55870831 2.429571609  3.7869934661 5.902818e+00  9.200772e+00
## 4  1.854390  0.07050839 0.004971433  0.0003505278 2.471515e-05  1.742625e-06
## 5  1.404975  0.12928774 0.016715318  0.0021610857 2.794019e-04  3.612324e-05
## 6  6.262046  1.71506499 2.941447909  5.0447743199 8.652116e+00  1.483894e+01
##             X6            X7           X8            X9          X10
## 1 3.099849e-02 -1.737390e-02 9.737646e-03 -5.457713e-03 3.058915e-03
## 2 1.487226e-04 -3.423260e-05 7.879575e-06 -1.813701e-06 4.174731e-07
## 3 1.434132e+01  2.235393e+01 3.484326e+01  5.431048e+01 8.465420e+01
## 4 1.228697e-07  8.663346e-09 6.108386e-10  4.306925e-11 3.036743e-12
## 5 4.670291e-06  6.038114e-07 7.806541e-08  1.009290e-08 1.304888e-09
## 6 2.544975e+01  4.364797e+01 7.485911e+01  1.283882e+02 2.201942e+02

We generated 100 observations. Only X, X², and X³ truly affect Y; the other powers are included to test whether model selection can identify the important predictors.

  1. Best subset selection
best_fit <- regsubsets(Y ~ ., data = sim_data, nvmax = 10,
                       method = "exhaustive")
best_info <- summary(best_fit)

cp_size <- which.min(best_info$cp)
bic_size <- which.min(best_info$bic)
r2_size <- which.max(best_info$adjr2)

# Best number of predictors under each criterion
data.frame(Criterion = c("Cp", "BIC", "Adjusted R2"),
           Number_of_predictors = c(cp_size, bic_size, r2_size))
##     Criterion Number_of_predictors
## 1          Cp                    3
## 2         BIC                    3
## 3 Adjusted R2                    6
# Coefficients of the selected models
cat("\nCoefficients of the Cp model:\n")
## 
## Coefficients of the Cp model:
print(round(coef(best_fit, id = cp_size), 3))
## (Intercept)          X1          X2          X3 
##       1.970       2.920      -2.092       1.020
cat("\nCoefficients of the BIC model:\n")
## 
## Coefficients of the BIC model:
print(round(coef(best_fit, id = bic_size), 3))
## (Intercept)          X1          X2          X3 
##       1.970       2.920      -2.092       1.020
cat("\nCoefficients of the adjusted R2 model:\n")
## 
## Coefficients of the adjusted R2 model:
print(round(coef(best_fit, id = r2_size), 3))
## (Intercept)          X1          X3          X4          X6          X8 
##       1.824       2.896       1.049      -3.608       1.989      -0.438 
##         X10 
##       0.034
par(mfrow = c(1, 3))

plot(1:10, best_info$cp, type = "b", xlab = "Model size",
     ylab = "Cp", main = "Cp (lower is better)")
abline(v = cp_size, col = "red", lty = 2)

plot(1:10, best_info$bic, type = "b", xlab = "Model size",
     ylab = "BIC", main = "BIC (lower is better)")
abline(v = bic_size, col = "red", lty = 2)

plot(1:10, best_info$adjr2, type = "b", xlab = "Model size",
     ylab = "Adjusted R2", main = "Adjusted R2 (higher is better)")
abline(v = r2_size, col = "red", lty = 2)

par(mfrow = c(1, 1))

Cp selects 3 predictors, BIC selects 3, and adjusted R² selects 6. The red lines show the best model sizes. Since our true model uses three powers, models that recover X, X², and X³ are closest to the data-generating process. Different criteria may select different extra terms because of random noise.

  1. Forward and backward stepwise selection
forward_fit <- regsubsets(Y ~ ., data = sim_data, nvmax = 10,
                          method = "forward")
backward_fit <- regsubsets(Y ~ ., data = sim_data, nvmax = 10,
                           method = "backward")

forward_info <- summary(forward_fit)
backward_info <- summary(backward_fit)

comparison8 <- data.frame(
  Method = c("Best subset", "Forward", "Backward"),
  Cp = c(cp_size, which.min(forward_info$cp),
         which.min(backward_info$cp)),
  BIC = c(bic_size, which.min(forward_info$bic),
          which.min(backward_info$bic)),
  Adjusted_R2 = c(r2_size, which.max(forward_info$adjr2),
                  which.max(backward_info$adjr2))
)
print(comparison8)
##        Method Cp BIC Adjusted_R2
## 1 Best subset  3   3           6
## 2     Forward  3   3           4
## 3    Backward  6   6           7
cat("\nForward model coefficients (chosen by BIC):\n")
## 
## Forward model coefficients (chosen by BIC):
print(round(coef(forward_fit, id = which.min(forward_info$bic)), 3))
## (Intercept)          X1          X2          X3 
##       1.970       2.920      -2.092       1.020
cat("\nBackward model coefficients (chosen by BIC):\n")
## 
## Backward model coefficients (chosen by BIC):
print(round(coef(backward_fit, id = which.min(backward_info$bic)), 3))
## (Intercept)          X1          X4          X5          X6          X7 
##       1.763       3.356      -2.612       0.553       0.936      -0.078 
##          X8 
##      -0.100

Compare the three rows. Best subset tests all predictor combinations, while forward selection adds predictors and backward selection removes them. Similar selected model sizes suggest that the methods are giving consistent results; differences are also possible.

  1. Lasso regression with cross-validation
x_matrix <- as.matrix(sim_data[, -1])

set.seed(123)
lasso_cv8 <- cv.glmnet(x_matrix, Y, alpha = 1, nfolds = 10)
#plot(lasso_cv8)

cat("Best lambda:", round(lasso_cv8$lambda.min, 5), "\n")
## Best lambda: 0.05072
cat("Lasso coefficients:\n")
## Lasso coefficients:
print(round(as.matrix(coef(lasso_cv8, s = "lambda.min")), 3))
##             lambda.min
## (Intercept)      1.930
## X1               2.888
## X2              -2.035
## X3               1.005
## X4               0.000
## X5               0.000
## X6               0.000
## X7               0.000
## X8               0.000
## X9               0.000
## X10              0.000

Cross-validation selects lambda = 0.0507 because it gives the smallest average validation error. Lasso shrinks coefficients and can make unimportant ones exactly zero. Ideally, the largest useful effects correspond to X, X², and X³.

  1. Repeat using Y = 2 + 3X⁷ + noise
Y7 <- 2 + 3 * X^7 + noise
sim_data7 <- sim_data
sim_data7$Y <- Y7

# Best subset selection (BIC)
best_fit7 <- regsubsets(Y ~ ., data = sim_data7, nvmax = 10)
best_size7 <- which.min(summary(best_fit7)$bic)
cat("Best subset model size (BIC):", best_size7, "\n")
## Best subset model size (BIC): 1
print(round(coef(best_fit7, id = best_size7), 3))
## (Intercept)          X7 
##       1.893       3.000
# Lasso selection
set.seed(123)
lasso_cv7 <- cv.glmnet(x_matrix, Y7, alpha = 1, nfolds = 10)
#plot(lasso_cv7)
cat("Lasso best lambda:", round(lasso_cv7$lambda.min, 5), "\n")
## Lasso best lambda: 4.58664
print(round(as.matrix(coef(lasso_cv7, s = "lambda.min")), 3))
##             lambda.min
## (Intercept)      2.129
## X1               0.000
## X2               0.000
## X3               0.000
## X4               0.000
## X5               0.000
## X6               0.000
## X7               2.912
## X8               0.000
## X9               0.000
## X10              0.000

Here X⁷ is the only true predictor (apart from the intercept). A successful method should give an important role to X⁷. Other powers may also enter because polynomial terms are correlated, especially for high powers of X.

Question 2

  1. Compare best subset, ridge, lasso, and PCR

We predict crim (per capita crime rate). We train each model using 70% of the observations and evaluate it on the remaining 30%. Lower test mean squared error (MSE) means better prediction on this split.

data("Boston", package = "MASS")

set.seed(321)
train_rows <- sample(1:nrow(Boston), size = round(0.7 * nrow(Boston)))
train <- Boston[train_rows, ]
test <- Boston[-train_rows, ]

cat("Training observations:", nrow(train), "\n")
## Training observations: 354
cat("Test observations:", nrow(test), "\n")
## Test observations: 152

Best subset selection

subset_fit <- regsubsets(crim ~ ., data = train,
                         nvmax = ncol(train) - 1)
subset_size <- which.min(summary(subset_fit)$bic)
subset_vars <- names(coef(subset_fit, id = subset_size))[-1]

# Refit the selected variables so we can predict new observations
subset_lm <- lm(reformulate(subset_vars, response = "crim"), data = train)
subset_pred <- predict(subset_lm, newdata = test)

cat("Selected predictors:", paste(subset_vars, collapse = ", "), "\n")
## Selected predictors: rad, black, lstat

Ridge and lasso regression

x_train <- model.matrix(crim ~ ., data = train)[, -1]
x_test <- model.matrix(crim ~ ., data = test)[, -1]

y_train <- train$crim

set.seed(321)
ridge_cv <- cv.glmnet(x_train, y_train, alpha = 0, nfolds = 10)
set.seed(321)
lasso_cv <- cv.glmnet(x_train, y_train, alpha = 1, nfolds = 10)

ridge_pred <- as.numeric(predict(ridge_cv, newx = x_test,
                                 s = "lambda.min"))
lasso_pred <- as.numeric(predict(lasso_cv, newx = x_test,
                                 s = "lambda.min"))

cat("Ridge lambda:", round(ridge_cv$lambda.min, 4), "\n")
## Ridge lambda: 0.4979
cat("Lasso lambda:", round(lasso_cv$lambda.min, 4), "\n")
## Lasso lambda: 0.0206

Principal components regression (PCR)

set.seed(321)
pcr_fit <- pcr(crim ~ ., data = train, scale = TRUE,
               validation = "CV", segments = 10)

# Choose the number of components with the lowest CV error.
# The first position is the intercept-only model, so we skip it.
pcr_errors <- as.numeric(RMSEP(pcr_fit, estimate = "CV")$val[1, 1, -1])
best_components <- which.min(pcr_errors)

pcr_pred <- as.numeric(predict(pcr_fit, newdata = test,
                               ncomp = best_components))
cat("Selected PCR components:", best_components, "\n")
## Selected PCR components: 13
  1. Compare test MSE and choose a model
mse <- function(actual, predicted) {
  mean((actual - predicted)^2)
}

results <- data.frame(
  Method = c("Best subset", "Ridge", "Lasso", "PCR"),
  Test_MSE = c(mse(test$crim, subset_pred),
               mse(test$crim, ridge_pred),
               mse(test$crim, lasso_pred),
               mse(test$crim, pcr_pred))
)

results <- results[order(results$Test_MSE), ]
print(transform(results, Test_MSE = round(Test_MSE, 3)), row.names = FALSE)
##       Method Test_MSE
##          PCR   76.724
##        Lasso   76.844
##        Ridge   78.029
##  Best subset   78.293
best_method <- results$Method[1]
cat("Model with lowest test MSE:", best_method, "\n")
## Model with lowest test MSE: PCR
library(ggplot2)

ggplot(results, aes(x = reorder(Method, Test_MSE), y = Test_MSE, fill = Method)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = round(Test_MSE, 3)), vjust = -0.2, size = 4) +
  labs(
    title = "Boston crime: model comparison",
    x = "Method",
    y = "Test MSE"
  ) +
  theme_minimal() +
  expand_limits(y = max(results$Test_MSE) * 1.15)

PCR has the smallest test MSE (76.724) in this train/test split, so it is our preferred model here. The comparison is based on unseen test observations, not the training error. Another random split could give slightly different results.

  1. Does the chosen model use all predictors?
all_vars <- colnames(x_train)

if (best_method == "Best subset") {
  used_vars <- subset_vars
} else if (best_method == "Lasso") {
  lasso_coefs <- as.matrix(coef(lasso_cv, s = "lambda.min"))
  used_vars <- rownames(lasso_coefs)[lasso_coefs[, 1] != 0]
  used_vars <- setdiff(used_vars, "(Intercept)")
} else {
  # Ridge uses all predictors; PCR builds components from all predictors.
  used_vars <- all_vars
}

cat("Number of predictors used:", length(used_vars),
    "out of", length(all_vars), "\n")
## Number of predictors used: 13 out of 13
cat("Predictors:", paste(used_vars, collapse = ", "), "\n")
## Predictors: zn, indus, chas, nox, rm, age, dis, rad, tax, ptratio, black, lstat, medv

The chosen method uses 13 out of 13 original predictors. Best subset includes only selected variables, and lasso can remove variables by setting coefficients to zero. Ridge keeps all predictors, while PCR forms components using the original predictors (even if only some components are retained). A smaller model can be easier to explain, but the main goal here is a low test MSE.