Linear regression is a supervised learning approach for predicting a quantitative output variable from one or several input variables. It is one of the most widely used techniques in statistics, and in machine learning it is one of the simplest approaches for predicting a quantitative outcome. It's simplicity also allows for easy interpretation of the importance of the input variables.
The predicted outputs, $\mathbf{\hat{y}}$, are created using a linear combination of the input variables:
$\mathbf{\hat{y}} = \beta_{0} + \beta_{1}\mathbf{x_{1}} + \dots + \beta_{P}\mathbf{x_{P}}$.
Our goal is to find the regression coefficients, $\beta_{0}, \beta_{1}, ..., \beta_{P}$, that minimize how far our predictions are from the known outputs $\mathbf{y}$. By far the most common criteria to minimize is the sum of squared residuals (each residual is the difference between the predicted outcome and the known outcome).
$\sum_{n = 1}^{N}(y_{n} - \hat{y}_{n})^2$.
Without further ado, let's look at linear regression with scikit-learn
(with a special guest appearance from the statsmodels
package).
scikit-learn
User GuideWe'll get started by importing the items we need from our various libraries.
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split, RandomizedSearchCV
from sklearn.utils.fixes import loguniform
from sklearn.metrics import mean_squared_error, r2_score
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
plt.rcParams['figure.figsize'] = [10, 6]
This example comes from ISLR in their Chapter 3 Lab.
The Boston Data contains information on census tracts in Boston. Our goal is to predict medv
, the median house value, using a variety of quantitative predictors about the census tracts.
These data are available from several sources, including the ISLR2
R package, whose datasets I recently exported to .csv files and uploaded to a GitHub repo.
data_df = pd.read_csv("https://raw.githubusercontent.com/LukeMoraglia/ISLR_datasets/master/data/Boston.csv")
data_df.head()
crim | zn | indus | chas | nox | rm | age | dis | rad | tax | ptratio | lstat | medv | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0.00632 | 18.0 | 2.31 | 0 | 0.538 | 6.575 | 65.2 | 4.0900 | 1 | 296 | 15.3 | 4.98 | 24.0 |
1 | 0.02731 | 0.0 | 7.07 | 0 | 0.469 | 6.421 | 78.9 | 4.9671 | 2 | 242 | 17.8 | 9.14 | 21.6 |
2 | 0.02729 | 0.0 | 7.07 | 0 | 0.469 | 7.185 | 61.1 | 4.9671 | 2 | 242 | 17.8 | 4.03 | 34.7 |
3 | 0.03237 | 0.0 | 2.18 | 0 | 0.458 | 6.998 | 45.8 | 6.0622 | 3 | 222 | 18.7 | 2.94 | 33.4 |
4 | 0.06905 | 0.0 | 2.18 | 0 | 0.458 | 7.147 | 54.2 | 6.0622 | 3 | 222 | 18.7 | 5.33 | 36.2 |
data_df.describe()
crim | zn | indus | chas | nox | rm | age | dis | rad | tax | ptratio | lstat | medv | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
count | 506.000000 | 506.000000 | 506.000000 | 506.000000 | 506.000000 | 506.000000 | 506.000000 | 506.000000 | 506.000000 | 506.000000 | 506.000000 | 506.000000 | 506.000000 |
mean | 3.613524 | 11.363636 | 11.136779 | 0.069170 | 0.554695 | 6.284634 | 68.574901 | 3.795043 | 9.549407 | 408.237154 | 18.455534 | 12.653063 | 22.532806 |
std | 8.601545 | 23.322453 | 6.860353 | 0.253994 | 0.115878 | 0.702617 | 28.148861 | 2.105710 | 8.707259 | 168.537116 | 2.164946 | 7.141062 | 9.197104 |
min | 0.006320 | 0.000000 | 0.460000 | 0.000000 | 0.385000 | 3.561000 | 2.900000 | 1.129600 | 1.000000 | 187.000000 | 12.600000 | 1.730000 | 5.000000 |
25% | 0.082045 | 0.000000 | 5.190000 | 0.000000 | 0.449000 | 5.885500 | 45.025000 | 2.100175 | 4.000000 | 279.000000 | 17.400000 | 6.950000 | 17.025000 |
50% | 0.256510 | 0.000000 | 9.690000 | 0.000000 | 0.538000 | 6.208500 | 77.500000 | 3.207450 | 5.000000 | 330.000000 | 19.050000 | 11.360000 | 21.200000 |
75% | 3.677083 | 12.500000 | 18.100000 | 0.000000 | 0.624000 | 6.623500 | 94.075000 | 5.188425 | 24.000000 | 666.000000 | 20.200000 | 16.955000 | 25.000000 |
max | 88.976200 | 100.000000 | 27.740000 | 1.000000 | 0.871000 | 8.780000 | 100.000000 | 12.126500 | 24.000000 | 711.000000 | 22.000000 | 37.970000 | 50.000000 |
Below is a pair plot of all the numeric variables. Our target, medv
, is the bottom row and rightmost column. Based on these scatterplots, we might expect rm
, the average number of rooms per dwelling, and lstat
, the percentage of the population with "lower status" (not really sure what that means), to be good linear predictors of medv
.
sns.pairplot(data_df)
plt.show()
Next, we'll split our data into the predictors X
and the target y
.
y = data_df["medv"]
X = data_df.drop("medv", axis=1)
X.head()
crim | zn | indus | chas | nox | rm | age | dis | rad | tax | ptratio | lstat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0.00632 | 18.0 | 2.31 | 0 | 0.538 | 6.575 | 65.2 | 4.0900 | 1 | 296 | 15.3 | 4.98 |
1 | 0.02731 | 0.0 | 7.07 | 0 | 0.469 | 6.421 | 78.9 | 4.9671 | 2 | 242 | 17.8 | 9.14 |
2 | 0.02729 | 0.0 | 7.07 | 0 | 0.469 | 7.185 | 61.1 | 4.9671 | 2 | 242 | 17.8 | 4.03 |
3 | 0.03237 | 0.0 | 2.18 | 0 | 0.458 | 6.998 | 45.8 | 6.0622 | 3 | 222 | 18.7 | 2.94 |
4 | 0.06905 | 0.0 | 2.18 | 0 | 0.458 | 7.147 | 54.2 | 6.0622 | 3 | 222 | 18.7 | 5.33 |
Then we'll use train_test_split
to get a training set and testing set.
X_train, X_test, y_train, y_test = train_test_split(
X,
y,
test_size=0.2,
random_state=42 # for reproducibility
)
We'll first use the LinearRegression()
class from scikit-learn
. There aren't really any hyperparameters for linear regression, so we'll use it out of the box. The call to the fit()
method trains lin_reg
using the training data.
lin_reg = LinearRegression()
lin_reg.fit(X_train, y_train)
LinearRegression()In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
LinearRegression()
Now that lin_reg
is our trained linear regression model, we can compute the predictions for the train and test data and see how well the model does.
y_pred_train = lin_reg.predict(X_train)
y_pred_test = lin_reg.predict(X_test)
Our outcome metrics are the mean squared error (the average of the squared residuals), and $R^2$. The $R^2$ is the squared correlation between $\mathbf{y}$ and $\mathbf{\hat{y}}$. It's also the percentage of variance of $\mathbf{y}$ that is explained by our predictors, so it has a minimum of 0 showing no prediction and a maximum of 1 showing perfect prediction.
mse_train = mean_squared_error(y_train, y_pred_train)
mse_test = mean_squared_error(y_test, y_pred_test)
print(f"Train MSE = {round(mse_train, 3)}")
print(f"Test MSE = {round(mse_test, 3)}")
Train MSE = 22.604 Test MSE = 22.778
r2_train = r2_score(y_train, y_pred_train)
r2_test = r2_score(y_test, y_pred_test)
print(f"Train R^2 = {round(r2_train, 3)}")
print(f"Test R^2 = {round(r2_test, 3)}")
Train R^2 = 0.74 Test R^2 = 0.689
Our $R^2$ shows us we have a decent model! We were able to explain nearly 69% of the variance in our test set medv
.
The next question is: which of our predictor variables are important for the prediction? For this we can look at the coefficients (the betas in our equation at the top).
pd.Series(lin_reg.coef_, index=lin_reg.feature_names_in_)
crim -0.127196 zn 0.030795 indus 0.025248 chas 3.049477 nox -18.183246 rm 4.240098 age -0.003594 dis -1.465507 rad 0.239990 tax -0.011275 ptratio -0.894845 lstat -0.537770 dtype: float64
scikit-learn
doesn't give additional information about these coefficients out of the box. Some other things that we'd like to know about each coefficient are the standard error, $t$ score, $p$ value, and confidence interval. We can get these pretty easily from the statsmodels
package.
statsmodels
¶The statsmodels
package makes it really easy to do linear regression. We first need to add a constant column to the X_train data for fitting the intercept ($\beta_0$). scikit-learn
did this automatically for us.
X_train_sm = sm.add_constant(X_train)
X_train_sm.head()
const | crim | zn | indus | chas | nox | rm | age | dis | rad | tax | ptratio | lstat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
477 | 1.0 | 15.02340 | 0.0 | 18.10 | 0 | 0.6140 | 5.304 | 97.3 | 2.1007 | 24 | 666 | 20.2 | 24.91 |
15 | 1.0 | 0.62739 | 0.0 | 8.14 | 0 | 0.5380 | 5.834 | 56.5 | 4.4986 | 4 | 307 | 21.0 | 8.47 |
332 | 1.0 | 0.03466 | 35.0 | 6.06 | 0 | 0.4379 | 6.031 | 23.3 | 6.6407 | 1 | 304 | 16.9 | 7.83 |
423 | 1.0 | 7.05042 | 0.0 | 18.10 | 0 | 0.6140 | 6.103 | 85.1 | 2.0218 | 24 | 666 | 20.2 | 23.29 |
19 | 1.0 | 0.72580 | 0.0 | 8.14 | 0 | 0.5380 | 5.727 | 69.5 | 3.7965 | 4 | 307 | 21.0 | 11.28 |
Then we create an instance of OLS
and fit it!
mod = sm.OLS(y_train, X_train_sm)
res = mod.fit()
The summary gives us those additional outputs I mentioned for each input variable. Most of the input variables have large $t$ values and small $p$ values, which demonstrate their importance to the prediction.
print(res.summary())
OLS Regression Results ============================================================================== Dep. Variable: medv R-squared: 0.740 Model: OLS Adj. R-squared: 0.732 Method: Least Squares F-statistic: 92.64 Date: Tue, 24 Jan 2023 Prob (F-statistic): 2.82e-106 Time: 15:44:31 Log-Likelihood: -1203.1 No. Observations: 404 AIC: 2432. Df Residuals: 391 BIC: 2484. Df Model: 12 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ const 36.9705 5.555 6.655 0.000 26.049 47.892 crim -0.1272 0.035 -3.600 0.000 -0.197 -0.058 zn 0.0308 0.016 1.886 0.060 -0.001 0.063 indus 0.0252 0.069 0.368 0.713 -0.110 0.160 chas 3.0495 0.973 3.133 0.002 1.136 4.963 nox -18.1832 4.317 -4.212 0.000 -26.671 -9.696 rm 4.2401 0.471 9.001 0.000 3.314 5.166 age -0.0036 0.015 -0.242 0.809 -0.033 0.026 dis -1.4655 0.231 -6.356 0.000 -1.919 -1.012 rad 0.2400 0.077 3.102 0.002 0.088 0.392 tax -0.0113 0.004 -2.604 0.010 -0.020 -0.003 ptratio -0.8948 0.147 -6.092 0.000 -1.184 -0.606 lstat -0.5378 0.056 -9.538 0.000 -0.649 -0.427 ============================================================================== Omnibus: 128.117 Durbin-Watson: 2.156 Prob(Omnibus): 0.000 Jarque-Bera (JB): 511.080 Skew: 1.355 Prob(JB): 1.05e-111 Kurtosis: 7.797 Cond. No. 1.15e+04 ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. [2] The condition number is large, 1.15e+04. This might indicate that there are strong multicollinearity or other numerical problems.
I love scikit-learn
for all types of analyses, but in the case of linear regression, statsmodels
provides a lot of good information very quickly.
We can also use qualitative predictors in linear regression. When our predictor only has two levels, we can create a single indicator or dummy variable, which just takes a value of 1 for one level and 0 for another. If we have more than two levels, we can dummy code using one less variables than the number of levels.
For example, if we had a predictor that was color with levels Blue, Green, and Red, we could dummy code this with two dummy variables. The first dummy variable could be a 1 if the observation is Blue and a 0 if not. The second is a 1 if Green and 0 if not. From this, a Blue observation is coded as [1 0], Green as [0 1], and Red as [0 0].
In this example, we only deal with the case of predictors with two levels for simplicity (by the way, did you notice we actually had a two-level dummy variable as one of our predictors in the Boston example?).
data_df = pd.read_csv("https://raw.githubusercontent.com/LukeMoraglia/ISLR_datasets/master/data/Carseats.csv")
data_df.head()
Sales | CompPrice | Income | Advertising | Population | Price | ShelveLoc | Age | Education | Urban | US | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 9.50 | 138 | 73 | 11 | 276 | 120 | Bad | 42 | 17 | Yes | Yes |
1 | 11.22 | 111 | 48 | 16 | 260 | 83 | Good | 65 | 10 | Yes | Yes |
2 | 10.06 | 113 | 35 | 10 | 269 | 80 | Medium | 59 | 12 | Yes | Yes |
3 | 7.40 | 117 | 100 | 4 | 466 | 97 | Medium | 55 | 14 | Yes | Yes |
4 | 4.15 | 141 | 64 | 3 | 340 | 128 | Bad | 38 | 13 | Yes | No |
data_df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 400 entries, 0 to 399 Data columns (total 11 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Sales 400 non-null float64 1 CompPrice 400 non-null int64 2 Income 400 non-null int64 3 Advertising 400 non-null int64 4 Population 400 non-null int64 5 Price 400 non-null int64 6 ShelveLoc 400 non-null object 7 Age 400 non-null int64 8 Education 400 non-null int64 9 Urban 400 non-null object 10 US 400 non-null object dtypes: float64(1), int64(7), object(3) memory usage: 34.5+ KB
For our example with the Carseat data, we want to predict Sales
and we're only going to use Price
, Urban
, and US
as our predictors. First let's graph these relationships.
fig, axes = plt.subplots(1, 3, figsize=(15, 5), sharey=True)
sns.scatterplot(ax=axes[0], x="Price", y="Sales", data=data_df)
sns.boxplot(ax = axes[1], x="Urban", y="Sales", data=data_df)
sns.boxplot(ax = axes[2], x="US", y="Sales", data=data_df)
plt.show()
It looks like Price
has a negative relationship with Sales
, and that stores in the US
have higher Sales
than those that aren't. Whether the store is in an Urban
environment or not doesn't seem to make a big difference in Sales
.
Next, we use the get_dummies
functon from pandas
to convert Urban
and US
to dummy variables.
X = data_df[['Price', 'Urban', 'US']]
X = pd.get_dummies(X, drop_first=True)
X.head()
Price | Urban_Yes | US_Yes | |
---|---|---|---|
0 | 120 | 1 | 1 |
1 | 83 | 1 | 1 |
2 | 80 | 1 | 1 |
3 | 97 | 1 | 1 |
4 | 128 | 1 | 0 |
From here, the process is just like our previous example.
y = data_df['Sales']
X_train, X_test, y_train, y_test = train_test_split(
X,
y,
test_size=0.2,
random_state=42
)
lin_reg = LinearRegression()
lin_reg.fit(X_train, y_train)
LinearRegression()In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
LinearRegression()
y_pred_train = lin_reg.predict(X_train)
y_pred_test = lin_reg.predict(X_test)
mse_train = mean_squared_error(y_train, y_pred_train)
mse_test = mean_squared_error(y_test, y_pred_test)
print(f"Train MSE = {round(mse_train, 3)}")
print(f"Test MSE = {round(mse_test, 3)}")
Train MSE = 5.636 Test MSE = 7.775
r2_train = r2_score(y_train, y_pred_train)
r2_test = r2_score(y_test, y_pred_test)
print(f"Train R^2 = {round(r2_train, 3)}")
print(f"Test R^2 = {round(r2_test, 3)}")
Train R^2 = 0.24 Test R^2 = 0.214
pd.Series(lin_reg.coef_, index=lin_reg.feature_names_in_)
Price -0.052808 Urban_Yes -0.109045 US_Yes 1.097797 dtype: float64
Our $R^2$ value is a bit lower than in the Boston
example, with only about 21% of the variance explained.
Just like before, we can get some more insight on the predictors using statsmodels
.
X_train_sm = sm.add_constant(X_train)
mod = sm.OLS(y_train, X_train_sm)
res = mod.fit()
print(res.summary())
OLS Regression Results ============================================================================== Dep. Variable: Sales R-squared: 0.240 Model: OLS Adj. R-squared: 0.232 Method: Least Squares F-statistic: 33.19 Date: Tue, 24 Jan 2023 Prob (F-statistic): 1.12e-18 Time: 15:44:33 Log-Likelihood: -730.73 No. Observations: 320 AIC: 1469. Df Residuals: 316 BIC: 1485. Df Model: 3 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ const 12.9002 0.706 18.268 0.000 11.511 14.290 Price -0.0528 0.006 -9.293 0.000 -0.064 -0.042 Urban_Yes -0.1090 0.292 -0.373 0.709 -0.684 0.466 US_Yes 1.0978 0.278 3.950 0.000 0.551 1.645 ============================================================================== Omnibus: 0.567 Durbin-Watson: 1.949 Prob(Omnibus): 0.753 Jarque-Bera (JB): 0.676 Skew: -0.007 Prob(JB): 0.713 Kurtosis: 2.775 Cond. No. 632. ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
As predicted from our visualizations at the start, the most important predictors of Sales
were Price
and US
.
I hope you enjoyed this overview of linear regression using scikit-learn
and statsmodels
. Look for more posts about machine learning algorithms either coming soon or already available!