Implementing ridge regression in R

Ridge regression is one of the most popular regularization methods in machine learning. Regularization refers to a set of techniques that reduce the impact of overfitting during model training. Ridge regression, also known as L2 regularization, is one of several types of regularization techniques for linear regression models.

Another well-known regularization technique is Lasso regression or L1 regularization. Regularization techniques are primarily meant to create models with lower variance at the cost of introducing a small amount of bias. This means that models will be able to better withstand smallchanges in the underlying properties of inputted data at the time of prediction relative to when the model was trained. This can help to enhance model accuracy, and address the consequences of multicollinearity in the models that they produce. Note that Ridge regression may also be applied in logistic regression for binary classification problems.

Overfitting and the bias-variance tradeoff

Overfitting refers to the tendency of machine learning models to learn “noise” in the data that they are trained on, which means that the model will start to rely on irrelevant information when computing an output. While it can be minimized, noise comes from imperfections in data collection processes, missing data and outlier values, as well as other sources, and is ultimately impossible to completely remove.

During model training, machine learning algorithms seek to find the correct parameter values to optimize an objective function by using training data to learn the relationships between different input variables. During this process, models can learn relationships that exist purely due to noise, and not because of valid relationships, which therefore results in a model that is “overfit” to its training data. This can lead to models that appear to have high accuracy during model training, but that produce poor results when deployed, since they will have different data inputted into them with different noise than what was used during training.

The inability to completely remove noise from data also introduces one of the principal problems that regularization aims to address: the trade-off between bias and variance in supervised machine learning models. Bias refers to the difference between a model’s predicted values, and what that value should actually be. Variance refers to the model’s ability to produce consistent outputs across different input datasets which may contain different levels of noise.

Minimizing bias during model training maximizes overall accuracy, however this comes at the cost of overfit models which are not generalizable to other data. However, minimizing variance is done by making models simpler, so that they do not rely so much on relationships that exist due to noise when making predictions. This results in the model not being able to fully capture the relationships between different variables and therefore achieve low overall accuracy again. While the specifics on how this is achieved depend on the individual technique, the overall goal of regularization is to find the optimal balance between bias and variance so as to produce models that are accurate while simultaneously being generalizable.

Key characteristics of ridge regression

Ridge regression aims to introduce a small amount of bias to achieve a significant reduction in variance. This is done by “shrinking” the values of regression coefficient estimates towards zero in a process called shrinkage. Shrinkage is achieved by introducing a regularization penalty parameter, lambda (λ), into the model’s loss function that minimizes the impact of each feature, reducing the contribution of noise to coefficient estimates. Note that coefficient values cannot be reduced completely to zero, unlike in lasso regression. Ordinary least squares (OLS) regression seeks to minimize the sum of the squared distances between the datapoints and a fitted hyperplane in a high dimensional space. Ridge modifies this process by minimizing the sum of the squared distances plus lambda times the sum of squares of the coefficients, thus introducing a small amount of bias:

The λ penalty added to the residual sum of squares (RSS) is proportional to the square of the magnitude of the coefficients. As λ increases, the shrinkage increases, and the coefficients are pushed closer to zero. In R’s glmnet package, the parameters alpha (α) and lambda (λ) are available to choose the regularization method and shrinkage penalty efficiently. Setting alpha (α) equal to 0 corresponds to implementing Ridge regression.

Ridge regression is instrumental in dealing with multicollinearity, where predictor variables are highly correlated. Although the ordinary least squares (OLS) estimates are unbiased in multicollinearity, their large variances can result in a substantial deviation of the observed values from the true values. Ridge regression enhances the reliability of predictions by introducing a degree of bias to the regression estimates, thereby reducing standard errors.

Ridge regression versus alternative regularization methods

Because Ridge regression does not allow coefficients to shrink completely to zero, it is useful to use when you are dealing with multicollinearity and when removing correlated variables is undesirable or impossible. However, Lasso regression makes coefficients equal to zero and is able to completely remove the effects of certain variables on the model. This makes Lasso ideal for feature selection, or the removal of variables, and for model interpretability. Lasso might also introduce significant bias for large coefficients by overly penalizing them, leading to too much underfitting.

Then there is Elastic net regression, which is another technique used to improve the precision and interpretability of statistical models. Elastic net regression serves as a middle ground between ridge and lasso regression by combining the L1 and L2 penalties.

Implementing ridge regression

This tutorial is a beginner’s guide to implementing ridge regression in R. You will pretend that you are a Data Scientist aiming to accurately predict the energy efficiency of a building.

This data set:

• Objective: The goal of this data set is to predict the energy efficiency of a building, which is quantified by a specific metric (for example, heating load and cooling load). This metric is influenced by various architectural and material features of the building.
• Features: The data set consists of 100 features (n_features = 100) that represent different aspects of the building's design and surrounding environment. These could include wall materials, insulation thickness, window type, orientation, surrounding vegetation, and local climate conditions, among others.
• Samples: There are 1,000 buildings (n_samples = 1000) in the data set, each with measurements or specifications for the 100 features.
• Noise: The data set includes a level of noise (noise_level = 5) to account for measurement errors or unobserved variables that might affect the energy efficiency metric.
• Sparsity and Relevance: Out of the 100 features, only 10 are truly relevant to predicting energy efficiency, mimicking the real-world scenario where only a subset of features significantly impacts the outcome. These relevant features have varying degrees of influence, represented by the true coefficients generated with a normal distribution.
• Multicollinearity: The features exhibit multicollinearity, where some features are highly correlated with each other. This is common in real-world data, as certain design aspects of buildings might be related (for example, the use of specific materials might correlate with insulation properties).

We’ll cover all the steps in implementing ridge regression in R, starting with exploratory data analysis, scaling the data, splitting the data set, fitting the ridge regression model, determining the optimal λ value, fitting the final model with the optimal λ value, and finally making predictions and evaluating the model.

Let’s begin!

Prerequisites

Create an IBM Cloud account so that you can create a watsonx.ai project.

A Jupyter notebook containing the R code samples used in this tutorial are available in this GitHub Repository for you to download and import into your watsonx project.

Steps

Step 1: Set up your environment

While you can choose from several tools, this tutorial walks you through how to set up an IBM account to use a Jupyter Notebook. Jupyter Notebooks are widely used within data science to combine code, text, images, and data visualizations to formulate a well-formed analysis.

1. Log in to watsonx.ai using your IBM Cloud account.

2. Create a watsonx.ai project.

3. Create a Jupyter Notebook and select an R runtime.

This step will open a Notebook environment where you can copy the code from this tutorial to tackle a simple regression problem. The code is available in the IBM Developer watsonx GitHub repository.

Step 2: Install and load the necessary packages

In this step, we will create our first code cell and install and load the necessary R packages:

• dplyr: A data manipulation package that simplifies and enhances the process of manipulating and summarizing data frames.
• ggplot2 : A data visualization package that enables users to create a wide range of static and dynamic plots.
• glmnet : An R package for fitting generalized linear models with regularization using elastic net regularization techniques.
• tidyr : A package designed to help you create “tidy” data.

#install.packages(c('glmnet', 'ggplot2', 'tidyr', ‘dplyr’))

library(glmnet)
library(ggplot2)
library(tidyr)
library(dplyr)

Note that the install.packages command is commented, assuming that the required libraries are already installed. If you encounter an error, uncomment this line and re-run the code cell.

Step 3: Generate the data set

In this step, we will generate a synthetic data set. More detailed explanations of each step can be found in the accompanying jupyter notebook.

set.seed(42) # Ensure reproducibility

n_features <- 100 #
n_samples <- 1000 # Number of samples
noise_level <- 5 # Noise level

# Generate true coefficients with sparsity and relevance
true_coef <- rep(0, n_features)
# Designate 10 features as relevant
true_coef[sample(1:n_features, 10)] <- rnorm(10, mean = 10, sd = 5)


# Generate design matrix X with correlated features
X <- matrix(rnorm(n_samples * n_features), nrow = n_samples, ncol = n_features)
for (i in 1:(n_features - 1)) {
  X[, i + 1] <- 0.5*X[, i] + rnorm(n_samples) # Introduce multicollinearity
}

# Generate target variable y with linear relationship and noise
y <- X %*% true_coef + rnorm(n_samples, mean = 0, sd = noise_level)

The resulting matrix can be turned into a data frame with R’s data.frame() function:

df <- data.frame(X)
df$Y <- y

Now that we have our data inside a data.frame object, it will be easy to manipulate it and to fit Ridge and other machine learning models. Let’s visualize the first few rows using the head() function:

head(df)

Step 4: Exploratory data analysis (EDA)

In this step, we will pretend that we don’t know how our data was generated already, and explore it how one would were this a real data set.

During EDA, you normally investigate multiple things such as:

• Distributions of your outcome variable and features
• Corrleations between features
• Whether your dataset is imbalanced for classification tasks
• Whether there are missing values
• Identification of outliers
• The meaning of your features and aquire more domain knowledge if necessary

First, let’s create a scatter plot to visualize the relationship between some of our features and our outcome variable. Scatterplots can help identify correlations between features and identify outliers:

# Select a subset of predictors for visualization
selected_predictors <- c("X1", "X4", "X9")

# Reshape the data to long format for selected predictors
data_long <- pivot_longer(df, cols = selected_predictors, names_to = "Predictor", values_to = "Value")

# Create the plot
ggplot(data_long, aes(x = Value, y = Y, color = Predictor)) +
  geom_point() +
  theme_minimal() +
  labs(title = "Relationship between Selected Predictors and Y",
       x = "Predictor Value",
       y = "Outcome Y") +
  scale_color_manual(values = c("X1" = "blue", "X4" = "red", "X9" = "green"))

Next, let’s create a histogram to visualize the distribution of one of our features. This can be done quickly with R’s hist() function. As you can see, this feature is clearly normally distributed:

hist(df$X40)

Finally, let's make a correlation matrix of some of our features. This is another useful technique for quickly identifying collinearity in your features. It will be hard to visualize 100 features, so we will only show 6 for this example:

Step 5: Scale the data

Feature scaling is crucial in regression to ensure a balanced contribution of features to the regularization term. Ridge regression attempts to minimize the sum of squared residuals between the observed data and a fitted hyperplane in a high-dimensional space, while also penalizing the sum of squared coefficients by a regularization term. This regularization term penalizes large coefficients, which helps prevent overfitting, especially when features are on different scales. Without proper scaling, features with larger scales can exert an outsized influence on the model.

Feature scaling is the process of standardizing scales across features. Standarization is one of the most common methods of feature scaling, and is done by substracting the mean from an observation and then dividing by its standard deviation. We do not need to perform this manually as the glmnet function will perform it for us when we fit out model.

Step 6: Split the data into training and test data sets

Now, we will split our data into training and test data sets. This is crucial to avoid overfitting during model training, and can help evaluate the performance of the model on unseen data, making sure that it is learning real information and not just relationships that appear because of noise.

set.seed(123)  # For reproducibility
train_index <- sample(1:n_samples, size = round(0.8 * n_samples))
test_index <- setdiff(1:n_samples, train_index)
X_train <- X[train_index, ]
y_train <- y[train_index]
X_test <- X[test_index, ]
y_test <- y[test_index]

Step 7: Determine the pptimal lambda value

Ordinary least squares (OLS) regression attempts to minimize the sum of the squared residuals between the training data and the best fit line. N is the number of observations, and i is the index of a given observation.

$RSS(\beta) = \sum_{i=1}^{n} (y_i - (X_i \beta))^2$

With Ridge, we will introduce a regularization penalty $\lambda$ into this equation. P refers to the number of parameters in the model.

$\quad RSS_{\text{Ridge}}(\beta) = \sum_{i=1}^{n} (y_i - (X_i \beta))^2 + \lambda \sum_{j=1}^{p} \beta_j^2$

There are many techniques to find the best value for lambda. One of the most common techniques is known as k-fold cross validation. This method splits your testing data into k parts, or folds, including 1 fold to be used as a testing fold for evaluation of model performance within the training set.

For a predefined set of lambda values, the model will be trained k-1 times and then tested on the testing fold where a metric recording indicating the model's performance will be recorded ("mean squared error" for regression). Afterwards, a different fold will be set aside as the test fold, and this process will repeat until each fold has been set aside for testing. After each fold is tested, an average of all of the folds metrics is computed, which is the final cross validation score for a given lambda. The next lambda in the set is tested, and the value of lambda with the best CV score is selected as the optimal value for the final model.

This approach ensures that the lambda picked results in a generalizable model, since it is tested on unseen data during hyperparameter selection before testing.

We'll perform cross validation with the cv.glmnet() function over a pre-specified range of 1000 values. When you don't specify a range of lambdas to test, glmnet will automatically define one for you.

lambda_values <- 10^seq(-3, 3, length = 1000)

# Perform cross-validation to select the best lambda
cv_ridge_model <- cv.glmnet(X_train, y_train, alpha = 0, lambda = lambda_values, 
                            nfolds = 10, standardize = TRUE)

# Best lambda value
best_lambda <- cv_ridge_model$lambda.min

Step 8: Fit the optimal Ridge model

Now, we’ll fit a standard OLS model with no regularization pentalty and compare it’s performance to our Ridge model. We’ll do this by calculating the root mean squared error (RMSE) on our test set using both model. This metric is calculated by squaring the difference of the model’s predicted values from the true value, and then taking it’s square root. We square the errors because we are only interested in the distance between the predicted and actual values, and we don’t want negative errors to “cancel out” the impact of positive ones. Taking the square root of this metric ensures that we are using the same units as our target variable.

A lower RMSE means that the model’s predictions had less overall error:

ridge_model <- glmnet(X_train, y_train, alpha = 0, lambda = best_lambda, standardize = TRUE)
# Fit OLS model
ols_model <- lm(y_train ~ X_train)

# Evaluate performance of ridge and OLS on test set
ols_pred <- predict(ols_model, data.frame(X_test))
ridge_pred <- predict(ridge_model, X_test)
ols_rmse <- sqrt(mean((y_test - ols_pred)^2))
ridge_rmse <- sqrt(mean((y_test - ridge_pred)^2))

Next, we’ll create a bar graph comparing the RMSE of our Ridge and OLS models. Our Ridge model clearly is superior for this dataset:

rmse_data <- data.frame(Model = c("OLS", "Ridge"), RMSE = c(ols_rmse, ridge_rmse))

ggplot(rmse_data, aes(x = Model, y = RMSE, fill = Model)) +
  geom_bar(stat = "identity", color = "black", fill = c("blue", "green"), width = 0.5) +
  geom_text(aes(label = RMSE), vjust = -0.3, color = "black", size = 3.5, position = position_dodge(width = 0.9)) +
  labs(title = "Comparison of RMSE between OLS and Ridge Models",
       x = "Model", y = "RMSE") +
  theme_minimal()

Summary

Ridge Regression can help to create robust and generalizable models by introducing a small amount of bias during model training. It is a technique grounded in the bias and variance tradeoff, and it can help to deal with noise and multicollinearity in training data sets without having to completely remove features.

Using R in watsonx.ai, it is simple to import data, explore and prepare it, and to train Ridge and other machine learning models with the glmnet package.

Try watsonx for free

Build an artificial intelligence (AI) strategy for your business on one collaborative AI and data platform called IBM watsonx, which combines new generative AI capabilities powered by foundation models and traditional machine learning into a powerful platform spanning the AI lifecycle. With watsonx.ai, you can train, validate, tune, and deploy models with ease and build AI applications in a fraction of the time with a fraction of the data.

Create a free IBM Cloud accout to try watsonx for free.

Next steps

Explore more articles and tutorials about watsonx on IBM Developer.

Also, to learn more about other supervised machine learning models that you can apply to classification and regression problems, see these tutorials in the Getting started with machine learning learning path.