course: (ML Series) Regression Approaches (Linear, Ridge, Lasso)
By Julien Hernandez Lallement, 2020-10-03, in category Course
By Julien Hernandez Lallement, 2020-10-03, in category Course
Generally speaking, regression methods are supervised learning techniques. They use continous scaled variables (independent variables) to predict the behavior of a dependent variable. They can use different equations that will fit straight lines (linear regression), polynomial functions (detecting interaction effects) or other functions to predict the dependent variable.
In this post, I will be focusing on linear regressions in low and high dimensions. In a subsequent post, I will discuss logistic regression which is quite useful for classification problems (not a classification algorithm per se though, only if a probability threshold is set).
Note that part of this post is a copy of the introduction to machine learning post, which used linear regression as an illustrative example
Regressions can be used in different scenarios:
Suppose that you obtain data for age and height of a bunch of people. You would be able to look at the linear relationship between the two. For a given person i, we can represent the age by X and the height by Y:
show_example()
We therefore want to find a relationship of the form:
where $w$ are the parameters that we want to play with in order to obtain the best fit of the data.
You can then more variables to the model, which would become a multi-variate model:
$$ y(\mathbf{w}, \mathbf{x}) = w_0 + w_1 x_1 (+ w_2 x_2 + ...) $$In some cases, a straight line will not have the best fit (this is something that you can test of course, but here a reasonable amount of common sense and more importantly data visualization will allow you to determine which model could fit best). In such cases, polynomial functions might allow you to get a better fit (see equation above).
Now, using this linear function, we want to measure how far away the actual data y is from the model. Question is, how do we define this line? How do we maximize the parameters so that the line captures the most variance in the data, maximizing the information we extract from each data point?
One could think of drawing many different lines, each better capturing a particular sub-population of the plotted data...
show_example(random_lines=3)
One way to capture a best fit for linear models is to use the residuals, which capture a measure of error between the actual data, and the plotted line.
show_example(plot_errors=True)
The black dotted lines represent the difference between what your model would predict, and the real values of the target you are aiming at predicting. Reducing these errors is what you should aim for, without falling into the issue of over-fitting which I discuss a bit later. In other words, we want to minimize the residuals.
For example, let's take the Mean Squared Error (MSE; AKA Residual Sum of Squares) cost function:
$$ RSS : J(w) = \frac{1}{n}\sum \limits _{i=1} ^{n} [y_i - w_0 - \sum \limits _{j=1} ^{p} (w_j - x_ij)]^2 $$Which you can write using the unwrapped equivalent of the $y$ function used above:
where:
n = number of training examples
y(w,x^i) = observed value
$yi$ = predictions
$w$ = parameter coefficient
This is basically a simple sum of the squares of the vertical distances between the real data points and the straight line of your model.
The constant $\frac{1}{n}$ in front plays an important role for high sample sizes.
The MSE quantifies the average deviation of predicted values from the actual correct values. Conceptually, I can be explained in terms of model’s performance on the training set (the dataset fed to the data before assessing its generalization to unseen data). As you can see from the equation, the cost is larger when the model is performing poorly on this data set.
It is a popular choice given that:
One objective of the algorithm you would be using is to minimize the cost (term J in the equation above), by testing different parameters.
Since you are looking for the lowest amount of error, an easy way to find the minimum analytically:
Which allows you to obtain both parameters $\theta_0$ & $\theta_1$
Sometimes, your data will be correlated, but the overall shape won't folow one of a straigt line. In such cases, one can fit a a polynomial regression, which will result in findind the coefficients of a polynomial of a specific degree d that best fits the data. As such, linear regressions are a type of polynomial regression where d=1.
Note that this is not bounded to Linear Regression.
Note as well that the higher the number of degrees, the higher is the danger of overfitting.
Below I will be using PolynomialFeatures from Sklearn as demonstration.
I will be using two different famous use cases: the ChickWeight & the Boston Housing datasets.
Below, I will use the regression approach on the famous Chick Weight dataset.
import pandas as pd
import numpy as np
import plotnine as p9
import os
os.chdir('/home/julien/website/datasets')
df = pd.read_csv('ChickWeight.csv')
df.head(10)
This dataset contains information about single chicks ("Chick" column), namely:
# Low case all column names
df.columns= [c.lower().replace(" ","") for c in df.columns]
def make_ml_data(dataf):
return (dataf
.dropna() # Just in case, drop NaNs
.assign(weight_float = lambda d: d['weight'].astype(float))
.assign(time_float = lambda d: d['time'].astype(float))
.assign(diet_str = lambda d: d['diet'].astype(str))
[['weight_float', 'time_float', 'diet_str']] )
ml_df = df.pipe(make_ml_data)
ml_df.head(10)
Above, I simply changed the format to be able to work with the data as numbers, and stored it in a new dataframe called ml_df.
Let's have a quick look at the data, per diet
(p9.ggplot()
+ p9.geom_point(data=df, mapping=p9.aes("time", "weight", color="chick"))
+ p9.facet_wrap('diet')
)
There seems to be indeed a rather linear relation in each diet between time and weight.
Now let's see if we can capture that using a regression approach.
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler, PolynomialFeatures, OneHotEncoder
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.model_selection import GridSearchCV
from sklego.preprocessing import ColumnSelector
I first will only use the weight to predict the age of the chick (time).
# Preprocessing pipeline
panda_grabber = Pipeline([
("union", FeatureUnion([
("continous", Pipeline([
("select", ColumnSelector(["weight_float"])),
("scale", StandardScaler())
]))
]))
])
# Fit and apply pipeline
panda_grabber.fit(ml_df,ml_df['diet_str']).transform(ml_df);
I define an ML pipeline.
# ML pipeline
pipeline = Pipeline([
("grab", panda_grabber),
("poly", PolynomialFeatures(interaction_only=True)),
("model", LinearRegression())
])
# Define parameters
param_poly = [PolynomialFeatures(interaction_only=True),
None]
# Define model to be used
param_model = [LinearRegression()]
# Define Grid Search approach
mod = GridSearchCV(pipeline,
#iid=True,
return_train_score=True,
param_grid={"model": param_model,
"poly": param_poly},
cv=20)
# Fit model
mod.fit(ml_df, ml_df['time_float']);
mod_results = pd.DataFrame(mod.cv_results_)
mod_results[["mean_train_score","mean_test_score"]]
Above, the results of the model fit.
We have two lines since we provided two parameters to the model: polynomial fit, and no additional parameters. In this case, since we are using only one parameter, there is no difference between the scores of both models.
We have a good accuracy (~70%)in predicting the age of the chick based on the weight provided.
We also see that the accuracy is comparable between training and testing set, which suggest that we are not suffering from under or over fitting.
Now let's try to add one additional feature to the model, the diet.
# Preprocessing pipeline
panda_grabber = Pipeline([
("union", FeatureUnion([
("continous", Pipeline([
("select", ColumnSelector(["weight_float"])),
("scale", StandardScaler())
])),
("coutinous", Pipeline([
("select", ColumnSelector(["diet_str"])),
("encode", OneHotEncoder(sparse=False))
]))
]))
])
# Fit and apply pipeline
panda_grabber.fit(ml_df,ml_df['time_float']).transform(ml_df);
# ML pipeline
pipeline = Pipeline([
("grab", panda_grabber),
("poly", PolynomialFeatures(interaction_only=True)),
("model", LinearRegression())
])
# Define parameters
param_poly = [PolynomialFeatures(interaction_only=True),
None]
# Define model to be used
param_model = [LinearRegression()]
# Define Grid Search approach
mod = GridSearchCV(pipeline,
#iid=True,
return_train_score=True,
param_grid={"model": param_model,
"poly": param_poly},
cv=20)
# Fit model
mod.fit(ml_df, ml_df['time_float']);
mod_results = pd.DataFrame(mod.cv_results_)
mod_results
We can see that now that we have added a novel feature to the model, the polynomial fit is different than the naive linear fit.
We can also see that the accuracy is slightly better for the polynomial fit.
mod_results[["mean_train_score","mean_test_score"]]
Adding the diet increases our accuracy by ~4%, which is not bad ;)
I am creating 10000 data points taken from a uniform distribution.
Since we are trying to predict the age of the chick, I am drawing data for weight, and for demonstration purpose, I will simply use one diet type
n = 10000
simulated_df = pd.DataFrame({
"weight_float": np.random.uniform(ml_df.weight_float.min(), ml_df.weight_float.max(), n),
#"time_float": np.random.uniform(ml_df.time_float.min(), ml_df.time_float.max()+5, n),
"diet_str": "1"
})
Now I can apply the model to this data to see what the predictions would be
pltr = simulated_df.assign(pred = lambda d: mod.predict(d))
Finally, let's plot the real data and compare the age given for diet = 1 and different weight
# Assign some noise to the x_axis for better visualization
pltr = pltr.assign(x_axis = lambda d: d.diet_str.astype(float) + 0.5 + np.random.rand(n,1).flatten())
(p9.ggplot()
+ p9.geom_point(data=pltr, mapping=p9.aes('x_axis', 'weight_float', color='pred'))
+ p9.geom_point(data=df.loc[df.diet_str == '1'], mapping=p9.aes('diet', 'weight', color='time'), size=5, alpha=0.5)
)
Visually, we can see a common pattern between the real data one the left, and what the model predicted for the simulated data points on the right.
In cases where the Linear Regression is over or under fitting the data, one can play with the cost function to control the way the coefficients are computed
Let's first look again at the cost function we defined previouly, the Residual Sum of Squares:
As discussed above, using this cost dunction will adjust the coefficients to fit the data best. Importantly, the fit will be done on the training dataset. Now, if there is noise in the training data, or it is of poor quality, then the estimated coefficients will cause the model to not generalize well to the future data (overfitting). This is where regularization can be useful since it will shrink or moduldate the learned estimates towards zero. In other wrods, regularization reduces the variance of the model, without substantial increasing its bias.
Ridge and Lasso regressions are variant of a naive linear model, that can help when we are suffering from over or under fitting
In a ridge regression, the cost function is modified by adding a penalty that will control for for multicolinearity. That penalty is equivalent to the square of the magnitude of the coefficients.
Adding the regularization term forces the model to select simpler models to fit the data and since we omit the $w_0$ term in the regularization term, the larger the $w_j$, the greater will we see a change in the prediction.
$\lambda$ is the tuning parameter, that determines the level of penalty of our model. This is how the Ridge regression technique prevents coefficients from rising too high.
When $\lambda$ = 0, no penalty is applied, and the estimates will be equal to least squares (as in naive linear regression). However, as $\lambda$ increases, the penalty grows, and the ridge regression coefficient estimates will approach 0. FYI, this is also known as the L2 norm.
Hence, the parameter $\lambda$ is critical for a good model fit (cross validation can be useful here).
For one dimensional linear regression, regularization terms might not be extremely useful, but they become interesting for higher order polynomials, where the different features might influence the model in an uncontrolled manner. In that case, the regularization term will reduce the size of some coefficients, controlling for unwanted effect on the fit.
Below, I show the effect of the regularization on the coefficients. I will be using the boston housing dataset.
import os
os.chdir('/home/julien/website/datasets')
import matplotlib.pyplot as plt
df = pd.read_csv('BostonHousing.csv')
from sklearn.linear_model import Ridge
from sklearn.model_selection import train_test_split
Below, a function that will fit different levels of lambda using a Ridge Regression, and plot the coefficients for each feature
def run_regularizations(data, alpha, model):
# Split data
X_train,X_test,y_train,y_test=train_test_split(data.drop(columns={'medv'}),
data.medv,
test_size=0.3,
random_state=3
)
# Plot parameters
marker_size=10
alpha=0.5
colors = ['red', 'blue', 'green', 'cyan', 'black']
markers = ['*', 'd','o', 'h', 's', 'p']
pos = -1
# Fit models
for lambdas in alpha_levels:
pos = pos + 1
mod = model(alpha=lambdas)
mod.fit(X_train, y_train)
#print('**'+str(model)+(' ')+'alpha = '+str(lambdas)+'**')
print('lambda level = ' +str(lambdas))
#print('')
print("Train Score: ", mod.score(X_train, y_train))
print("Test Score: ", mod.score(X_test, y_test))
print('')
print('Number of feature used: ', np.sum(mod.coef_!=0))
print('');print('')
# Plot
plt.plot(mod.coef_,
alpha=alpha,
linestyle='none',
marker=markers[pos],
markersize=marker_size,
color=colors[pos],
label=[str(model) + str(lambdas)]
)
plt.xlabel('Feature',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.xticks(ticks = range(0,13),labels=data.columns[:-1],rotation=90)
plt.legend(fontsize=13,
bbox_to_anchor=(1.05, 1),
loc='upper left',
)
plt.show()
I will test two different values of $\lambda$, which in Sklearn is referred to as alpha.
alpha_levels=[0.01,2,100]
run_regularizations(data=df, alpha=alpha_levels, model=Ridge)
Here, we can see how the coefficient level approaches 0 with increasing values of $\lambda$. Note that this plot is for demonstration purpose, since a naive Linear Regression already provides a good fit, and non regularization is needed on this dataset.
It should be noted that the Ridge Regression, due to the use of the L2 norm, will shrink the coefficients for least important predictors. However, and importantly, these coefficients will never be equal to 0, and all predictors will be kept in the model.
However, the L1 norm used in the Lasso Regression will force some coefficient to be equal to 0 when the parameter $\lambda$ is sufficiently large.
Therefore, the Lasso Regression can also performs feature selection, which is a advantage over the Ridge Regression. This can be quite useful if one wants to reduce the number of features used by the model.
Mathematically, this aspect of Lasso Regressio is enabled by the shape of the cost function, which will cause some feature to vanish when the coefficient takes a certain value relative to the cost function of a linear regression. For more information, check this book.
Lasso Regression (Least Absolute Shrinkage and Selection Operator) is another variation of naive linear regressions. Instead of squaring the coefficients, we take the absolute value, as shown below:
This variation differs from a ridge regression in that it will penalize mainly high coefficients. FYI, this is known as the L1 norm.
Note that simimlarly as for the Ridge Regression, the $\lambda$ parameter is crucial in the Lasso Regression.
Let s be any value that can be taken by $\lambda$.
Altogether, the Ridge Regression can be thought of as solving an equation, where summation of squares of coefficients is less than or equal to s, while the Lasso Regression can be thought of as an equation where summation of modulus of coefficients is less than or equal to s. See here for quote and more information on Ridge and Lasso.
from sklearn.linear_model import Lasso
Again, I will test two different values of $\lambda$.
In addition to the same code used for the Ridge Regression, I will be looking at the number of features selected by the Lasso Regression, depending on the value of $\lambda$.
In addition, I will be using a simple Linear Regression, as comparison.
alpha_levels=[0.1,2,5,30,100]
run_regularizations(data=df, alpha=alpha_levels, model=Lasso)
Note how the coefficients are able to take a 0 value, hence resulting in a decreasing number of features used for an increasing value of $\lambda$.
This is selection feature performed by a Lasso Regression.
As for the Ridge Regression above, this plot is for demonstration purposes. I recommend using this algorithm on dataset with a higher number of features, since it might help to visualize the power of Lasso's feature selection inherent functions. He avoid to do here to no include too many different datasets.
Since the Lasso and Ridge regression do not seem to add value, I will propose a pipeline using only Linear Regression.
corr = df.corr().abs()
corr.style.background_gradient(cmap='coolwarm')
Based on the correlation coefficient, I should probably be working with lstat, tax, rm & ptratio in the model.
But for this demonstration, I will use all features provided in the dataset.
X = df.drop(columns={'medv'})
y = df.medv
# Scale data
scaler = StandardScaler()
scaler.fit(X)
scaled_X = scaler.transform(X)
#Split data
X_train,X_test,y_train,y_test=train_test_split(scaled_X,
y,
test_size=0.3,
random_state=3
)
# Load model
lr = LinearRegression()
# Fit model
lr.fit(X_train,y_train)
# Show predictions vs true values
plt.scatter(x=(lr.predict(X_test)),y=y_test);
plt.xlabel('Predictions');
plt.ylabel('True values');
Let's look at each feature and how it influences the predicted value
lr.fit(scaled_X, y)
sorted(list(zip(X.columns, lr.coef_)), key=lambda coeff: abs(coeff[1]), reverse=True)
The coefficients above represent the average change in the housing price for one unit of change in the feature while holding other features in the model constant. For example, holding all else constant, a unit increase in LSTAT decreases our target (housing price) by 3.469 and a unit increase in RM increases our target by 2.67.
Be careful though! We are looking at correlative data, so while looking at variables such as rm and lstat might be a good place to start, we cannot conclude that decreasing and increase rm and lstat, respectively, affects housing price.
You will soon realize how important it will be to obtain confidence interval around your estimates. One way to do that is to use a bootstrap. I have already described this method in other posts related to permutation.
A bootstrap is a random sample with replacement of our data and this sample is of the same size of the original data. This is a way of generating multiple views of the same data. Let's create 1000 bootstraps of our data. This example is taken from here.
from sklearn.utils import resample
n_bootstraps = 1000
bootstrap_X = []
bootstrap_y = []
for _ in range(n_bootstraps):
sample_X, sample_y = resample(scaled_X, y)
bootstrap_X.append(sample_X)
bootstrap_y.append(sample_y)
Then on each of these datasets, we can fit our model and get our coefficents:
linear_regression_model = LinearRegression()
coeffs = []
for i, data in enumerate(bootstrap_X):
linear_regression_model.fit(data, bootstrap_y[i])
coeffs.append(linear_regression_model.coef_)
coef_df = pd.DataFrame(coeffs, columns=X.columns)
plt.style.use("classic")
coef_df.plot(kind='box')
plt.xticks(rotation=90);
coef_df['lstat'].plot(kind='hist');
We can now say with confidence that the LSTAT coefficient is negative, and almost certaintly located between ~-1.2 and ~-5.5.
That was it for the linear regression. I hope you enjoyed it, don't hesitate in contacting me if you see any mistakes or have any question :)
Cheers!