Goal

Before the break, we started talking about ridge regression. Today, we are going to motivate when and why we might want to consider using this technique. We will also see how to use R to train these models.

The Data

We are going to work with a data set from your textbook called the Hitters data set. This data set provides information about major league baseball players with the goal of predicting \(Y\) = the salary of the players in thousands of U.S. dollars.

To load the data, run the following code:

# Load the library 
library(ISLR)

# Load the data 
data(Hitters)

# Remove columns we don't need for today
Hitters <-Hitters[,-c(14,15,19,20)]

# Scale the matrix (we'll talk later in the lab about why!)
Hitters <- as.data.frame(scale(Hitters))

Our data set has \(n = 322\) rows with 16 possible features. All of these features are numeric. Note: Categorical features are allowed in Ridge Regression! We just aren’t using them for our example today.

If you open the data and look at the features, you will notice that they all seem to contain relatively small numbers. This is because each feature has been scaled to convert it into a z-score. This means a value of 1 for Hits means that for this player, their number of total hits this season is 1 standard deviation above the average number of hits across all players. A value of -2 means that for this player, their number of total hits this season is 2 standard deviations below the average number of hits across all players.

This may seem like an odd scale to use for features. However, for shrinkage techniques like ridge regression, this is very useful. Our technique will involve using the numeric value of the coefficients as a penalty, so it is important that all our features are on the same scale.

# Set a random seed 
set.seed(100)

# Create the betas 
betas <- rnorm(ncol(Hitters), mean =5, sd = 4)

# Add on an intercept
betas <-c(150, betas)

# Create Y
Y <- as.matrix(cbind(1, Hitters))%*%(as.matrix(betas)) +
      rnorm(nrow(Hitters),0,25)

LSLR: Checking for Independence

Models like LSLR assume that we do not have strong relationships among the features that we will use to predict \(Y\). The model treats each feature as though it were independent. This means that before we use LSLR, it is important to check to see if is reasonable to assume that the features are independent.

One way to check this is to create a correlogram, which is basically a plot checking the correlation between each feature and every other feature. This isn’t a perfect tool, as it only checks for linear relationships among numeric features, but is a useful one.

To create a correlogram in R, we use the following.

# Load the library 
library(corrplot)

# Create the plot 
corrplot(cor(Hitters), method = "circle",type ="upper")

The plot you will see is an upper triangular matrix. I specified that in the code using type = "upper". If you prefer a lower triangular matrix, you can use type = "lower".

Each box in the matrix represents the correlation between two features. For instance, the 2nd box on the 1st row represents the correlation between AtBat and Hits. The correlation is represented with a dark blue circle to indicate a strong, positive correlation. The larger the circle and the darker the color, the higher the correlation between two features is. This means that the correlation between AtBat and Hits is close to 1.

You will notice a line of solid blue circles along the diagonal of the plot, as every feature is perfectly correlated with itself.

This plot has a lot of different options if you don’t like the circles. See what happens if you change method to method = "number", method = "ellipse", or method = "pie". This website gives a nice overview of how much you can personalize this plot: http://www.sthda.com/english/wiki/visualize-correlation-matrix-using-correlogram

Whichever version of the plot you prefer, you should notice a complicated structure of relationships among the features. A lot of features are correlated with other features, which means that we are violating an assumption of linear regression. What does this mean for our ability to estimate the relationship between the features and a response? Let’s explore.

LSLR: Training the Model

Because we created the response variable ourselves, we actually know the true value of \(\boldsymbol{\beta}\). This means we can compare our estimated coefficients from LSLR to the true values of the parameters.

Estimates$Difference <- ( Estimates$LSLR - Estimates$Truth)  

Right now, our LSLR model is not able to estimate some of the coefficients well. We were expecting that because some of our conditions were violated, but now we actually see the problem in action.

The trick now is to figure out how we might estimate the coefficients differently in order to improve our results.

Ridge Regression

In LSLR, we choose the value of \(\boldsymbol{\hat{\beta}}\) that gives us the smallest value possible for the RSS.

\[RSS = \sum_{i=1}^n (y_i - \hat{y}_i)^2\]

However, when our features are correlated, we see that choosing the value of \(\boldsymbol{\hat{\beta}}\) that makes the RSS its absolute smallest yields estimates that do not capture the true value of \(\boldsymbol{\beta}\) very well.

There are many, many different possible values of \(\boldsymbol{\hat{\beta}}\). Each one yields a potentially different value of the RSS if we choose it. The value of \(\boldsymbol{\hat{\beta}}\) that we got from LSLR yields the smallest possible value of the RSS, but the estimates were not reflective of the true parameters.

Specifically, what we noticed is that many of our coefficient estimates from LSLR are too large in absolute value. For example, the coefficient for CHits (the cumulative number of hits over the season) is supposed to be 1.70, but with LSLR we obtain an estimate of -46.68, a value over 20 times larger in absolute value than the truth!

Because we are noticing that the coefficients are too large, we are going to change how we estimate \(\boldsymbol{\hat{\beta}}\). Instead of choosing the value of \(\boldsymbol{\hat{\beta}}\) that minimizes the RSS, we will choose the value that minimizes:

\[RSS + \lambda \sum_{j} \hat{\beta}_j^2,\]

This means we are now choosing estimates of \(\boldsymbol{\hat{\beta}}\) based on two components:

\[\underbrace{RSS}_{\text{Data Component}} + \underbrace{\lambda \sum_{j} \hat{\beta}_j^2}_{\text{Penalty Component}}\] The penalty component is called a shrinkage penalty. This is because the penalty term gets larger as the square of the coefficients gets larger.

Because we want to minimize this entire term (the data component + the penalty component), we can encourage the choice of coefficients that have smaller values than we get without the penalty component. Since the problem we have currently with LSLR is that our coefficients are too big, this is a reasonable choice! Estimating the coefficients in this way is called using ridge regression.

Think about it like a budget. Our goal is to make the RSS small. In LSLR, our budget is unlimited - we can choose any coefficients we want to make the RSS as small as possible. However, in Ridge, we have a budget. We still want to make the RSS small, but we are restricted in how large we can let our coefficients be in order to achieve that goal.

The key to ridge regression is choosing the tuning parameter \(\lambda\), which is a scalar where \(\lambda \geq 0\). This tuning parameter allows us to control how much shrinkage we want on the coefficient estimates. In other words, it allows us to control the balance between the data component and the penalty component; this controls the budget. A larger \(\lambda\) means a smaller budget.

library(glmnet)

# Train Ridge
ridge.model <- glmnet(Hitters, Y, alpha = 0, lambda = )

# Store the coefficients 
coefficients(ridge.model)

sum(coefficients( model )[-1]^2)

Ridge: Choosing The Tuning Parameter

We want to choose a value of \(\lambda\) that shrinks the coefficients, but we do not want to shrink them so much that the relationship between the features and \(Y\) is no longer represented by those coefficients. To help us choose an appropriate value of \(\lambda\), we are going to return to a concept that is very familiar to us now - cross validation.

The process begins as we did when choosing the tuning parameter \(K\) in KNN - we will try a whole bunch of different values of the tuning parameter \(\lambda\). For each value of \(\lambda\), we use cross validation (specifically 10-fold CV today) to estimate the test RMSE we would get with that choice of \(\lambda\).

If we choose a value of \(\lambda\) that starts to break the relationship between the features and the response, the test RSME will increase. This means the test RMSE is a good tool to us to help us select \(\lambda\). We choose the value of \(\lambda\) that gives us the smallest validation RMSE.

In summary, we start with a list of possible choices for \(\lambda\). For each \(\lambda\) on the list, we

• 1. Solve for the ridge estimator \(\boldsymbol{\hat{\beta}}_{Ridge}\).
• 2. Use 10-fold CV to estimate the test RMSE with that choice of \(\boldsymbol{\hat{\beta}}_{Ridge}\).

We then ultimately select the value of \(\lambda\) for which the validation RMSE is the smallest.

# Set a random seed 
set.seed(100)

# Run 10 fold-CV for a variety of choices of lambda
cv.ridge.mod <- cv.glmnet(as.matrix(Hitters), Y , alpha = 0 , lambda = )

seq( from = , to = , by = )

ridgePlot <- function(ridge.mod, metric, title){
  suppressMessages(library(ggplot2))
  
  smallestLambda <- ridge.mod$lambda[which.min(ridge.mod$cvm)] 
  
  if(metric == "MSE"){
  g1 <- ggplot( data.frame(ridge.mod$lambda), aes( x = ridge.mod$lambda, y = (ridge.mod$cvm))) + geom_point() + geom_vline( xintercept = smallestLambda, col = "blue" , lty = 2 ) +  labs(caption = paste("Validation MSE values for Different Tuning Parameters. Smallest MSE at lambda = ", smallestLambda), title = title, y = "Validation MSE", x = "Tuning Parameter")
  
  }
  
  if(metric == "RMSE"){
  g1 <- ggplot( data.frame(ridge.mod$lambda), aes( x = ridge.mod$lambda, y = sqrt(ridge.mod$cvm))) + geom_point() + geom_vline( xintercept = smallestLambda, col = "blue" , lty = 2 ) +  labs(caption = paste("Validation RMSE values for Different Tuning Parameters. Smallest RMSE at lambda = ", smallestLambda), title = title, y = "Validation RMSE", x = "Tuning Parameter")

  }
  
  g1
}

ridgePlot(cv.ridge.mod, metric = , title = )

cv.ridge.mod$lambda.min

min(sqrt(cv.ridge.mod$cvm))

sqrt(cv.ridge.mod$cvm[which(cv.ridge.mod$lambda==0)])

# Train Ridge
ridge.final <- glmnet(Hitters, Y, alpha = 0, lambda = )

Estimates$Ridge <- matrix(coefficients(ridge.final))[,1]  

Estimates <- data.frame("Truth" = betas)
Estimates$LSLR <- m1$coefficients
Estimates$Ridge <- matrix(coefficients(ridge.final))[,1]
Estimates$LSLRvsTruth <- Estimates$LSLR - Estimates$Truth
Estimates$RidgevsTruth <- Estimates$Ridge - Estimates$Truth
Estimates$CloserToTruth <- ifelse( abs(Estimates$LSLRvsTruth) < abs(Estimates$RidgevsTruth), "LSLR", "Ridge")