Loading required R packages
tidyverse for easy data manipulation and visualization
corrplot for correlations among the variables visualization
caret for easy machine learning workflow
glmnet for computing penalized regression
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(corrplot))
suppressPackageStartupMessages(library(caret))
suppressPackageStartupMessages(library(glmnet))
Preparing the data
# Load the data
data("Boston", package = "MASS")
str(Boston)
## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
# data visualization
p.cor<-cor(Boston)
corrplot.mixed(p.cor)
We’ll randomly split the data into training set (80% for building a predictive model) and test set (20% for evaluating the model). Make sure to set seed for reproducibility.
# Split the data into training and test set
set.seed(123)
training.samples <- Boston$medv %>% #medv is the response variable
createDataPartition(p = 0.8, list = FALSE)
# createDataPartition() from "caret" package, split data proportionally
train.data <- Boston[training.samples, ]
test.data <- Boston[-training.samples, ]
input matrix x: This should be created using the function model.matrix() allowing to automatically transform any qualitative variables (if any) into dummy variables, which is important because glmnet() can only take numerical, quantitative inputs. After creating the model matrix, we remove the intercept component at index = 1.
# Predictor variables
x <- model.matrix(medv~., train.data)[,-1]
# Response variable
y <- train.data$medv
cv.glmnet function
cv.glmnet(x, y, family, lambda, type.measure, nfolds, foldid, grouped, keep, parallel, alpha, ...)
x: input matrix
y: response variable
family: 分布函数。因变量为连续分布(gaussian),因变量为二项分布(binomial)
type.measure: mse(gaussian models) auc(two-class logistic regression), class(binomial and multinomial logistic regression)
nfolds: number of folds - default is 10
k-fold cross-validated ridge regression
# Find the best lambda using cross-validation
set.seed(123)
cv <- cv.glmnet(x, y, alpha = 0)
# Display the best lambda value (when MSE lowest)
cv$lambda.min
## [1] 0.7581162
# Fit the final model on the training data
ridge <- glmnet(x, y, alpha = 0, lambda = cv$lambda.min)
# Display regression coefficients
coef(ridge)
## 14 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 28.696325628
## crim -0.072846256
## zn 0.034166905
## indus -0.057447836
## chas 2.491226902
## nox -11.092319943
## rm 3.981320544
## age -0.003137072
## dis -1.192961128
## rad 0.140678482
## tax -0.006099569
## ptratio -0.863997247
## black 0.009365355
## lstat -0.479141586
# Make predictions on the test data
x.test <- model.matrix(medv ~., test.data)[,-1]
predictions <- ridge %>% predict(x.test) %>% as.vector()
# Model performance metrics
data.frame(
RMSE = RMSE(predictions, test.data$medv),
Rsquare = R2(predictions, test.data$medv)
)
## RMSE Rsquare
## 1 4.983643 0.671074
k-fold cross-validated LASSO regression
# Find the best lambda using cross-validation
set.seed(123)
cv <- cv.glmnet(x, y, alpha = 1)
# Display the best lambda value (when MSE lowest)
cv$lambda.min
## [1] 0.008516101
# Fit the final model on the training data
lasso <- glmnet(x, y, alpha = 1, lambda = cv$lambda.min)
# Display regression coefficients
coef(lasso)
## 14 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 36.905385212
## crim -0.092224281
## zn 0.048424768
## indus -0.008413932
## chas 2.286235925
## nox -16.796507506
## rm 3.811860949
## age .
## dis -1.596031128
## rad 0.285461530
## tax -0.012401417
## ptratio -0.950408780
## black 0.009646870
## lstat -0.528799085
# Make predictions on the test data
x.test <- model.matrix(medv ~., test.data)[,-1]
predictions <- lasso %>% predict(x.test) %>% as.vector()
# Model performance metrics
data.frame(
RMSE = RMSE(predictions, test.data$medv),
Rsquare = R2(predictions, test.data$medv)
)
## RMSE Rsquare
## 1 4.989791 0.6707754
k-fold cross-validated elastic net regression
The elastic net regression can be easily computed using the caret workflow, which invokes the glmnet package.
We use caret to automatically select the best tuning parameters alpha and lambda. The caret packages tests a range of possible alpha and lambda values, then selects the best values for lambda and alpha, resulting to a final model that is an elastic net model.
Here, we’ll test the combination of 10 different values for alpha and lambda. This is specified using the option tuneLength.
The best alpha and lambda values are those values that minimize the cross-validation error
For cross-validated ridge and lasso model, similar codes are used.
# Build the model using the training set
set.seed(123)
elastic <- train(
medv ~., data = train.data, method = "glmnet",
trControl = trainControl("cv", number = 10),
tuneLength = 10
)
# Best tuning parameter
elastic$bestTune
## alpha lambda
## 93 1 0.0170322
# Coefficient of the final model. You need to specify the best lambda
coef(elastic$finalModel, elastic$bestTune$lambda)
## 14 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 36.256402548
## crim -0.088577749
## zn 0.046967102
## indus -0.008012063
## chas 2.269725178
## nox -16.361460313
## rm 3.829573241
## age .
## dis -1.561700859
## rad 0.271508808
## tax -0.011838778
## ptratio -0.944725503
## black 0.009594568
## lstat -0.529927972
# Make predictions on the test data
x.test <- model.matrix(medv ~., test.data)[,-1]
predictions <- elastic %>% predict(x.test)
# Model performance metrics
data.frame(
RMSE = RMSE(predictions, test.data$medv),
Rsquare = R2(predictions, test.data$medv)
)
## RMSE Rsquare
## 1 4.986551 0.6710845
repeated k-fold cross-validated elastic net regression
The process of splitting the data into k-folds can be repeated a number of times, this is called repeated k-fold cross validation.
The final model error is taken as the mean error from the number of repeats.
We generally recommend the (repeated) k-fold cross-validation to estimate the prediction error rate. It can be used in regression and classification settings.
Another alternative to cross-validation is the bootstrap resampling methods, which consists of repeatedly and randomly selecting a sample of n observations from the original data set, and to evaluate the model performance on each copy.
# Build the model using the training set
set.seed(123)
elastic <- train(
medv ~., data = train.data, method = "glmnet",
trControl = trainControl(method = "repeatedcv", number = 10, repeats = 3), #10-fold cross validation with 3 repeats
tuneLength = 10
)
# Best tuning parameter
elastic$bestTune
## alpha lambda
## 54 0.6 0.0393466
# Coefficient of the final model. You need to specify the best lambda
coef(elastic$finalModel, elastic$bestTune$lambda)
## 14 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 35.513125388
## crim -0.085096843
## zn 0.045515676
## indus -0.011769111
## chas 2.270041356
## nox -15.920104479
## rm 3.850861006
## age .
## dis -1.529248662
## rad 0.253320815
## tax -0.011020254
## ptratio -0.936756271
## black 0.009520818
## lstat -0.529225087
# Make predictions on the test data
x.test <- model.matrix(medv ~., test.data)[,-1]
predictions <- elastic %>% predict(x.test)
# Model performance metrics
data.frame(
RMSE = RMSE(predictions, test.data$medv),
Rsquare = R2(predictions, test.data$medv)
) #RMSE lower
## RMSE Rsquare
## 1 4.984612 0.6712767
References
Penalized Regression Essentials: Ridge, Lasso & Elastic Net
Cross-Validation Essentials in R