CS-109A Introduction to Data Science

Lab 5: Regularization and Cross-Validation¶

Harvard University
Fall 2018
Instructors: Pavlos Protopapas and Kevin Rader
Lab Instructor: Keivn Rader (today at least)
Authors: Kevin Rader, Rahul Dave, David Sondak, Will Claybaugh, Pavlos Protopapas

Table of Contents¶

0. Learning Goals
1. Review of regularized regression
2. Bootstrapped Sampling Distributions and Confidence Intervals
3. Ridge regression for Simple Regression
4. Ridge regression with polynomial features on a grid
5. Cross-validation --- Multiple Estimates
6. Cross-validation --- Finding the best regularization parameter

Learning Goals¶

In this lab, you will work with some noisy data. You will use simple linear and ridge regressions to fit linear, high-order polynomial features to the dataset. You will attempt to figure out what degree polynomial fits the dataset the best and ultimately use cross validation to determine the best polynomial order. Finally, you will automate the cross validation process using sklearn in order to determine the best regularization paramter for the ridge regression analysis on your dataset.

By the end of this lab, you should:

• Really understand regularized regression principles.
• Have a good grasp of working with ridge regression through the sklearn API
• Understand the effects of the regularization (a.k.a penalization or tuning) parameter on fits from ridge regression
• Understand the ideas behind cross-validation
  • Why is it necessary?
  • Why is it important?
  • Basic implementation details.
• Be able to use sklearn objects to automate the cross validation process.

This lab corresponds to lectures 5, 6, and 7 and maps to homework 4 (and beyond).

Part 1: Review of regularized regression¶

We briefly review the idea of regularization as introduced in lecture. Recall that in the ordinary least squares problem we find the regression coefficients $\boldsymbol{\beta}\in\mathbb{R}^{m}$ that minimize the loss function \begin{align*} L(\boldsymbol{\beta}) = \frac{1}{n} \sum_{i=1}^n \|y_i - \boldsymbol{\beta}^T \mathbf{x}_i\|^2. \end{align*} Recall that we have $n$ observations. Here $y_i$ is the response variable for observation $i$ and $\mathbf{x}_i\in\mathbb{R}^{m}$ is a vector from the predictor matrix corresponding to observation $i$.

The general idea behind regularization is to penalize the loss function to account for possibly very large values of the coefficients $\boldsymbol{\beta}$. Instead of minimizing $L(\boldsymbol{\beta})$, we minimize the regularized loss function \begin{align*} L_{\text{reg}}(\boldsymbol{\beta}) = L(\boldsymbol{\beta}) + \lambda R(\boldsymbol{\beta}) \end{align*} where $R(\boldsymbol{\beta})$ is a penalty function and $\lambda$ is a scalar that weighs the relative importance of this penalty. In this lab we will explore one regularized regression model: ridge regression. In ridge regression, the penalty function is the sum of the squares of the parameters, which is written as \begin{align*} L_{\text{ridge}}(\boldsymbol{\beta}) = \frac{1}{n} \sum_{i=1}^n \|y_i - \boldsymbol{\beta}^T \mathbf{x}_i\|^2 + \lambda \sum_{j=1}^m \beta_{j}^{2}. \end{align*}

In lecture, you also learned about LASSO regression in which the penalty function is the sum of the absolute values of the parameters. This is written as, \begin{align*} L_{\text{LASSO}}(\boldsymbol{\beta}) = \frac{1}{n} \sum_{i=1}^n \|y_i - \boldsymbol{\beta}^T \mathbf{x}_i\|^2 + \lambda \sum_{j=1}^m |\beta_j|. \end{align*}

In this lab, we will show how these optimization problems can be solved with sklearn to determine the model parameters $\boldsymbol{\beta}$. We will also show how to choose $\lambda$ appropriately via cross-validation.

Dataset¶

We will work with a synthetic dataset contained in data/noisypopulation.csv. The data were generated from a specific function $f\left(x\right)$ (the actual mathematical form will remain a mystery). Noise was added to the function to generate synthetic, noisy (aka, real) observations via $y = f\left(x\right) + \epsilon$ where $\epsilon$ was drawn from a random distribution. Even if you could make observations at every single value of $x$, the true function may still be obscured. Of course, the data you actually observe are a subset of all the possible observations in the population. In this lab, we will refer to observations at every single value of $x$ as the population and the subset of observations as in-sample y or simply the observations.

The dataset contains three columns:

1. f is the true function value
2. x is the predictor
3. y is the measured response.

In this lab, we will try out some regression methods to fit the data and see how well our model matches the true function f.

Let's take a quick look at the dataset. We will plot the true function value and the population.

It is often the case that you just can't make observations at every single value of $x$. We will simulate this situation by making a random choice of $60$ points from the full $200$ points. We do it by choosing the indices randomly and then using these indices as a way of getting the appropriate samples.

Note: If you are not familiar with the numpy sort method or the numpy random.choice() method, then please take a moment to look them up in the numpy documentation.

Moving on, let's get the $60$ random samples from our dataset.

Let's take one more look at our data to see which points we've selected.

Now we do our favorite thing and split the sample data into training and testing sets.

Note that here we are actually getting indices instead of the actual training and test set. This is okay and is another way of generating train-test splits.

Great! At this point we've explored our data a little bit, selected a sample of the dataset, and done a train-test split on the sample dataset.

Let's move on to the data analysis. We'll begin with ridge regression. In particular we'll do ridge regression on a single predictor and compare it with simple linear regression.

To start, let's fit the old classic, linear regression.

Part 2: Bootstrapping¶

But wait! Unlike statsmodels, we don't get confidence intervals for the betas. Fortunately, we can bootstrap to build the confidence intervals

1. In the code below, two key steps of bootstrapping are missing. Fill in the code to draw sample indices with replacement and to fit the model to the bootstrap sample. You'll need numpy's np.random.choice. Here's the function documentation in case you need it.
2. Visualize the results, and use numpy's np.percentile: function documentation.

From the above, we find that the bootstrap $90\%$ confidence interval is well away from $0$. We can confidently say that $\beta_{1}$ is not secretly $0$ and we're being fooled by randomness.

Part 3: Ridge regression for Simple Regression¶

To begin, we'll use sklearn to do simple linear regression on the sampled training data. We'll then do ridge regression with the same data, setting the penalty parameter $\lambda$ to zero. What happens when we set $\lambda = 0$ for Ridge's Penalty factor?

We will store the regression coefficients in a dataframe for easy comparison. The cell below provides some code to set up the dataframe ahead of time. Notice that we don't know the actual values in the pandas series, so we just set them to NaN. We will overwrite these later.

We start with simple linear regression to get the ball rolling.

We will use the above $\boldsymbol\beta$ coefficients as a benchmark for comparision to the ridge method. The same coefficients can be obtained with ridge regression, which we demonstrate now.

For reference, here is the ridge regression documentation: sklearn.linear_model.Ridge.

The snippet of code below implements the ridge regression with $\lambda = 0$.

Note: The weight $\lambda$ is referred to as alpha in the documentation.

Remark: $\lambda$ goes by many names including, but not limited to: regularization parameter, penalization parameter, penalty factor, tuning parameter, shrinking parameter, and weight. Regardless of these names, it is a hyperparameter. That is, you set it before you begin the training process. An algorithm can be very sensitive to its hyperparameters and we will discuss how a method for selecting the "correct" hyperparameter values later in this lab.

The beta coefficients for linear and ridge regressions coincide for $\lambda = 0$, as expected. We plot the data and fits.

Make a plot with of the ridge regression predictions with $\lambda = 0, 5, 10, 100$. Be sure to include a legend.

What happens for very large $\lambda$ (e.g. $\lambda \to \infty$)?

Your plot should look something like the following plot (doesn't have to be exact):

Part 3 Recap¶

That was nice, but we were just doing simple linear regression. We really want to do more interesting regression problems like multilinear regression. We will do so in the next section.

Part 4: Ridge regression with polynomial features on a grid¶

Now we'll make a more complex model by adding polynomial features. Instead of building the linear model $y = \beta_0 + \beta_1 x$, we build a polynomial model $y = \beta_0 + \beta_1 x + \beta_2 x^2 + \ldots \beta_d x^d$ for some $d$ to be determined. This regression will be linear though, since we'll be treating $x^2, \ldots, x^d$ themselves as predictors in the linear model.

The design matrix $\mathbf{X}$ contains columns corresponding to $1, x, x^2, \ldots, x^d$. To build it, we use sklearn. (The particular design matrix is also known as the Vandermonde matrix). For example, if we have three observations

\begin{align*} \left\{\left(x_{1}, y_{1}\right), \left(x_{2}, y_{2}\right), \left(x_{3}, y_{3}\right)\right\} \end{align*}
and we want polynomial features up to and including degree $4$, the design matrix looks like

\begin{align*} X = \begin{bmatrix} x_1^0 & x_1^1 & x_1^2 & x_1^3 & x_1^4\\ x_2^0 & x_2^1 & x_2^2 & x_2^3 & x_2^4\\ x_3^0 & x_3^1 & x_3^2 & x_3^3 & x_3^4\\ \end{bmatrix} = \begin{bmatrix} 1& x_1^1 & x_1^2 & x_1^3 & x_1^4\\ 1 & x_2^1 & x_2^2 & x_2^3 & x_2^4\\ 1 & x_3^1 & x_3^2 & x_3^3 & x_3^4\\ \end{bmatrix}. \end{align*}

1. Make a toy vector called toy, where
   \begin{align*} \mathrm{toy} = \begin{bmatrix} 0 \\ 2 \\ 5 \\ \end{bmatrix}. \end{align*}
2. Build the feature matrix up to (and including) degree $4$. Confirm that the entries in the matrix are what you'd expect based on the above discussion.

Note: You may use sklearn to build the matrix using PolynomialFeatures().

We now continue working with our data.

The code provided below is missing a few lines and it's missing many comments. Do the following:

1. Comment every line of the code
   • Normally, you won't do such excessive commenting. In this case, we want to make sure you understand every single line since you didn't actually write this code.
2. Fill in the missing lines
   • Create a ridge regression object at each $\lambda$ value in the list
   • Perform the ridge regression using the fit method from the newly created ridge regression object
   • Make a prediction on the grid and store the results in ypredict_ridge.

Note: We're not giving you an example figure here since we gave you most of the code.

Warning! Make sure you understand the entire code! There are many nice things in there.

As you can see, as we increase $\lambda$ from 0 to 1, we start out overfitting, then doing well, and then our fits develop a mind of their own irrespective of data, as the penalty term dominates.

YOUR DISCUSSION HERE

Part 4 Recap¶

We did a ridge regression on our dataset where the features were the polynomial terms. We also assessed the impact of the regularization parameter on the solution.

Part 5: Cross-validation --- Finding the best penalization parameter¶

In order to determine the critical value of $\lambda$ that gives us our best predictive model, which we'll refer to as $\lambda^*$, we will assess this via cross-validation. To do this we use the concept of a meta-estimator from scikit-learn.

Model selection is supported by two distinct meta-estimators:

1. GridSearchCV
2. RandomizedSearchCV The input to these meta-estimators is an estimator, which has some hyperparameters (e.g. $\lambda$) that need to be optimized, and a set of hyperparameter settings to search through.

The concept of a meta-estimator allows us to wrap, for example, cross-validation, or methods that build and combine simpler models or schemes. For example:

est = Ridge()
    parameters = {"alpha": [1e-8, 1e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 1e-2, 1e-1, 1.0]}
    gridclassifier = GridSearchCV(est, param_grid=parameters, cv=4, scoring="neg_mean_squared_error")

The GridSearchCV replaces the manual iteration over the folds using KFolds and the averaging we just did, doing it all for us. It takes a hyperparameter grid in the shape of a dictionary as input, and sets $\lambda$ to the values you want to try, one by one. It then trains the model using cross-validation, and gets the error for each value of the hyperparameter $\lambda$. Finally it compares the errors for the different $\lambda$'s, and picks the best choice model.

Here is a helper function that we will use to get the best Ridge regression.

We also output the mean cross-validation error at different $\lambda$ (with a negative sign, as scikit-learn likes to maximize negative error which is equivalent to minimizing error).

Now we make a log-log plot of -fit_scores versus fit_lambdas.

SK-learn's cross_val_score: Easier Cross Validation¶

GridSearchCV is an important tool when you are searching over many hyperparameters (and believe us, you will be), but when you only need to get CV scores for a particular model, some students find cross_val_score more intuitive.

We can loop over particular models and get scores for each (equivalent to GridSearchCV over the given parameter settings).

Built-in Cross Validation: RidgeCV and LassoCV¶

Some sklearn models have built-in, automated cross validation to tune their hyper parameters.

Important note:¶

1. For any tool more automated than literally using k_fold, just setting cv=5 will NOT shuffle your data by default.
2. To force shuffling, explicitly pass a KFold object (with shuffling turned on) to the cv argument
3. You may prefer a strategy where you shuffle the rows of your data at the outset of analysis

Part 5b: Refitting on full training set¶

At this point, we have determined the best penalization parameter for the ridge regression on our current dataset using cross validation. Let's refit the estimator on the training set and calculate and plot the test set error and the polynomial coefficients. Notice how many of these coefficients have been pushed to lower values or 0.