library("tidyverse")
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library("glmnet")
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Loaded glmnet 4.1-8
library("caret")
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift

1 Feature Selection

data(Hitters, package = "ISLR")
Hitters <- na.omit(Hitters)

1.1 Forward selection

Implement a forward stepwise selection procedure to find the best five features. The response feature here is Salary which is numeric, so this is a regression and not a classification problem.

selectFeatureDirect <- function(train, test, cls.train, cls.test, features) {
    ## identify a feature to be selected
    current.smallest.mse <- Inf
    selected.i <- NULL
    for (i in i:ncol(train)) {
        current.f <- colnames(train)[i]

        if (current.f %in% c(features, "Salary")) {
            next
        }
        model <- lm(data = cbind(train[, c(features, current.f, "Salary")]), Salary ~
            .)

        test.mse <- mean((cls.test - predict.lm(model, test[, c(features, current.f,
            "Salary")]))^2)

        if (test.mse < current.smallest.mse) {
            current.smallest.mse <- test.mse
            selected.i <- colnames(train)[i]
        }
    }
    return(selected.i)
}
set.seed(5003)
inTrain <- createDataPartition(Hitters$Salary, p = 0.6)[[1]]

train <- Hitters[inTrain, ]
test <- Hitters[-inTrain, ]
salary.train <- Hitters$Salary[inTrain]
salary.test <- Hitters$Salary[-inTrain]

features.direct <- NULL  #store the features

# best 5 feature
for (i in 1:5) {
    selected.i <- selectFeatureDirect(train, test, salary.train, salary.test, features.direct)
    print(selected.i)
    features.direct <- c(features.direct, selected.i)
}
## [1] "CRBI"
## [1] "Hits"
## [1] "PutOuts"
## [1] "Division"
## [1] "Walks"

1.2 Lasso regression

1.2.1 Split data and fit

Split the dataset into 60% train and 40% test. Perform a lasso regression using the glmnet package.

inTrain <- createDataPartition(Hitters$Salary, p = 0.6)[[1]]

train <- Hitters[inTrain, c("CRBI", "Hits", "PutOuts", "Division", "Walks", "Salary")]  #using the best 5 features with the response feature Salary
test <- Hitters[-inTrain, c("CRBI", "Hits", "PutOuts", "Division", "Walks", "Salary")]
salary.train <- Hitters$Salary[inTrain]
salary.test <- Hitters$Salary[-inTrain]

# lasso regression model
x <- model.matrix(Salary ~ ., data = train)[, -1]
y <- salary.train
lasso.model <- glmnet(x, y, alpha = 1)

1.2.2 Inspect coefficients as function of λ

Plot the lasso regression coefficients as a function of λ. Use cross-validation to find the best λ value (You can use the inbuilt CV function in glmnet). Support your answer with a visualization of the selection of λ.

plot(lasso.model, xvar = "lambda", label = TRUE)

cv.lasso <- cv.glmnet(x, y, alpha = 1)
best.lambda <- cv.lasso$lambda.min

cat("Best lambda value:", best.lambda)
## Best lambda value: 8.911837
plot(cv.lasso)