16.4. Regularization

We just saw how cross-validation can help find a dimension for a fitted model that balances under- and overfitting. Rather than selecting the dimension of the model, we can build a model with all of the features, but restrict the size of the coefficients. We keep from overfitting by adding to the MSE a penalty term on the size of the coefficients. The penalty, called a regularization term, is \( \lambda \sum_{j = 1}^{p} \theta_j^2 \). We fit the model by minimizing the combination of mean squared error plus this penalty:

\[ \frac{1}{n} \sum_{i=1}^{n}(y_i - \mathbf{x}_i \boldsymbol{\theta})^2 ~+~ \lambda \sum_{j = 1}^{p} \theta_j^2 \]

When the regularization parameter, \(\lambda\), is large, it penalizes large coefficients. (We typically choose it by cross-validation.)

Penalizing the square of the coefficients is called \( L_2 \) regularization, or ridge regression. Another popular regularization penalizes the absolute size of the coefficients:

\[ \begin{aligned} \frac{1}{n} \sum_{i=1}^{n}(y_i - \mathbf{x}_i \boldsymbol{\theta})^2 ~+~ \lambda \sum_{j = 1}^{p} |\theta_j| \end{aligned} \]

This \( L_1 \) regularized linear model is also called lasso regression (lasso stands for Least Absolute Shrinkage and Selection Operator).

To get an idea about how regularization works, let’s think about the extreme cases: when \( \lambda \) is really large and when it’s close to 0 (\( \lambda \) is never negative). With a big regularization parameter, the coefficients are heavily penalized, so they shrink. On the other hand, when \( \lambda \) is tiny, the coefficients aren’t restricted. In fact, when \( \lambda \) is 0, we’re back in the world of ordinary least squares. A couple of issues crop up when we think about controlling the size of the coefficients through regularization:

• We do not want to regularize the intercept term. This way, a large penalty fits a constant model.

• When features have very different scales, the penalty can impact them differently, with large-valued features being penalized more than others. To avoid this, we standardize all of the features to have mean 0 and variance 1 before fitting the model.

Let’s look at an example with 35 features.

16.5. Example: A Market Analysis

A market research project for a pharmaceutical company wants to model consumer interest in purchasing a cold sore health-care product. They gather data from 1,023 consumers. Each consumer is asked to rate on a 10-point scale 35 factors according to whether the factor matters to them when considering purchasing a cold sore treatment. They also rate their interest in purchasing the product.1

We begin by reading in the data:

ma_df = pd.read_csv('data/market-analysis.csv')

Table 16.1 lists the 35 factors and provides their correlation to the outcome, their interest in purchasing the product.

Table 16.1 Codebook

	Corr	Description		Corr	Description
x1	0.70	provides soothing relief	x19	0.54	has a non-messy application
x2	0.58	moisturizes cold sore blister	x20	0.70	good for any stage of a cold
x3	0.69	provides long-lasting relief	x21	0.49	easy to apply/take
x4	0.70	provides fast-acting relief	x22	0.52	package keeps from contamination
x5	0.72	shortens duration of a cold	x23	0.57	easy to dispense a right amount
x6	0.68	stops the virus from spreading	x24	0.63	worth the price it costs
x7	0.67	dries up cold sore	x25	0.57	recommended most by pharamacists
x8	0.72	heals fast	x26	0.54	recommended by doctors
x9	0.72	penetrates deep	x27	0.54	FDA approved
x10	0.65	relieves pain	x28	0.64	a brand I trust
x11	0.61	prevents cold	x29	0.60	clinically proven
x12	0.73	prevents from getting worse	x30	0.68	a brand I would recommend
x13	0.57	medicated	x31	0.74	an effective treatment
x14	0.61	prescription strength	x32	0.37	portable
x15	0.63	repairs damaged skin	x33	0.37	discreet packaging
x16	0.67	blocks virus from spreading	x34	0.55	helps conceal cold sores
x17	0.42	contains SPF	x35	0.63	absorbs quickly
x18	0.57	non-irritating

Based on their labels alone, some of these 35 features appear to measure similar aspects of desirability. We can compute the correlations between the explanatory variables to confirm this:

	y	x1	x2	x3	x4	x5	x6	...	x29	x30	x31	x32	x33	x34	x35
y	1.00	0.70	0.58	0.69	0.70	0.72	0.68	...	0.60	0.68	0.74	0.37	0.37	0.55	0.63
x1	0.70	1.00	0.69	0.74	0.81	0.75	0.73	...	0.63	0.68	0.79	0.42	0.33	0.57	0.71
x2	0.58	0.69	1.00	0.63	0.65	0.63	0.63	...	0.60	0.62	0.66	0.43	0.36	0.55	0.63
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
x33	0.37	0.33	0.36	0.34	0.32	0.32	0.37	...	0.39	0.41	0.35	0.34	1.00	0.35	0.36
x34	0.55	0.57	0.55	0.67	0.60	0.59	0.64	...	0.58	0.48	0.58	0.16	0.35	1.00	0.54
x35	0.63	0.71	0.63	0.70	0.72	0.69	0.68	...	0.65	0.67	0.72	0.46	0.36	0.54	1.00

36 rows × 36 columns

We see, for example, that the last feature, “absorbs quickly,” is highly correlated with the first three: “provides soothing relief,” “moisturizes,” and “provides long-lasting relief.”

Before we start to fit regularized models, we set up the design matrix and outcome vector and split the data into train and test sets. We place 200 observations in the test set:

y = ma_df["y"]
X = ma_df.drop(columns=["y"])

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=200, random_state=42
)

As mentioned earlier, we need to standardize the features. We standardize the train set, and then, when we go to evaluate the model, we use the train set standardization on the test set. The StandardScaler method helps us with this process:

from sklearn.preprocessing import StandardScaler

scalerX = StandardScaler().fit(X_train) 
X_train_scaled = scalerX.transform(X_train)
X_test_scaled = scalerX.transform(X_test)

We confirm that the means of the 35 features in the train set are all 0 and the SDs are all 1:

np.allclose(X_train_scaled.mean(axis=0), 0)

True

np.allclose(X_train_scaled.std(axis=0), 1)

True

It’s important to note that this is not the case for the test set, because we use the averages and SDs from the train set to standardize the test set features:

X_test_scaled.mean(axis=0)

array([-0.17, -0.12, -0.15, ..., -0.12, -0.1 , -0.24])

To perform a lasso regression, we use the Lasso method in scikit-learn. Let’s see how the coefficients for the 35 features in the data frame change with \(\lambda\). We set up a range of values for the regularization parameter and fit the lasso for each (this parameter is referred to as alpha in Lasso):

from sklearn.linear_model import Lasso

coefs = []
mses = []
alphas = np.arange(0.01, 2, 0.01)

for a in alphas:
    model = Lasso(alpha=a)
    model.fit(X_train_scaled, y_train)
    coefs.append(model.coef_)
    mses.append(mean_squared_error(y_test, model.predict(X_test_scaled)))

For each feature, we can overlay a line plot of the coefficient against \( \lambda \):

As \( \lambda \) increases, the model fitting is penalized more heavily and many coefficients shrink to 0. The lefthand side of the plot shows greater model complexity and corresponds to small \( \lambda \). Notice that for \( \lambda = 0.5 \), nearly all of the coefficients are 0. The exceptions are, from largest to smallest: x12, x30, x31, x5, x7, x1, x9, x28, x20.

We can also plot the MSE as a function of \( \lambda \) to see how it changes with the increase in the penalty:

px.line(x=alphas, y=mses,
        labels={"x": "Lambda", "y": "MSE"},
        width=350, height=250)

When \( \lambda \) reaches about 1.25, the penalty is so large that the lasso regression is simply fitting a constant model to the data, and the MSE doesn’t change. We are again faced with a model selection question, but this time it’s in the form of deciding on \( \lambda \). We can use cross-validation to help us.

We use the LassoCV method to perform 5-fold cross-validation to choose \(\lambda\):

from sklearn.linear_model import LassoCV

lasso_cv_model = LassoCV(
    alphas=np.arange(0.01, 1, 0.01), cv=5, max_iter=100000
)

Notice that we have specified a maximum number of iterations because the minimization uses numerical optimization (see Chapter 20) to solve for the coefficients, and we have placed a cap on the number of iterations to run to reach the specified tolerance for the optimal parameter. We’re ready to fit the model on the train set:

lasso_best = lasso_cv_model.fit(X_train_scaled, y_train)

Cross-validation has chosen the following regularization parameter:

lasso_best.alpha_

0.04

Next, let’s compute the MSE for predictions in the test set using this cross-validated model:

y_test_pred = lasso_best.predict(X_test_scaled)
mean_squared_error(y_test, y_test_pred)

1.0957444655340487

Let’s find how many coefficients are not 0:

sum(np.abs(lasso_best.coef_) > 0)

16

We next demonstrate how the ridge regression coefficients change with the penalty. Ridge regression uses an \( L_2 \) penalty, but otherwise the implementation is similar to lasso regression. As before, we try a range of \( \lambda \) values and then plot the coefficients to see how they change:

from sklearn.linear_model import Ridge

coefsR = []
alphasR = np.arange(1, 1001, 25)

for a in alphasR:
    ridge = Ridge(alpha=a)
    ridge.fit(X_train_scaled, y_train)
    coefsR.append(ridge.coef_)

This plot shows a very different shape compared to the lasso coefficient plot. The coefficients do not disappear entirely. The plots do have an important similarity in that as the penalty increases, the coefficients shrink. In a way, with ridge regression, we are using all of the variables a little bit. For large \( \lambda \), many coefficients, while not 0, are quite small. For example, we can count the number of coefficients larger than 0.05 for \( \lambda = 600 \):

sum(abs(coefsR[24, :]) > 0.05)

13

Similar to lasso, many of the features have tiny coefficients.

Again, to select the best \( \lambda \) we can turn to cross-validation and use RidgeCV. We omit the details here because the approach is the same as LassoCV.

Using \( L_1 \) or \( L_2 \) regularization allows us to avoid model overfitting by penalizing large coefficients. \( L_1 \) regularization has the advantage of zeroing out coefficients; it performs feature selection by discarding the features with a coefficient of 0. This is particularly useful when working with high-dimensional data with many features. A model that only uses a few features to make a prediction runs much faster than a model that requires many calculations. Since unneeded features tend to increase model variance without decreasing bias, we can sometimes increase the accuracy of other models by using lasso regression to select a subset of features to use.

Note

Sometimes we prefer one type of regularization over the other because it maps more closely to the domain we are working with. For example, if we know that the phenomenon we are trying to model results from many small factors, we might prefer ridge regression because it won’t discard these factors. On the other hand, some outcomes result from a few highly influential features. We prefer lasso regression in these situations because it discards unneeded features. Both \(L_2\) and \(L_1\) serve as useful throttles to navigate between under- and overfitting.

Regularization, train-test split, and cross-validation all have the goal of reducing over-fitting. The problem with over-fitting comes from using data to both fit a model and estimate the model’s error in predicting new observations. In the next section, we provide further intuition for this idea.

---

1

Stan Lipovetsky and W. Michael Conklin (2015) Predictor relative importance and matching regression parameters. Journal of Applied Statistics, 42:5, 1017-1031.