Key Word(s): model selection, cross validation, ridge, lasso, regularization
CS109A Introduction to Data Science
Standard Section 4: Regularization and Model Selection¶
Harvard University
Fall 2019
Instructors: Pavlos Protopapas, Kevin Rader, and Chris Tanner
Section Leaders: Marios Mattheakis, Abhimanyu (Abhi) Vasishth, Robbert (Rob) Struyven
#RUN THIS CELL
import requests
from IPython.core.display import HTML
styles = requests.get("http://raw.githubusercontent.com/Harvard-IACS/2018-CS109A/master/content/styles/cs109.css").text
HTML(styles)
For this section, our goal is to get you familiarized with Regularization in Multiple Linear Regression and to start thinking about Model and Hyper-Parameter Selection.
Specifically, we will:
- Load in the King County House Price Dataset
- Perform some basic EDA
- Split the data up into a training, validation, and test set (we'll see why we need a validation set)
- Scale the variables (by standardizing them) and seeing why we need to do this
- Make our multiple & polynomial regression models (like we did in the previous section)
- Learn what regularization is and how it can help
- Understand ridge and lasso regression
- Get an introduction to cross-validation using RidgeCV and LassoCV
# Data and Stats packages
import numpy as np
import pandas as pd
pd.set_option('max_columns', 200)
# Visualization packages
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
import warnings
warnings.filterwarnings("ignore")
EDA: House Prices Data From Kaggle¶
For our dataset, we'll be using the house price dataset from King County, WA. The dataset is from Kaggle.
The task is to build a regression model to predict the price, based on different attributes. First, let's do some EDA.
# Load the dataset
house_df = pd.read_csv('../data/kc_house_data.csv')
house_df = house_df.sample(frac=1, random_state=42)[0:4000]
print(house_df.shape)
print(house_df.dtypes)
house_df.head()
Now let's check for null values and look at the datatypes within the dataset.
house_df.info()
house_df.describe()
Let's choose a subset of columns here. NOTE: The way I'm selecting columns here is not principled and is just for convenience. In your homework assignments (and in real life), we expect you to choose columns more rigorously.
bedroomsbathroomssqft_livingsqft_lotfloorssqft_abovesqft_basementlatlongprice: Our response variable
cols_of_interest = ['bedrooms', 'bathrooms', 'sqft_living', 'sqft_lot', 'floors', 'sqft_above', 'sqft_basement',
'lat', 'long', 'price']
house_df = house_df[cols_of_interest]
# Convert house price to 1000s of dollars
house_df['price'] = house_df['price']/1000
Let's see how the response variable (price) is distributed
fig, ax = plt.subplots(figsize=(12,5))
ax.hist(house_df['price'], bins=100)
ax.set_title('Histogram of house price (in 1000s of dollars)');
# This takes a bit of time but is worth it!!
# sns.pairplot(house_df);
Train-Validation-Test Split¶
Up until this point, we have only had a train-test split. Why are we introducing a validation set? What's the point?
This is the general idea:
- Training Set: Data you have seen. You train different types of models with various different hyper-parameters and regularization parameters on this data.
- Validation Set: Used to compare different models. We use this step to tune our hyper-parameters i.e. find the optimal set of hyper-parameters (such as $k$ for k-NN or our $\beta_i$ values or number of degrees of our polynomial for linear regression). Pick your best model here.
- Test Set: Using the best model from the previous step, simply report the score e.g. R^2 score, MSE or any metric that you care about, of that model on your test set. DON'T TUNE YOUR PARAMETERS HERE!. Why, I hear you ask? Because we want to know how our model might do on data it hasn't seen before. We don't have access to this data (because it may not exist yet) but the test set, which we haven't seen or touched so far, is a good way to mimic this new data.
Let's do 60% train, 20% validation, 20% test for this dataset.
from sklearn.model_selection import train_test_split
# first split the data into a train-test split and don't touch the test set yet
train_df, test_df = train_test_split(house_df, test_size=0.2, random_state=42)
# next, split the training set into a train-validation split
# the test-size is 0.25 since we are splitting 80% of the data into 20% and 60% overall
train_df, val_df = train_test_split(train_df, test_size=0.25, random_state=42)
print('Train Set: {0:0.2f}%'.format(100*train_df.size/house_df.size))
print('Validation Set: {0:0.2f}%'.format(100*val_df.size/house_df.size))
print('Test Set: {0:0.2f}%'.format(100*test_df.size/house_df.size))
Modeling¶
In the last section, we went over the mechanics of Multiple Linear Regression and created models that had interaction terms and polynomial terms. Specifically, we dealt with the following sorts of models.
$$ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_M x_M $$Let's adopt a similar process here and get a few different models.
Creating a Design Matrix¶
From our model setup in the equation in the previous section, we obtain the following:
$$ Y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}, \quad X = \begin{bmatrix} x_{1,1} & x_{1,2} & \dots & x_{1,M} \\ x_{2,1} & x_{2,2} & \dots & x_{2,M} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n,1} & x_{n,2} & \dots & x_{n,M} \\ \end{bmatrix}, \quad \beta = \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_M \end{bmatrix}, \quad \epsilon = \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{bmatrix}, $$$X$ is an n$\times$M matrix: this is our design matrix, $\beta$ is an M-dimensional vector (an M$\times$1 matrix), and $Y$ is an n-dimensional vector (an n$\times$1 matrix). In addition, we know that $\epsilon$ is an n-dimensional vector (an n$\times$1 matrix).
X = train_df[cols_of_interest]
y = train_df['price']
print(X.shape)
print(y.shape)
Scaling our Design Matrix¶
Warm-Up Exercise¶
Warm-Up Exercise: for which of the following do the units of the predictors matter (e.g., trip length in minutes vs seconds; temperature in F or C)? A similar question would be: for which of these models do the magnitudes of values taken by different predictors matter?
(We will go over Ridge and Lasso Regression in greater detail later)
- k-NN (Nearest Neighbors regression)
- Linear regression
- Lasso regression
- Ridge regression
Solutions
- kNN: yes. Scaling affects distance metric, which determines what "neighbor" means
- Linear regression: no. Multiply predictor by $c$ -> divide coef by $c$.
- Lasso: yes: If we divided coef by $c$, then corresponding penalty term is also divided by $c$.
- Ridge: yes: Same as Lasso, except penalty divided by $c^2$.
Standard Scaler (Standardization)¶
Here's the scikit-learn implementation of the standard scaler. What is it doing though? Hint: you may have seen this in STAT 110 or another statistics course multiple times.
$$ z = \frac{x-\mu}{\sigma} $$In the above setup:
- $z$ is the standardized variable
- $x$ is the variable before standardization
- $\mu$ is the mean of the variable before standardization
- $\sigma$ is the standard deviation of the variable before standardization
Let's see an example of how this works:
from sklearn.preprocessing import StandardScaler
x = house_df['sqft_living']
mu = x.mean()
sigma = x.std()
z = (x-mu)/sigma
# reshaping x to be a n by 1 matrix since that's how scikit learn likes data for standardization
x_reshaped = np.array(x).reshape(-1,1)
z_sklearn = StandardScaler().fit_transform(x_reshaped)
# Plotting the histogram of the variable before standardization
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(24,5))
ax = ax.ravel()
ax[0].hist(x, bins=100)
ax[0].set_title('Histogram of sqft_living before standardization')
ax[1].hist(z, bins=100)
ax[1].set_title('Manually standardizing sqft_living')
ax[2].hist(z_sklearn, bins=100)
ax[2].set_title('Standardizing sqft_living using scikit learn');
# making things a dataframe to check if they work
pd.DataFrame({'x': x, 'z_manual': z, 'z_sklearn': z_sklearn.flatten()}).describe()
Min-Max Scaler (Normalization)¶
Here's the scikit-learn implementation of the standard scaler. What is it doing though?
$$ x_{new} = \frac{x-x_{min}}{x_{max}-x_{min}} $$In the above setup:
- $x_{new}$ is the normalized variable
- $x$ is the variable before normalized
- $x_{max}$ is the max value of the variable before normalization
- $x_{min}$ is the min value of the variable before normalization
Let's see an example of how this works:
from sklearn.preprocessing import MinMaxScaler
x = house_df['sqft_living']
x_new = (x-x.min())/(x.max()-x.min())
# reshaping x to be a n by 1 matrix since that's how scikit learn likes data for normalization
x_reshaped = np.array(x).reshape(-1,1)
x_new_sklearn = MinMaxScaler().fit_transform(x_reshaped)
# Plotting the histogram of the variable before normalization
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(24,5))
ax = ax.ravel()
ax[0].hist(x, bins=100)
ax[0].set_title('Histogram of sqft_living before normalization')
ax[1].hist(x_new, bins=100)
ax[1].set_title('Manually normalizing sqft_living')
ax[2].hist(x_new_sklearn, bins=100)
ax[2].set_title('Normalizing sqft_living using scikit learn');
# making things a dataframe to check if they work
pd.DataFrame({'x': x, 'x_new_manual': x_new, 'x_new_sklearn': x_new_sklearn.flatten()}).describe()
The million dollar question
Should I standardize or normalize my data? This, this and this are useful resources that I highly recommend. But in a nutshell, what they say is the following:
Pros of Normalization
- Normalization (which makes your data go from 0-1) is widely used in image processing and computer vision, where pixel intensities are non-negative and are typically scaled from a 0-255 scale to a 0-1 range for a lot of different algorithms.
- Normalization is also very useful in neural networks (which we will see later in the course) as it leads to the algorithms converging faster.
- Normalization is useful when your data does not have a discernible distribution and you are not making assumptions about your data's distribution.
Pros of Standardization
- Standardization maintains outliers (do you see why?) whereas normalization makes outliers less obvious. In applications where outliers are useful, standardization should be done.
- Standardization is useful when you assume your data comes from a Gaussian distribution (or something that is approximately Gaussian).
Some General Advice
- We learn parameters for standardization ($\mu$ and $\sigma$) and for normalization ($x_{min}$ and $x_{max}$). Make sure these parameters are learned on the training set i.e use the training set parameters even when normalizing/standardizing the test set. In sklearn terms, fit your scaler on the training set and use the scaler to transform your test set and validation set (don't re-fit your scaler on test set data!).
- The point of standardization and normalization is to make your variables take on a more manageable scale. You should ideally standardize or normalize all your variables at the same time.
- Standardization and normalization is not always needed and is not an automatic thing you have to do on any data science homework!! Do so sparingly and try to justify why this is needed.
Interpreting Coefficients
A great quote from here
[Standardization] makes it so the intercept term is interpreted as the expected value of 𝑌𝑖 when the predictor values are set to their means. Otherwise, the intercept is interpreted as the expected value of 𝑌𝑖 when the predictors are set to 0, which may not be a realistic or interpretable situation (e.g. what if the predictors were height and weight?)
Standardizing our Design Matrix¶
features = ['bedrooms', 'bathrooms', 'sqft_living', 'sqft_lot', 'floors', 'sqft_above', 'sqft_basement',
'lat', 'long']
X_train = train_df[features]
y_train = np.array(train_df['price']).reshape(-1,1)
X_val = val_df[features]
y_val = np.array(val_df['price']).reshape(-1,1)
X_test = test_df[features]
y_test = np.array(test_df['price']).reshape(-1,1)
scaler = StandardScaler().fit(X_train)
# This converts our matrices into numpy matrices
X_train_t = scaler.transform(X_train)
X_val_t = scaler.transform(X_val)
X_test_t = scaler.transform(X_test)
# Making the numpy matrices pandas dataframes
X_train_df = pd.DataFrame(X_train_t, columns=features)
X_val_df = pd.DataFrame(X_val_t, columns=features)
X_test_df = pd.DataFrame(X_test_t, columns=features)
display(X_train_df.describe())
display(X_val_df.describe())
display(X_test_df.describe())
scaler = StandardScaler().fit(y_train)
y_train = scaler.transform(y_train)
y_val = scaler.transform(y_val)
y_test = scaler.transform(y_test)
One-Degree Polynomial Model¶
import statsmodels.api as sm
from statsmodels.regression.linear_model import OLS
model_1 = OLS(np.array(y_train).reshape(-1,1), sm.add_constant(X_train_df)).fit()
model_1.summary()
Two-Degree Polynomial Model¶
def add_square_terms(df):
df = df.copy()
cols = df.columns.copy()
for col in cols:
df['{}^2'.format(col)] = df[col]**2
return df
X_train_df_2 = add_square_terms(X_train)
X_val_df_2 = add_square_terms(X_val)
# Standardizing our added coefficients
cols = X_train_df_2.columns
scaler = StandardScaler().fit(X_train_df_2)
X_train_df_2 = pd.DataFrame(scaler.transform(X_train_df_2), columns=cols)
X_val_df_2 = pd.DataFrame(scaler.transform(X_val_df_2), columns=cols)
print(X_train_df.shape, X_train_df_2.shape)
# Also check using the describe() function that the mean and standard deviations are the way we want them
X_train_df_2.head()
model_2 = OLS(np.array(y_train).reshape(-1,1), sm.add_constant(X_train_df_2)).fit()
model_2.summary()
Three-Degree Polynomial Model¶
# generalizing our function from above
def add_square_and_cube_terms(df):
df = df.copy()
cols = df.columns.copy()
for col in cols:
df['{}^2'.format(col)] = df[col]**2
df['{}^3'.format(col)] = df[col]**3
return df
X_train_df_3 = add_square_and_cube_terms(X_train)
X_val_df_3 = add_square_and_cube_terms(X_val)
# Standardizing our added coefficients
cols = X_train_df_3.columns
scaler = StandardScaler().fit(X_train_df_3)
X_train_df_3 = pd.DataFrame(scaler.transform(X_train_df_3), columns=cols)
X_val_df_3 = pd.DataFrame(scaler.transform(X_val_df_3), columns=cols)
print(X_train_df.shape, X_train_df_3.shape)
# Also check using the describe() function that the mean and standard deviations are the way we want them
X_train_df_3.head()
model_3 = OLS(np.array(y_train).reshape(-1,1), sm.add_constant(X_train_df_3)).fit()
model_3.summary()
N-Degree Polynomial Model¶
# generalizing our function from above
def add_higher_order_polynomial_terms(df, N=7):
df = df.copy()
cols = df.columns.copy()
for col in cols:
for i in range(2, N+1):
df['{}^{}'.format(col, i)] = df[col]**i
return df
N = 8
X_train_df_N = add_higher_order_polynomial_terms(X_train,N)
X_val_df_N = add_higher_order_polynomial_terms(X_val,N)
# Standardizing our added coefficients
cols = X_train_df_N.columns
scaler = StandardScaler().fit(X_train_df_N)
X_train_df_N = pd.DataFrame(scaler.transform(X_train_df_N), columns=cols)
X_val_df_N = pd.DataFrame(scaler.transform(X_val_df_N), columns=cols)
print(X_train_df.shape, X_train_df_N.shape)
# Also check using the describe() function that the mean and standard deviations are the way we want them
X_train_df_N.head()
model_N = OLS(np.array(y_train).reshape(-1,1), sm.add_constant(X_train_df_N)).fit()
model_N.summary()
You can also create a model with interaction terms or any other higher order polynomial term of your choice. Note: Can you see how creating a function that takes in a dataframe and a degree and creates polynomial terms up until that degree can be useful? This is what we have you do in your homework!
from sklearn.metrics import r2_score
x = [1,2,3,N]
models = [model_1, model_2, model_3, model_N]
X_trains = [X_train_df, X_train_df_2, X_train_df_3, X_train_df_N]
X_vals = [X_val_df, X_val_df_2, X_val_df_3, X_val_df_N]
r2_train = []
r2_val = []
for i,model in enumerate(models):
y_pred_tra = model.predict(sm.add_constant(X_trains[i]))
y_pred_val = model.predict(sm.add_constant(X_vals[i]))
r2_train.append(r2_score(y_train, y_pred_tra))
r2_val.append(r2_score(y_val, y_pred_val))
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, r2_train, 'o-', label=r'Training $R^2$')
ax.plot(x, r2_val, 'o-', label=r'Validation $R^2$')
ax.set_xlabel('Number of degree of polynomial')
ax.set_ylabel(r'$R^2$ score')
ax.set_title(r'$R^2$ score vs polynomial degree')
ax.legend();
We notice a big difference between training and validation R^2 scores: seems like we are overfitting. Introducing: regularization.
What about Multicollinearity?¶
There's seemingly a lot of multicollinearity in the data. Take a look at this warning that we got when showing our summary for our polynomial models:
What is multicollinearity? Why do we have it in our dataset? Why is this a problem?
Does regularization help solve the issue of multicollinearity?
What does Regularization help with?¶
We have some pretty large and extreme coefficient values in our most recent models. These coefficient values also have very high variance. We can also clearly see some overfitting to the training set. In order to reduce the coefficients of our parameters, we can introduce a penalty term that penalizes some of these extreme coefficient values. Specifically, regularization helps us:
- Avoid overfitting. Reduce features that have weak predictive power.
- Discourage the use of a model that is too complex
Big Idea: Reduce Variance by Increasing Bias¶
Image Source: here
Ridge Regression¶
Ridge Regression is one such form of regularization. In practice, the ridge estimator reduces the complexity of the model by shrinking the coefficients, but it doesn’t nullify them. We control the amount of regularization using a parameter $\lambda$. NOTE: sklearn's ridge regression package represents this $\lambda$ using a parameter alpha. In Ridge Regression, the penalty term is proportional to the L2-norm of the coefficients.
Lasso Regression¶
Lasso Regression is another form of regularization. Again, we control the amount of regularization using a parameter $\lambda$. NOTE: sklearn's lasso regression package represents this $\lambda$ using a parameter alpha. In Lasso Regression, the penalty term is proportional to the L1-norm of the coefficients.
Some Differences between Ridge and Lasso Regression¶
- Since Lasso regression tend to produce zero estimates for a number of model parameters - we say that Lasso solutions are sparse - we consider to be a method for variable selection.
- In Ridge Regression, the penalty term is proportional to the L2-norm of the coefficients whereas in Lasso Regression, the penalty term is proportional to the L1-norm of the coefficients.
- Ridge Regression has a closed form solution! Lasso Regression does not. We often have to solve this iteratively. In the sklearn package for Lasso regression, there is a parameter called
max_iterthat determines how many iterations we perform.
Why Standardizing Variables was not a waste of time¶
Lasso regression puts constraints on the size of the coefficients associated to each variable. However, this value will depend on the magnitude of each variable. It is therefore necessary to standardize the variables.
Let's use Ridge and Lasso to regularize our degree N polynomial¶
Exercise: Play around with different values of alpha. Notice the new $R^2$ value and also the range of values that the predictors take in the plot.
from sklearn.linear_model import Ridge
# some values you can try out: 0.01, 0.1, 0.5, 1, 5, 10, 20, 40, 100, 200, 500, 1000, 10000
alpha = 100
ridge_model = Ridge(alpha=alpha).fit(X_train_df_N, y_train)
print('R squared score for our original OLS model: {}'.format(r2_val[-1]))
print('R squared score for Ridge with alpha={}: {}'.format(alpha, ridge_model.score(X_val_df_N,y_val)))
fig, ax = plt.subplots(figsize=(18,8), ncols=2)
ax = ax.ravel()
ax[0].hist(model_N.params, bins=10, alpha=0.5)
ax[0].set_title('Histogram of predictor values for Original model with N: {}'.format(N))
ax[0].set_xlabel('Predictor values')
ax[0].set_ylabel('Frequency')
ax[1].hist(ridge_model.coef_.flatten(), bins=20, alpha=0.5)
ax[1].set_title('Histogram of predictor values for Ridge Model with alpha: {}'.format(alpha))
ax[1].set_xlabel('Predictor values')
ax[1].set_ylabel('Frequency');
from sklearn.linear_model import Lasso
# some values you can try out: 0.00001, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.5, 1, 2, 5, 10, 20
alpha = 0.01
lasso_model = Lasso(alpha=alpha, max_iter = 1000).fit(X_train_df_N, y_train)
print('R squared score for our original OLS model: {}'.format(r2_val[-1]))
print('R squared score for Lasso with alpha={}: {}'.format(alpha, lasso_model.score(X_val_df_N,y_val)))
fig, ax = plt.subplots(figsize=(18,8), ncols=2)
ax = ax.ravel()
ax[0].hist(model_N.params, bins=10, alpha=0.5)
ax[0].set_title('Histogram of predictor values for Original model with N: {}'.format(N))
ax[0].set_xlabel('Predictor values')
ax[0].set_ylabel('Frequency')
ax[1].hist(lasso_model.coef_.flatten(), bins=20, alpha=0.5)
ax[1].set_title('Histogram of predictor values for Lasso Model with alpha: {}'.format(alpha))
ax[1].set_xlabel('Predictor values')
ax[1].set_ylabel('Frequency');
Model Selection and Cross-Validation¶
Here's our current setup so far:
So we try out 10,000 different models on our validation set and pick the one that's the best? No! Since we could also be overfitting the validation set!
One solution to the problems raised by using a single validation set is to evaluate each model on multiple validation sets and average the validation performance. This is the essence of cross-validation!
Image source: here
from sklearn.linear_model import RidgeCV
from sklearn.linear_model import LassoCV
alphas = (0.001, 0.01, 0.1, 10, 100, 1000, 10000)
# Let us do k-fold cross validation
k = 4
fitted_ridge = RidgeCV(alphas=alphas).fit(X_train_df_N, y_train)
fitted_lasso = LassoCV(alphas=alphas).fit(X_train_df_N, y_train)
print('R^2 score for our original OLS model: {}\n'.format(r2_val[-1]))
ridge_a = fitted_ridge.alpha_
print('Best alpha for ridge: {}'.format(ridge_a))
print('R^2 score for Ridge with alpha={}: {}\n'.format(ridge_a, fitted_ridge.score(X_val_df_N,y_val)))
lasso_a = fitted_lasso.alpha_
print('Best alpha for lasso: {}'.format(lasso_a))
print('R squared score for Lasso with alpha={}: {}'.format(lasso_a, fitted_lasso.score(X_val_df_N,y_val)))
We can also look at the coefficients of our CV models.
Final Step: report the score on the test set for the model you have chosen to be the best.