Regularization Techniques

Similar to the backwards elimination algorithm and the forward selection algorithm, regularization techniques are another way to go about determining which explanatory variables may lead to overfitting and/or don't bring enough "predictive power" to the model, and should therefore be left out.

Introducing a Penalty Term into a Linear Regression

Similarity to Adjusted R^2 Metric

Similar to the adjusted R^2 metric, which we tried to maximize in our backwards elimination algorithm and the forward selection algorithm, these regularization techniques also introduce a type of penalty term that penalizes models that have too many explanatory variables that don't bring enough predictive power to the model.

In our adjusted R^2 equation, the penalty term was $(\frac{n-1}{n-p-1})$

$$Adjusted \ R^2=1-(\frac{SSE}{SST})\cdot(\frac{n-1}{n-p-1})$$

Adding More Slopes

Similarly, by increasing the number of slopes $p$, the adjusted R^2 will be encouraged to decrease. However, if we were to add a slope and the resulting decrease in the new model's error (SSE) were to be "large enough", then the adjusted R^2 would actually increase overall, and thus:

1. this resulting model would be deemed "more parsimonious",
2. this new slope would be considered bringing enough "predictive power" to the model, and
3. this new slope would not be considered overfitting.

However, unfortunately, the quest of trying to find the linear regression model with the highest adjusted R^2 in the backwards elimination algorithm and the forward selection algorithms involved having to fit multiple models, each time checking the adjusted R^2 of the test models to see if the adjusted R^2 value got any better.

Objective Function

Modifying the Objective Function with a Penalty Term

So, rather than having to fit multiple basic linear regression models (and thus having to optimize multiple basic linear regression objective functions in search of the model that gives you the highest adjusted R^2), wouldn't it be great if we could just incorporate a penalty term directly into linear regression objective function?

This is the idea behind how a regularized linear regression model works. In a regularized linear regression model, we take the objective function from our basic (ie. nonregularized) linear regression model and add a type of penalty term.

Goal of the Penalty Term

So what specifically should we set this penalty term to be? Remember, the overall goal of our feature selection techniques is to determine which explanatory variables may cause our model to overfit to the training dataset, and therefore to determine which explanatory should we leave out of the model.

Goal 1: Interpretation of what to leave out

Question: So using just our regularized linear regression optimization problem above, what kind of optimal slope results could we get back that signals to us that the corresponding explanatory variable (or indicator variable) should be left out of the model?

Answer: If a returned optimal slope value is set $\hat{\beta}^*_i=0$, this serves as pretty clear signal that the model thinks the corresponding explanatory variable (or indicator variable) $x_i$ should be left out!

Goal 2: Leave out variables that lead to overfitting

Questions: Now that we know how our fitted regularized linear regression model will indicate to us which variables should be left out, how should we set up this penalty term to encourage the model to only indicate leaving out the variables that lead to "too much" overfitting? Furthermore, how do we define "too much" overfitting?

Answers: In sections 8.2, 8.3, and 8.4 we'll talk about three of the most popular penalty term functions that are used that attempt to satisfy these two goals.

• 8.2. L1 Penalty Term (ie. LASSO Regression)
• 8.3. L2 Penalty Term (ie. Ridge Regression)
• 8.4 Combination of L1 and L2 Penalty Terms (ie. Elastic Net Regression)

LASSO Regression (ie. L1 Penalty Term)

In the regularized LASSO linear regression optimization problem shown below, the penalty term that is added is

We call this a L1 penalty term, because what we call the L1 norm of a given list (ie. vector) of values is the sum of the absolute values of the values in this list. (The LASSO stands for Least Absolute Shrinkage and Selection Operator).

Also in this penalty term, we introduce a parameter $\lambda\geq 0$ that we need to choose ahead of time. This $\lambda$ value can be anything you want as long as it is not negative.

The Ideal Scenario

To help us see why adding this penalty term would help to reduce overfitting, let's think about what the absolute best intercept and slope solutions to this optimization would be. We can see that the full objective function now is comprised of the sum of two terms.

The Penalty Term Balancing Act

However, note that this "ideal scenario" that we just described is unrealistic. It's unlikely that we'll be able to fit a model where ALL the slopes $\hat{\beta}_j=0$ AND there is absolutely NO model error (ie. $SSE=0$). In addition, what we actually see is that these terms end up competing with one another. That is, choosing intercepts and slopes that lead to a decrease of one term, often leads to an increase in the other term, and vice versa.

The Impact of $\lambda$

Now let's think about the impact that our preselected $\lambda$ parameter has on this balancing act.

If we were to set $\lambda$ to be as low as possible (=0), then we can see that this penalty term would actually have no impact at all on the linear regression intercepts and slopes that were chosen. With $\lambda=0$, then we can see that our objective function actually just becomes the same as our nonregularized linear regression objective function. Thus in this case, we are fitting just a basic linear regression model. The slopes values face no penalty for being as large as they need to be in order to minimize the SSE.

Which slope magnitudes get to stay large?

But does this mean that all of the $|\hat{\beta}_j|$ terms need to decrease equally? Or is it possible that the optimization problem will cleverly keep some slopes high that correspond to the variables that contribute the most predictive power to the model thereby leading to much larger reductions to the SSE? And could it also keep other slopes low that correspond to the variables that do not contribute a lot independently to the predictive power of the model (or in other words, may be causing the model to overfit)? Luckily, in practice we see that the answer to these questions is yes!

Thus the resulting optimal slopes that are set low, are the ones that bring the least predictive power to the model and thus are the most likely variable to overfit the model.

Clearer Slope Interpretation with LASSO Regression

So what we learned is that the less predictive power that a given variable independently brings to a regression model, the more likely it's corresponding slope is going to be set low. But then we might ask: how do we know if a resulting slope value is low enough to warrant leaving it out of the model because it does not bring enough predictive power to the model?

One of the benefits of using the L1 penalty term in regression regularization (as opposed to the other two penalty terms that we'll discuss in 8.3 and 8.4) is how the absolute value function interacts with our Calculus techniques that are used to solve for these optimal intercept and slope solutions. In general, if a slope is encouraged by the optimization problem to be comparatively small, then the optimal solution values for these small slopes can actually force these small slope values to be set exactly equal to 0.

These small slope values being set exactly equal to 0, gets rid of a lot of the interpretation ambiguity. If a slope is set to be 0, then the LASSO regression model is suggesting that this slope's corresponding variable can be left out of the model because it is not bringing enough predictive power to the model, based on the $\lambda$ value that you chose.

Nonregularized Linear Regression for Predicting Tumor Size - for Comparison

Before we fit a LASSO linear regression model with our scaled training features matrix and target array, let's also fit a basic nonregularized linear regression model for comparison purposes. To do so, let's use the LinearRegression() that we used in Data Science Discovery.

Recall that we first instantiate an object with this LinearRegression() function. And then we use the .fit() function with our features matrix ($X_{train}$) and target array ($y_{train}$) to actually fit the model.

from sklearn.linear_model import LinearRegression

lin_reg_mod = LinearRegression() 
lin_reg_mod.fit(X_train, y_train)

LinearRegression()

LinearRegression()

And let's display the slopes for this fitted linear regression model below by using the .coef_ attribute. Because we scaled our gene explanatory variables, we are technically able to interpret the corresponding slope magnitudes as an indication of how important that particular corresponding variable is when it comes to predicting the tumor size in this particular model, in the context of the other existing variables in the model.

df_slopes = pd.DataFrame(lin_reg_mod.coef_.T, columns = ['lin_reg_mod'], index=X_train.columns)
df_slopes.sort_values(by=['lin_reg_mod'])

Issue #1: Biased Slopes due to Multicollinearity

However, remember that our features matrix showed many pairs of variables that were collinear. Thus, it is possible that our slope interpretations could still be biased due to this multicollinearity.

Issue #2: Could we be Overfitting?

Remember that one of our other goals is to use our predictive model to attain good predictions for new datasets. Let's assess the test dataset $R^2$ of this linear regression model before moving on to our LASSO linear regression model.

We can use the .score() function along with our test features matrix $X_{test}$ and target array $y_{test}$ to calculate this test R^2 of the model.

lin_reg_mod.score(X_test, y_test)

-0.3071595344334954

Negative R^2?

Interestingly enough we got a $R^2$ value that may not be what we were expecting. The R^2 of the test dataset is negative for this model! If the R^2 value was between 0 and 1, this value would represent the percent of tumor size variability that is explained by the model.

However, when the R^2 value is negative, one downside is that we lose the interpretability. Another huge downside, is that this implies that this model is a not a good fit for the test dataset. In general, the lower the R^2 value, the worse the fit.

LASSO Linear Regression for Predicting Tumor Size

Now, let's actually apply our LASSO Linear regression model to this dataset.

Recall that our LASSO linear regression objective function required for us to preselect a value for $\lambda$. Let's try out a couple different values for $\lambda\in [0.025, 0.05, 0.1, 0.5]$ and see how this impacts our results and research goals.

We can use the Lasso() function to fit a linear LASSO regression model. This Lasso() function operates in a similar way to the LinearRegression() function. First, we create a Lasso() object with specified parameters. Then we take this object and .fit() it with the features matrix and the target array.

• The alpha parameter represents the $\lambda$ parameter.
• The algorithm that is used to attempt to find the optimal intercept and slopes to this LASSO regression objective function may not always arrive at an optimal solution in a fixed number of steps. Therefore, in the interest of time, we can set the max_iter parameter to be the maximum number of iterations that this algorithm should take before just stopping early and returning the best solution that it currently has.

Fitting Four LASSO Models

Model 1: $\lambda=0.025$

from sklearn.linear_model import Lasso
lasso_mod_025 = Lasso(alpha=0.025, max_iter=1000)
lasso_mod_025.fit(X_train, y_train)

Lasso(alpha=0.025)

Lasso(alpha=0.025)

Model 2: $\lambda=.05$

lasso_mod_05 = Lasso(alpha=.05, max_iter=1000)
lasso_mod_05.fit(X_train, y_train)

Lasso(alpha=0.05)

Lasso(alpha=0.05)

Model 3: $\lambda=.1$

lasso_mod_1 = Lasso(alpha=.1, max_iter=1000)
lasso_mod_1.fit(X_train, y_train)

Lasso(alpha=0.1)

Lasso(alpha=0.1)

Model 4: $\lambda=.5$

lasso_mod_5 = Lasso(alpha=.5, max_iter=1000)
lasso_mod_5.fit(X_train, y_train)

Lasso(alpha=0.5)

Lasso(alpha=0.5)

Inspecting the Slopes

Next, let's inspect and compare the slopes of the models that we've created so far. We can similarly using the .coef attribute to extract each of our 4 model slopes.

By comparing the slopes of the four LASSO linear regression models, we can see what we would expect. As we increase $\lambda$, the number of non-zero slopes decreased. In fact, when $\lambda=0.5$, there were only seven genes with non-zero slopes. This result may suggest that these genes are the most important when it comes to predicting breast tumor size.

df_slopes = pd.DataFrame({'lin_reg_mod': lin_reg_mod.coef_.T, 
                         'lasso_mod_025': lasso_mod_025.coef_.T,
                         'lasso_mod_05': lasso_mod_05.coef_.T,
                         'lasso_mod_1': lasso_mod_1.coef_.T,
                         'lasso_mod_5': lasso_mod_5.coef_.T}, index=X_train.columns)
df_slopes

Which model is best?

But could it be possible that some of these LASSO models have left out too many explanatory variables and thus may be underfitting?

Given that one of our goals is to build a model that will yield good predictions for new breast tumors, let's inspect the test dataset R^2 of each of these 4 LASSO models.

Nonregularized Linear Regression Model: Test Dataset R^2

lin_reg_mod.score(X_test, y_test)

-0.3071595344334954

LASSO Model 1: $\lambda=0.025$ Test Dataset R^2

lasso_mod_025.score(X_test, y_test)

0.1539346136658848

LASSO Model 2: $\lambda=0.05$ Test Dataset R^2

lasso_mod_05.score(X_test, y_test)

0.2786886369663709

LASSO Model 3: $\lambda=0.1$ Test Dataset R^2

lasso_mod_1.score(X_test, y_test)

0.2490838650434536

LASSO Model 4: $\lambda=0.5$ Test Dataset R^2

lasso_mod_5.score(X_test, y_test)

0.1445914387347873

We can see that the linear regression model that achieves the best test dataset R^2=0.279 is the LASSO model that used $\lambda=0.05$. Ideally, a more complete analysis would try out many other values of $\lambda$ to see if we can find a test dataset R^2 that was even higher.

Can we trust the non-zero slope interpretations in our chosen LASSO model?

Ideally, we'd be able to trust our interpretations of our non-zero slopes in this resulting LASSO model that we chose as well. Recall that many of our 50 gene explanatory variables in the nonregularized linear regression model were collinear with each other, which should make us more skeptical about how we interpret the resulting slopes.

If we were to remove the explanatory variables that had zero slopes in our chosen LASSO model, would the remaining explanatory variables still be collinear?

Below is a list of the 37 gene names that had non-zero slopes in our best LASSO model.

nonzero_slope_genes = df_slopes[np.abs(df_slopes['lasso_mod_05'])>0.00001].index
nonzero_slope_genes

Index(['X159', 'X960', 'X980', 'X986', 'X1023', 'X1028', 'X1064', 'X1092',
'X1109', 'X1141', 'X1169', 'X1173', 'X1179', 'X1193', 'X1203', 'X1206',
'X1208', 'X1219', 'X1264', 'X1272', 'X1329', 'X1351', 'X1362', 'X1416',
'X1417', 'X1430', 'X1444', 'X1470', 'X1506', 'X1514', 'X1529', 'X1563',
'X1595', 'X1609', 'X1637', 'X1656', 'X1683'],
dtype='object')

And below is the correlation matrix of each of these 37 remaining genes.

X_train_remaining = X_train[nonzero_slope_genes]
X_train_remaining.corr()

sns.heatmap(X_train_remaining.corr(),vmin=-1, vmax=1, cmap='RdBu')
plt.show()

Unfortunately, there is still a high degree of multicollinearity in the remaining explanatory variables. Therefore, we should also be cautious when interpreting the slopes of our chosen LASSO model as well.

Ridge Regression (ie. L2 Penalty Term)

Similar to LASSO regression, we introduce a parameter $\lambda\geq 0$ that we must preselect.

The Ideal Scenario

Similar to the LASSO linear regression, the absolute best intercept and slope solutions to this optimization would be ones in which:

The Penalty Term Balancing Act

However, similarly, this "ideal scenario" that we just described is unrealistic. We also see that choosing intercepts and slopes that lead to a decrease of one term, often leads to an increase in the other term, and vice versa.

The Impact of $\lambda$

For similar reasons described when we used $\lambda$ in the LASSO linear regression, we see the following impact of $\lambda$ on our results.

$\lambda$ is Set Low

$\lambda$ is Set High

Which slope magnitudes get to stay large?

Similar to LASSO linear regression, the resulting optimal slopes that are set "low", are the ones that bring the least predictive power to the model and thus are the most likely variable to overfit the model.

Less Clear Slope Interpretation with Ridge Regression

Unfortunately, one of the downsides of using the L2 penalty in ridge regression, is that unlike using the L1 penalty, slopes that are "comparatively small" will NOT be set exactly equal to 0. Therefore, the resulting slopes found with ridge regression provide much less of a clear indication as to which explanatory variables should be left out of the model because they don't bring enough predictive power to the model.

Without more careful inspection, it is not always clear if a ridge regression slope is small because it does not bring enough predictive power to the model, or if it is small for other reasons. This also begs the question: "what does it mean for a slope magnitude to be considered 'small'"?

Benefits of Ridge Regression

However, despite being harder to interpret, ridge regression does have some benefits over LASSO regression.

1. In the presence of multicollinearity, ridge regression slopes can be more trusted than those that would have been returned by an nonregularized linear regression model.

   For instance, suppose again we were trying to predict the height of a person $y$, using their left foot size $x_1$ and right foot size $x_2$. Again, it's likely that $x_1$ and $x_2$ are highly collinear. Thus, it's possible that your basic linear regression model may end up with slopes like $\hat{\beta}_1=1000$ and $\hat{\beta}_2=-999$ that are extremely large in magnitude with opposing signs, the negative slope $\hat{\beta}_2=-999$ being misleading in particular given that there is most likely a positive relationship between right foot size and height.

   Why would the model do this? Because $x_1$ and $x_2$ are so highly correlated, their two large opposing signs in this scenario could easily cancel this "largeness" out. Or in other words, their joint effort together in...
   $$\hat{\beta}_1x_1+\hat{\beta}_2x_2 = 100x_1-999(value\ very\ close\ to\ x_1)\approx 1x_1$$
   ...may actually end up leading to a much more realistic, modest predicted impact on $y$.

   Because ridge regression discourages slopes with large magnitudes, this type of scenario would be less likely to happen. Furthermore, because ridge regression does not "zero out" slopes like LASSO does, the predicted impact of these collinear variables on the response variable is more likely to get evenly distributed in the slopes.

2. Ridge Regression can be less compuationally complex in larger datasets
   Furthermore, due to some mathematical properties that the L2 norm has, but the L1 doesn't, finding the optimal slopes and intercept for a ridge regression model can take less time than it does for the LASSO model especially for larger datasets.

Linear Ridge Regression for Predicting Tumor Size

Now, let's actually apply our linear ridge regression model to this dataset.

Recall that our linear ridge regression objective function required for us to preselect a value for $\lambda$. Let's try out the following values $\lambda\in [1,2.5,5,10]$ and see how this impacts our results and research goals.

We can similarly use the Ridge() function to fit a linear ridge regression model. This Ridge() function operates in a similar way to the Lasso() function.

Fitting Four Linear Ridge Regression Models

Model 1: $\lambda=1$

from sklearn.linear_model import Ridge
ridge_mod_1 = Ridge(alpha=1, max_iter=1000)
ridge_mod_1.fit(X_train, y_train)

Ridge(alpha=1, max_iter=1000)

Ridge(alpha=1, max_iter=1000)

Model 2: $\lambda=2.5$

ridge_mod_25 = Ridge(alpha=2.5, max_iter=1000)
ridge_mod_25.fit(X_train, y_train)

Ridge(alpha=2.5, max_iter=1000)

Ridge(alpha=2.5, max_iter=1000)

Model 3: $\lambda=5$

ridge_mod_5 = Ridge(alpha=5, max_iter=1000)
ridge_mod_5.fit(X_train, y_train)

Ridge(alpha=5, max_iter=1000)

Ridge(alpha=5, max_iter=1000)

Model 4: $\lambda=10$

ridge_mod_10 = Ridge(alpha=10, max_iter=1000)
ridge_mod_10.fit(X_train, y_train)

Ridge(alpha=10, max_iter=1000)

Ridge(alpha=10, max_iter=1000)

Inspecting the Slopes

Next, let's inspect and compare the slopes of the models that we've created so far.

df_slopes = pd.DataFrame({'lin_reg_mod': lin_reg_mod.coef_.T, 
                         'ridge_mod_1': ridge_mod_1.coef_.T,
                         'ridge_mod_25': ridge_mod_25.coef_.T,
                         'ridge_mod_5': ridge_mod_5.coef_.T,
                         'ridge_mod_10': ridge_mod_10.coef_.T}, index=X_train.columns)
df_slopes

By comparing the slopes of the four linear ridge regression models, we can see what we would expect. As we increase $\lambda$, the squared sum of the slopes decreases more and more.

(df_slopes**2).sum()

lin_reg_mod 1989.465059
ridge_mod_1 454.498812
ridge_mod_25 226.579597
ridge_mod_5 120.894160
ridge_mod_10 61.376527
dtype: float64

We can see that many slopes individually decreased as we increased $\lambda$, however there were some that increased, like X1109.

Also notice, unfortunately, none of our slopes were zeroed out as in the case with the LASSO model. Thus, the task of choosing which genes may be overfitting the model and thus we should delete is much more difficult.

plt.figure(figsize=(14,5))
plt.plot(df_slopes['lin_reg_mod'], label='lin_reg_mod')
plt.plot(df_slopes['ridge_mod_1'], label='ridge_mod_1')
plt.plot(df_slopes['ridge_mod_25'], label='ridge_mod_25')
plt.plot(df_slopes['ridge_mod_5'], label='ridge_mod_5')
plt.plot(df_slopes['ridge_mod_10'], label='ridge_mod_10')
plt.axhline(0, color='black', linestyle='--')
plt.ylabel('Slope')
plt.xlabel('Gene Explanatory Variable')
plt.xticks(rotation=90)
plt.legend(bbox_to_anchor=(1,1))
plt.show()

We might further inspect the ridge regression slopes by dividing each slope by it's original linear regression model slope to measure the factor by which it changed.

df_slopes_normalized = df_slopes.div(df_slopes['lin_reg_mod'], axis=0)
df_slopes_normalized

We see that as the value of $\lambda$ is increased, there are some slopes (like X980) that dramatically decline in magnitude compared to their original magnitude in the linear regression model. On the other hand, there were other slopes (like X1109) that dramatically increased in magnitude compared to their original magnitude in the linear regression model.

plt.figure(figsize=(14,5))
plt.plot(df_slopes_normalized['ridge_mod_1'], label='ridge_mod_1')
plt.plot(df_slopes_normalized['ridge_mod_25'], label='ridge_mod_25')
plt.plot(df_slopes_normalized['ridge_mod_5'], label='ridge_mod_5')
plt.plot(df_slopes_normalized['ridge_mod_10'], label='ridge_mod_10')
plt.axhline(0, color='black', linestyle='--')
plt.ylabel('Percent of Linear Regression Slope')
plt.xlabel('Gene Explanatory Variable')
plt.xticks(rotation=90)
plt.ylim([-1,2])
plt.legend(bbox_to_anchor=(1,1))
plt.show()

Which model is best?

Finally, let's determine which ridge regression model is best for the purpose of predicting the size of new breast tumors, by again comparing the test data R^2 of each model.

Given that one of our goals is to build a model that will yield good predictions for new breast tumors, let's inspect the test dataset R^2 of each of these 3 LASSO models.

Nonregularized Linear Regression Model: Test Dataset R^2

lin_reg_mod.score(X_test, y_test)

-0.3071595344334954

Ridge Regression Model 1: $\lambda=1$ Test Dataset R^2

ridge_mod_1.score(X_test, y_test)

0.12825147516877422

Ridge Regression Model 2: $\lambda=2.5$ Test Dataset R^2

ridge_mod_25.score(X_test, y_test)

0.24810734542753

Ridge Regression Model 3: $\lambda=5$ Test Dataset R^2

ridge_mod_5.score(X_test, y_test)

0.27410328541870965

Ridge Regression Model 4: $\lambda=10$ Test Dataset R^2

ridge_mod_10.score(X_test, y_test)

0.2607143483417401

We can see that the linear ridge regression model that achieves the best test dataset R^2=0.274 is the model that used $\lambda=5$. Ideally, a more complete analysis would try out many other values of $\lambda$ to see if we can find a test dataset R^2 that was even higher. Notice, that this ridge regression model was not quite as high as the best LASSO model that we found with a test dataset R^2=0.279.

Can we trust the non-zero slope interpretations in our chosen ridge regression model?

Ideally, we'd be able to trust our interpretations of our non-zero slopes in this resulting ridge regression model that we chose as well. Because ridge regression is more effective at downplaying the impacts of multicollinearity (which we know this dataset has), we may be able to trust the interpretive value of the slopes in our ridge regression more.

Notice how the magnitude of some slopes dramatically changed in the two models below.

df_slopes = pd.DataFrame({'lin_reg_mod': lin_reg_mod.coef_.T, 
                         'ridge_mod_5': ridge_mod_5.coef_.T}, index=X_train.columns)
df_slopes

Elastic Net Regression (ie. L1 and L2 Penalty Term)

We can see that this is a combination of the:

Similar to LASSO and ridge regression, we introduce a parameter $\lambda\geq 0$ that we must preselect. We also now use an $0\leq \alpha\leq 1$ parameter that we must also preselect.

The Ideal Scenario

Similar to the LASSO and ridge linear regression, the absolute best intercept and slope solutions to this optimization would be ones in which:

The Penalty Term Balancing Act

However, similarly, this "ideal scenario" that we just described is unrealistic. We also see that choosing intercepts and slopes that lead to a decrease of one term, often leads to an increase in the other term, and vice versa.

The Impact of $\lambda$

For similar reasons described when we used $\lambda$ in the LASSO linear regression, we see the following impact of $\lambda$ on our results.

$\lambda$ is Set Low

$\lambda$ is Set High

Which slope magnitudes get to stay large?

Similar to LASSO and ridge regression, the resulting optimal slopes that are set "low", are the ones that bring the least predictive power to the model and thus are the most likely variable to overfit the model.

The Impact of the $\alpha$ Parameter

Now, let's explore the impact that $\alpha$ has on our resulting slopes.

• $\alpha$ is Set High

  When $\alpha$ is set higher, the penalty coefficient in front of the L1 term $\left[|\hat{\beta}_1|+|\hat{\beta}_2|+...+|\hat{\beta}_{50}|\right]$ becomes higher. Thus, the optimization problem is going to try harder to keep this term in particular low. Thus, your results slopes will tend to look more like those that were returned by LASSO regression (ie. zeroing slopes out).

• $\alpha$ is Set High

  Alternativly, When $\alpha$ is set lower, the penalty coefficient in front of the L2 term $\left[\hat{\beta}_1^2+\hat{\beta}_2^2+...+\hat{\beta}_{50}^2\right]$ becomes higher. Thus, the optimization problem is going to try harder to keep this term in particular low. Thus, your results slopes will tend to look more like those that were returned by ridge regression (ie. less zeroing slopes out and more focus on balancing the weight in collinear slopes).

Less Clear Slope Interpretation Depending on $\alpha$

The smaller the value of $\alpha$ is set, ie. the more the returned slopes tend to resemble a ridge regression solution, the less the small slopes will be set exactly equal to 0. Thus, the less the slopes have interpretive value about which variables should be removed from the model because they do not bring enough predictive power.

Benefits of Elastic Net Regression

The benefits of elastic net regression is that it essentially can provide a "middle ground" solution between the type of results that you might get with LASSO and ridge regression. By toggling the $\alpha$ parameter, you can have your results resemble the type of solutions that you might get with LASSO vs. ridge regression to varying degrees.

Elastic Net Linear Regression for Predicting Tumor Size

Finally, let's actually apply our elastic netlinear regression model to this dataset.

Recall that our linear ridge regression objective function required for us to preselect a value for $\lambda$. Let's try out a series of different $\lambda$ values to see which model yields the best test dataset performance.

Also recall, that we need to preselect our $\alpha$ value as well. A more complete analysis would also try out many different combinations of $\alpha$ and $\lambda$ values to assess which model gave us the best test data results. However, for now we will just try $\alpha=0.7$. Because this $\alpha$ is higher and closer to 1, we might expect our resulting slopes to resemble more what we might get with a LASSO model.

We can similarly use the ElasticNet() function to fit a linear ridge regression model. This ElasticNet() function operates in a similar way to the Lasso() function. Let's be careful about the parameters here.

• the alpha parameter in the ElasticNet() function below represents the $\lambda$ parameter in the objective function above
• the l1_ratio parameter in the ElasticNet() function below represents the $\alpha$ parameter in the objective function above.

Fitting Four Elastic Net Linear Regression Models

Model 1: $\lambda=0.01$

from sklearn.linear_model import ElasticNet
en_mod_1 = ElasticNet(alpha=.01, l1_ratio=0.7)
en_mod_1.fit(X_train, y_train)
en_mod_1.score(X_test, y_test)

0.06547371156125181

Model 2: $\lambda=0.025$

en_mod_2 = ElasticNet(alpha=.025, l1_ratio=0.7)
en_mod_2.fit(X_train, y_train)
en_mod_2.score(X_test, y_test)

0.23638055648260659

Model 3: $\lambda=0.05$

en_mod_3 = ElasticNet(alpha=.05, l1_ratio=0.7)
en_mod_3.fit(X_train, y_train)
en_mod_3.score(X_test, y_test)

0.28513850610195246

Model 4: $\lambda=0.1$

en_mod_4 = ElasticNet(alpha=.1, l1_ratio=0.7)
en_mod_4.fit(X_train, y_train)
en_mod_4.score(X_test, y_test)

0.2639254089412595

We can see that the elastic net model that used $\lambda=0.05$ and $\alpha=0.7$ achieved the highest test dataset R^2=0.285. Also note that this is the highest breast tumor test dataset R^2 that we have achieved so far out of all explored models!

Inspecting the Slopes

Next, let's inspect and compare the slopes of the models that we've created so far. Notice how we do see some zeroed out slopes in our results, regardless of what $\lambda$ value that we chose. This makes sense as our $\alpha$ was set larger, and therefore it's results will look more similar to what we would get with a LASSO model.

Also take note, just like we observed with the LASSO models, the higher the value of $\lambda$, the more slopes that were zeroed out.

df_slopes = pd.DataFrame({'lin_reg_mod': lin_reg_mod.coef_.T, 
                         'en_mod_1': en_mod_1.coef_.T,
                         'en_mod_2': en_mod_2.coef_.T,
                         'en_mod_3': en_mod_3.coef_.T,
                         'en_mod_4': en_mod_4.coef_.T}, index=X_train.columns)
df_slopes

Can we trust the non-zero slope interpretations in our chosen elastic net model?

Ideally, we'd be able to trust our interpretations of our non-zero slopes in this resulting elastic net model that we chose which also happened to have the best test R^2.

Because we set $\alpha$ to be high, this means that our results more closely resembled a LASSO model solution. Thus, we'd also want to inspect whether the remaining explanatory variables (ie. those with non-zero slopes) were collinear. Unfortunately, we do see that quite a few of these remaining explanatory variables are still collinear. Thus we should be more cautious about interpreting the slopes in this model.

nonzero_slope_genes = df_slopes[np.abs(df_slopes['en_mod_3'])>0.00001].index
X_train_remaining = X_train[nonzero_slope_genes]

sns.heatmap(X_train_remaining.corr(),vmin=-1, vmax=1, cmap='RdBu')
plt.show()

On the other hand if we had set our $\alpha$ to be lower, that would mean that our results would closely resemble a ridge regression model solution, which inherently tends to mitigate the selection of slopes with unreliable interpretations due to multicollinearity.

Conclusion and Potential Shortcomings

Goal Assessment

Goal 1: Best Predictions for New Breast Tumors

Thus, if our first and foremost goal was to select a model which we think would perform the best predicting new breast tumor sizes, then we would select our model that achieved the highest test data R^2=0.285. This was our elastic net model with $\lambda=0.05$ and $\alpha=0.7$. However, as we mentioned above we cannot quite trust our slopes interpretations due the extent to which the elastic net model did not assist in dealing with the multicollinearity.

Goal 2: More Trust in Slope Interpretations

However, if our goal was to achieve a model in which we'd be able to trust the interpretations of our slopes more (in the unfortunate presence of highly multicollinar explanatory variables), then we might consider selecting one of our ridge regression models. In addition, the test data R^2=0.274 of our best ridge regression model was not too much lower than that of our best final model.

Analysis Shortcomings

In our analysis above, we created a single training and test dataset by randomly splitting the dataset using a random state of 102.

Let's perform this same analysis again but using a different random state = 207 when making this split.

#Creating a new set of training and test datasets
df_train_other, df_test_other = train_test_split(df, test_size=0.1, random_state=207)
X_train_other = df_train_other.drop(['size'], axis=1)
X_test_other = df_test_other.drop(['size'], axis=1)
y_train_other=df_train_other['size']
y_test_other=df_test_other['size']

#Fitting the linear regression model
lin_reg_mod = LinearRegression() 
lin_reg_mod.fit(X_train_other, y_train_other)
print('Linear Regression Test R^2:',lin_reg_mod.score(X_test_other, y_test_other))

Linear Regression Test R^2: -0.2213718065663819

Unfortuately, we see that by just choosing a different random state to split up the training and test dataset our analysis shows us different results.

1. All of our elastic net model test data R^2 values are negative, indicating that all of our elastic net trained models were poor fits for the test dataset.
2. The best elastic net model in this case the one with $\lambda=.1, \alpha=0.7$, because it had the highest test data R^2. Notice how this is a different $\lambda$ parameter than what was chosen in our first analysis.

#Fitting 4 elastnic net models with the same parameters as before
en_mod = ElasticNet(alpha=.01, l1_ratio=0.7)
en_mod.fit(X_train_other, y_train_other)
print('Elastic Net - alpha=0.01 Test R^2:',en_mod.score(X_test_other, y_test_other))

en_mod = ElasticNet(alpha=.025, l1_ratio=0.7)
en_mod.fit(X_train_other, y_train_other)
print('Elastic Net - alpha=0.025 Test R^2:',en_mod.score(X_test_other, y_test_other))

en_mod = ElasticNet(alpha=.05, l1_ratio=0.7)
en_mod.fit(X_train_other, y_train_other)
print('Elastic Net - alpha=0.05 Test R^2:',en_mod.score(X_test_other, y_test_other))

en_mod = ElasticNet(alpha=.1, l1_ratio=0.7)
en_mod.fit(X_train_other, y_train_other)
print('Elastic Net - alpha=0.1 Test R^2:',en_mod.score(X_test_other, y_test_other))

Elastic Net - alpha=0.01 Test R^2: -0.14750397218825695
Elastic Net - alpha=0.025 Test R^2: -0.11125207416411964
Elastic Net - alpha=0.05 Test R^2: -0.0856548712042815
Elastic Net - alpha=0.1 Test R^2: -0.04788700566442472

This highlights the fact that the act of selecting a single training and test dataset based on a single random state and basing your model selection decisions off can have consequences. Doing so introduces an element of randomness to your model selections. Notice how we get two sets of different slopes, particularly in which different sets of slopes are set to 0.

Based on the differing results that we observed above, you should be more skeptical that your final selected elastic net model from 8.4 will indeed be the best when it comes to predicting new breast tumor sizes.

df_slopes = pd.DataFrame({'best_en_mod_old_rs': en_mod_3.coef_.T,
                         'best_en_mod_new_rs': en_mod.coef_.T}, index=X_train.columns)
df_slopes