Practical 3
Diego Garrido Martín
2020-03-20
1. Introduction
In this practical, you will learn to diagnose the problems of under- and overfitting, as well as to use cross-validation and regularization in both regression and classification problems in real datasets. Try to understand all the code. During the hands-on, you will find a series of short questions (Q) and tasks (T) that you will have to complete. Remember that they can appear as possible questions in the weekly questionaires.
Required libraries
Before starting, install and load the required libraries:
library(ggplot2)
library(caret)
library(glmnet)
library(kernlab)
library(randomForest)
library(xgboost)The caret R package
The caret package (short for Classification And REgression Training) is a set of functions that attempt to streamline the process for creating predictive models. It is probably the most popular package for machine learning in R. The package contains tools for: data splitting, pre-processing, feature selection, model training and evaluation, etc. You can learn more about it here. In this practical you will use it for a variety of tasks.
2. Datasets
Boston housing data
This dataset comes with the MASS R package, and contains data from 506 houses in the suburbs of Boston, including their prices. You can obtain a detailed description of the dataset by typing ?Boston
Diabetes data
This dataset is a corrected version of the original one from the National Institute of Diabetes and Digestive and Kidney Diseases. It contains data from 768 female patients and 9 clinical variables: number of times pregnant (pregnant), plasma glucose concentration (glucose), diastolic blood pressure (pressure), triceps skin fold thickness (thickness), 2-Hour serum insulin (insulin), body mass index (mass), diabetes pedigree function (pedigree), age (age) and diabetes diagnosis (diabetes)
Diabetes multi-class data
This dataset (Reaven and Miller, 1979) presents three different groups of patients (1: overt diabetes,2: chemical diabetes, 3: normal) and 5 clinical variables measured on them: relative weight (relwt), fasting plasma glucose (glufast), area under the plasma glucose curve for the oral glucose tolerance test (glutest), area under the plasma insulin curve for the oral glucose tolerance test (insulin) and steady state plasma glucose response (steady). It contains data from 144 patients.
3. Understanding underfitting and overfitting in regression
Let’s try to exemplify first how an underfitting problem would look like. You will use the Boston dataset in the MASS package to predict the median house value (mdev), based on the predictor variable lstat (percentage of lower status of the population). In this case, instead of employing lm directly, we will use the caret package to perform a 5-fold cross-validation (CV). Here, as we do not have any extra parameter, CV will be used simply to train 5 different models and select the best one according to its RMSE.
# Load the data
data("Boston", package = "MASS")
# Visualize and explore the data before!
# Split the data into training (70%) and test set (30%)
set.seed(123)
training.examples <- createDataPartition(Boston$medv, p = 0.7, list = FALSE)
training.set <- Boston[training.examples, ]
test.set <- Boston[-training.examples, ]
# Define cross-validation schema
cv <- trainControl(## 5-fold CV
method = "cv",
number = 5)
# Train the 'lm' model with caret
model <- train(medv ~ lstat,
data = training.set,
method = "lm",
trControl = cv)
model # caret selects the best model usign CV, and reports the corresponding RMSE## Linear Regression
##
## 356 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 286, 284, 285, 285, 284
## Resampling results:
##
## RMSE Rsquared MAE
## 6.061445 0.5592319 4.456393
##
## Tuning parameter 'intercept' was held constant at a value of TRUE# Visualize the scatter plot of medv vs lstat in the training set and plot lm
ggplot(training.set, aes(x = lstat, y = medv) ) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ x, se = F) +
theme_bw()
T1: Calculate the error (RMSE) in the training set and compare it with the one in the CV set reported above (you do not need to use caret for that, you can use lm as in practical 2). Are they similar? Or very different? Why?
While visualizing the best fit is a possible way to debug your learning algorithm (linear regression in this case), it is not always easy to plot the data and the model, especially when you have many predictors. Learning curves can help you debug a learning algorithm even when it is not easy to visualize the data:
# This function trains 'lm' on increasing sample sizes and computes the training and CV error (it allows to plot a learning curve)
set.seed(123)
LC <- learning_curve_dat(dat = training.set, # Our training set (now using all predictors)
proportion = (1:20)/20, # 20 sets of increasing sample sizes
outcome = "medv", # Our response variable
method = "lm", # Our algorithm (lm)
trControl = cv, # Our CV schema (defined above)
verbose = F) # Do not print progress
## Preprocessing to make the plot (don't worry about it)
LC_train <- LC[LC$Data == "Training", c("RMSE", "Training_Size") ]
LC_cv <- LC[LC$Data == "Resampling", c("RMSE", "Training_Size") ]
LC_cv <- data.frame("Training_Size" = unique(LC_cv$Training_Size),
"RMSE" = tapply(LC_cv$RMSE, LC_cv$Training_Size, median))
##
# Plot the learning curve
ggplot(data = LC_train, aes(x = Training_Size, y = RMSE)) +
geom_smooth(method = "loess", formula = y ~ x, aes(col = "Training"), se = F) +
geom_smooth(method = "loess", formula = y ~ x, data = LC_cv, aes(col = "CV"), se = F) +
coord_cartesian(ylim = c(0, 15)) + # Zoom in
xlab("Training set size (n)") +
theme_bw() 
As you can observe, with a low training set size, the training error is low but the CV error is high. In contrast, when there is a large training set size, both training and CV errors are similar. This reflects a high bias problem in the model: the linear regression model is too simple and thus unable to fit well our dataset. Let’s use polynomial regression to fit a better model for this dataset.
Q1: Why fitting a polynomial regression may help in this scenario? Getting more data could help in this case? Why?
# Train a polynomial regression model (order 10) with caret
set.seed(123)
model <- train(medv ~ poly(lstat, 10), data = training.set,
method = "lm",
trControl = cv)
model # caret selects the best model usign CV, and reports the corresponding RMSE## Linear Regression
##
## 356 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 286, 285, 285, 284, 284
## Resampling results:
##
## RMSE Rsquared MAE
## 6.148648 0.5793997 4.121792
##
## Tuning parameter 'intercept' was held constant at a value of TRUET2: Calculate the error (RMSE) in the training set and compare it with the one in the CV set reported here (you do not need to use caret for that, you can use lm). Are they similar? Or very different? Why?
Note that poly simply builds a hypothesis function of the form:
\[ h(x) = \beta_0 + \beta_1lstat + \beta_2lstat^2 +\beta_3lstat^3 + \ldots +\beta_{10}lstat^{10} \]
You can visualize the fit using geom_smooth:
# Visualize the polynomial regression model
ggplot(training.set, aes(lstat, medv) ) +
geom_point() +
stat_smooth(method = lm, formula = y ~ poly(x, 10), se = F) +
theme_bw()
You should see that the polynomial fit is very complex, and tries to follow the points almost perfectly, likely resulting in a low training error. This is an indicator that the polynomial regression model is overfitting the training data and will not generalize well. To better understand the problem, let’s have a look at the corresponding learning curve:
# This function fits 'lm' on increasing sizes of the training set and the CV set and computes the RMSE (essentially, it allows to plot a learning curve)
set.seed(123)
LC <- learning_curve_dat(dat = data.frame(medv = training.set$medv,
poly(training.set$lstat, 10)), # Our training set (polynomial hypothesis)
proportion = (1:20)/20, # 20 sets of increasing sample sizes
outcome = "medv", # Our response variable
method = "lm", # Our algorithm (lm)
trControl = cv, # Our CV schema (defined above)
verbose = F) # Do not print progress
## Preprocessing to make the plot (don't worry about it)
LC_train <- LC[LC$Data == "Training", c("RMSE", "Training_Size") ]
LC_cv <- LC[LC$Data == "Resampling", c("RMSE", "Training_Size") ]
LC_cv <- data.frame("Training_Size" = unique(LC_cv$Training_Size),
"RMSE" = tapply(LC_cv$RMSE, LC_cv$Training_Size, median))
##
# Plot the learning curve
ggplot(data = LC_train, aes(x = Training_Size, y = RMSE)) +
geom_smooth(method = "loess", formula = y ~ x, aes(col = "Training"), se = F) +
geom_smooth(method = "loess", formula = y ~ x, data = LC_cv, aes(col = "CV"), se = F) +
coord_cartesian(ylim = c(0, 15)) + # Zoom in
xlab("Training set size (n)") +
theme_bw()
As you can see, the training error is low, but the CV error is higher. There is a gap between the training and CV errors, indicating a high variance (overfitting) problem.
Q2: Why our new model is overfitting the data?
To solve our overfitting problem, one option is to use regularization. We can do that with the package glmnet.
Q3: How else could we solve an overfitting problem?
set.seed(123)
model <- train(medv ~ poly(lstat, 10),
data = training.set,
method = "glmnet",
trControl = cv,
)
model # caret selects the best model across CV sets as well as alpha and lambda## glmnet
##
## 356 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 286, 285, 285, 284, 284
## Resampling results across tuning parameters:
##
## alpha lambda RMSE Rsquared MAE
## 0.10 0.01329559 6.039915 0.5864598 4.097869
## 0.10 0.13295595 5.757884 0.6083730 4.036869
## 0.10 1.32955947 5.480063 0.6523224 3.937398
## 0.55 0.01329559 6.007313 0.5891535 4.092816
## 0.55 0.13295595 5.465427 0.6426481 3.943983
## 0.55 1.32955947 5.589746 0.6479949 4.002864
## 1.00 0.01329559 5.954563 0.5929676 4.082567
## 1.00 0.13295595 5.418438 0.6484620 3.923645
## 1.00 1.32955947 5.792818 0.6276331 4.144359
##
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 1 and lambda = 0.1329559.# Visualize the polynomial regression model
ggplot(training.set, aes(lstat, medv) ) +
geom_point() +
geom_line(data = data.frame(medv = predict(model, training.set),
lstat = training.set$lstat),
col = "blue", lwd = 1) +
theme_bw()
# To evaluate the generalization RMSE we use the test set
RMSE(pred = predict(model, newdata = test.set), obs = test.set$medv)## [1] 5.119524Note that you have trained several linear regression models with regularization, using CV to select the best regularization parameters (\(\alpha\) and \(\lambda\)) and the test set to compute the generalization error.
Q4: Explore the model object. Which kind of regularization have you used for the final model (Ridge, Lasso, Elastic net)? Justify your answer. Tip: try to understand the role of \(\alpha\).
Q5: Compare the final polynomial fit with the one obtained before regularization. What do you see? What do you think that would happen if we had a very large \(\lambda\) parameter (e.g. 1000)?
4. Regularized classification using caret
Let’s repeat the classification of the diabetes patients that you did in practical 2, now using caret:
# Load data
diabetes <- read.csv("diabetes.csv", header = TRUE) # read comma-separated file
# Split in training and test
set.seed(123)
training.examples <- createDataPartition(diabetes$diabetes, p = 0.8, list = FALSE)
training.set <- diabetes[training.examples, ]
test.set <- diabetes[-training.examples, ]
# Train model with 5-fold CV ('cv' object defined above)
set.seed(123)
model <- train(diabetes ~ .,
data = training.set,
method = "glm",
family = "binomial",
trControl = cv,
)
# Evaluate model
prediction <- predict(model, newdata = test.set)
confusionMatrix(reference = test.set$diabetes, data = prediction, positive = "pos")## Confusion Matrix and Statistics
##
## Reference
## Prediction neg pos
## neg 44 11
## pos 8 15
##
## Accuracy : 0.7564
## 95% CI : (0.646, 0.8465)
## No Information Rate : 0.6667
## P-Value [Acc > NIR] : 0.05651
##
## Kappa : 0.4356
##
## Mcnemar's Test P-Value : 0.64636
##
## Sensitivity : 0.5769
## Specificity : 0.8462
## Pos Pred Value : 0.6522
## Neg Pred Value : 0.8000
## Prevalence : 0.3333
## Detection Rate : 0.1923
## Detection Prevalence : 0.2949
## Balanced Accuracy : 0.7115
##
## 'Positive' Class : pos
## You can also use glmnet instead of glm to train logistic regression with regularization:
# Train model with 5-fold CV ('cv' object defined above)
set.seed(123)
model <- train(diabetes ~ .,
data = training.set,
method = "glmnet",
family = "binomial",
trControl = cv,
)
# Evaluate model
prediction <- predict(model, newdata = test.set)
confusionMatrix(reference = test.set$diabetes, data = prediction, positive = "pos")## Confusion Matrix and Statistics
##
## Reference
## Prediction neg pos
## neg 44 12
## pos 8 14
##
## Accuracy : 0.7436
## 95% CI : (0.6321, 0.8358)
## No Information Rate : 0.6667
## P-Value [Acc > NIR] : 0.09127
##
## Kappa : 0.4
##
## Mcnemar's Test P-Value : 0.50233
##
## Sensitivity : 0.5385
## Specificity : 0.8462
## Pos Pred Value : 0.6364
## Neg Pred Value : 0.7857
## Prevalence : 0.3333
## Detection Rate : 0.1795
## Detection Prevalence : 0.2821
## Balanced Accuracy : 0.6923
##
## 'Positive' Class : pos
## Q6: Explore the model object. Which kind of regularization are we using (Ridge, Lasso, Elastic net)? Justify your answer. Tip: try to understand the role of \(\alpha\).
Q7 Is regularization improving the classification? Why/Why not?
5. Other algorithms available in caret
Beyond linear and logistic regression, with caret you have plenty of algorithms. For example, in the case of a classification problem:
- Linear Discriminant Analysis
lda - Support Vector Machine with Polynomial Kernel
svmPoly - Random forest
rf - Neural network
nnet - Bosteed trees
xgbTree
T3: Browse the available models in caret and find whether the methods listed above can be used only for classification, only for regression, or for both.
T4: In addition, train each of the algorithms listed above to classify diabetes patients and evaluate them. To do that, simply fill the gap in the code below. Remember to set the same seed before running each method.
# Train model with 5-fold CV
set.seed(123)
model <- train(diabetes ~ .,
data = training.set,
method = " ", # Add here the method e.g. lda
trControl = cv,
)
model
# Evaluate model
prediction <- predict(model, newdata = test.set)
confusionMatrix(reference = test.set$diabetes, data = prediction, positive = "pos")Q8: Which one provides a better performance? Order them from best to worst. Which metric did you use for that? Which algorithm is slower?
Interestingly, some of these methods allow you to perform directly multiclass clasification, such as Linear Discriminant Analysis or Neural Networks. Check that by using lda in our multi-class diabetes dataset:
# Load and inspect the data
diabetes2 <- read.csv("diabetes2.csv", header = TRUE) # read comma-separated file
# Remove the 'id' column (it is just an individual identifier)
diabetes2$id <- NULL
diabetes2$group <- as.factor(diabetes2$group)
# Split the data into training and test set
set.seed(123)
training.examples <- createDataPartition(diabetes2$group, p = 0.7, list = FALSE)
training.set <- diabetes2[training.examples, ]
test.set <- diabetes2[-training.examples, ]
set.seed(123)
model <- train(group ~ .,
data = training.set,
method = "lda",
trControl = cv,
)
# Evaluate model
prediction <- predict(model, newdata = test.set)
confusionMatrix(reference = test.set$group, data = prediction, positive = "pos")## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 2 3
## 1 6 0 0
## 2 3 7 1
## 3 0 3 21
##
## Overall Statistics
##
## Accuracy : 0.8293
## 95% CI : (0.6794, 0.9285)
## No Information Rate : 0.5366
## P-Value [Acc > NIR] : 8.511e-05
##
## Kappa : 0.7098
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 1 Class: 2 Class: 3
## Sensitivity 0.6667 0.7000 0.9545
## Specificity 1.0000 0.8710 0.8421
## Pos Pred Value 1.0000 0.6364 0.8750
## Neg Pred Value 0.9143 0.9000 0.9412
## Prevalence 0.2195 0.2439 0.5366
## Detection Rate 0.1463 0.1707 0.5122
## Detection Prevalence 0.1463 0.2683 0.5854
## Balanced Accuracy 0.8333 0.7855 0.8983Q9: Is Linear Discriminant Analysis giving a better result than your multi-class classification attempt in practical 2?
Remember to answer the short questions (Q) and complete the tasks (T). They can appear as possible questions in the weekly questionaires.