Interpreting a Linear Regression Model

This page continues from Linear Regression: we train a linear regression model on the housing data, then interpret it — what its weights, bias, and RMSE mean in dollars, and which features matter most.

Linear Regression on the Housing Data

We use the California housing dataset from the Matplotlib lecture: housing.csv.zip (also on Kaggle). Unzip it into your ml4biz folder and run the code below in a notebook there. We use scikit-learn to run linear regression on it:

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

df = pd.read_csv('housing.csv')

# for simplicity, only use population, median_income, and housing_median_age to train the model
X_housing = df[['population', 'median_income', 'housing_median_age']]
y_housing = df['median_house_value']

# split train and test sets
X_train_housing, X_test_housing, y_train_housing, y_test_housing = train_test_split(
    X_housing, y_housing, test_size=0.2, random_state=42)

# train the linear regression model
lin_reg_housing = LinearRegression()
lin_reg_housing.fit(X_train_housing, y_train_housing)

# the fitted model is
# median_house_value = -17768 + 2.975 * population + 43399 * median_income + 1828 * housing_median_age
# 1.82807894e+03 is scientific notation for 1828.07894: https://www.mathsisfun.com/numbers/scientific-notation.html
print(f'The intercept is {lin_reg_housing.intercept_} and coefficient is {lin_reg_housing.coef_}')

# calculate MSE and RMSE on the test set
# RMSE is measured on the same scale with the same units as y - roughly the "average prediction error"
y_pred_housing = lin_reg_housing.predict(X_test_housing)
lin_mse_housing = mean_squared_error(y_test_housing, y_pred_housing)
lin_rmse_housing = np.sqrt(lin_mse_housing)

print(f'Median house value mean is {df.median_house_value.mean():.3f}, std is: {df.median_house_value.std():.3f}')
print(f'The RMSE is {lin_rmse_housing:.3f}, which is ~{lin_rmse_housing/df.median_house_value.std():.2f} std')

# train the final model using ALL data
# the intercept and coefficients are slightly different from the training data result
lin_reg_housing.fit(X_housing, y_housing)
print(f'The intercept is {lin_reg_housing.intercept_} and coefficient is {lin_reg_housing.coef_}')

Output:
The intercept is -17768.502521458664 and coefficient is [2.97530968e+00 4.33989999e+04 1.82807894e+03]
Median house value mean is 206855.817, std is: 115395.616
The RMSE is 81323.636, which is ~0.70 std
The intercept is -17362.509966188372 and coefficient is [3.16069821e+00 4.32273443e+04 1.82943214e+03]

A typical prediction error is about $81,000, about 0.7 of the standard deviation of house values.

Interpreting the Model

The model trained on the training set is:

predicted median_house_value = 2.975 × population
                             + 43399 × median_income
                             + 1828 × housing_median_age
                             - 17768

The Weights

Each weight $w_j$ says how much the prediction changes when feature $j$ increases by one unit, holding the other features fixed. To read a weight, you need the feature's unit:

• median_income is measured in tens of thousands of dollars (a value of 3.5 means $35,000). So a district whose median income is $10,000 higher is predicted to have houses worth $43,399 more.
• housing_median_age is in years: each extra year adds $1,828 to the predicted value.
• population is in people: each extra person adds $2.98, or about $2,975 per 1,000 people.

Three things to keep in mind:

• A small weight does not mean an unimportant feature. The size of a weight depends on the feature's unit. Population's weight is tiny because one person is a tiny change; measured in thousands of people, the same weight would be 2,975. To compare weights fairly, standardize the features first, so one unit means one standard deviation for every feature.
• A weight is an association, not a cause. The model says districts with higher incomes have more expensive houses. It does not say that raising people's income would raise house prices.
• Correlated features make weights unstable. If two features carry nearly the same information, the model can trade weight between them without changing its predictions, so the individual weights become hard to interpret. When the Inverse Fails shows the extreme case of a duplicated feature.

The Bias

The bias $b = -17768$ is the prediction when every feature is 0: a district with no people, no income, and brand-new houses. No such district exists, so a negative house value here is not a problem. The bias sets the height of the line (or plane) so that predictions fit the data where the data actually is. It is only meaningful on its own when "all features are 0" is a realistic case.

The RMSE and a Baseline

RMSE is measured in the same unit as the target, so the model's typical prediction error is about $81,000. Is that good? Compare it with a model-free baseline: ignore the features and predict the average house value for every district.

from sklearn.metrics import mean_absolute_error

# model-free baseline: predict the training set's average house value for every district
y_pred_baseline = np.full(len(y_test_housing), y_train_housing.mean())
baseline_rmse = np.sqrt(mean_squared_error(y_test_housing, y_pred_baseline))
print(f'Baseline RMSE is {baseline_rmse:.3f}')

# MAE: the average absolute error, without squaring
print(f'MAE is {mean_absolute_error(y_test_housing, y_pred_housing):.3f}')

Output:
Baseline RMSE is 114485.635
MAE is 60457.311

• The baseline's RMSE is $114,486, about the standard deviation of house values. That is why comparing RMSE with the standard deviation is useful: our three features cut the typical error by about 29%, from $114,486 to $81,324.
• RMSE squares the errors before averaging, so large errors count more. The average absolute error (MAE) is only $60,457; a few badly predicted districts push RMSE up to $81,324.
• Whether $81,000 is good enough depends on the decision. It may be fine for comparing neighborhoods, but too large for pricing one house. A model with more features would do better.

There is no general "good" RMSE or MAE, because both are in the target's unit: an error of 5 would be useless when predicting a product's star rating (1 to 5), but excellent when predicting a website's daily visits (tens of thousands). Always judge them against a baseline, the target's typical size, or the cost of a mistake.

Feature Importance

Which feature matters most? Linear regression has no feature_importances_ like a decision tree, and the raw weights are not importance scores: as The Weights shows, their size depends on each feature's unit. There are two common ways to measure importance.

1. Standardized weights. Standardize the features before fitting. Then every feature is measured in the same unit, one standard deviation, and the weights become comparable:

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# standardize the features, then fit: each weight is now the change in the
# prediction when that feature increases by one standard deviation
std_pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('lin_reg', LinearRegression()),
])
std_pipeline.fit(X_train_housing, y_train_housing)

for name, weight in zip(X_housing.columns, std_pipeline['lin_reg'].coef_):
    print(f'{name}: {weight:.0f}')

Output:
population: 3383
median_income: 82642
housing_median_age: 23038

A district whose median income is one standard deviation higher is predicted to be worth $82,642 more, compared with $23,038 for house age and $3,383 for population. The larger the absolute value, the more important the feature, and the sign shows the direction.

2. Permutation importance. Take the trained model and shuffle one feature's column in the test set: randomly reorder its values while leaving every other column, and the target, where they were. For example, shuffling median_income:

The column has the same values as before, but each one now sits next to the wrong house: district A looks poor and district D looks rich. Then predict again with the same model, without retraining:

• If the model relies on the feature, the predictions go badly wrong and the RMSE jumps.
• If the model barely uses the feature, scrambling it changes almost nothing.

Shuffling, rather than deleting the column, keeps the same trained model and the same kinds of values. It destroys only the relationship between the feature and the price, which is exactly what we want to measure. Because one random shuffle could leave some values near their original places, we shuffle several times and average:

from sklearn.inspection import permutation_importance

# shuffle each feature 10 times on the test set and average the increase in RMSE
result = permutation_importance(
    std_pipeline, X_test_housing, y_test_housing,
    scoring='neg_root_mean_squared_error', n_repeats=10, random_state=42)

for name, increase in zip(X_housing.columns, result.importances_mean):
    print(f'{name}: RMSE increases by {increase:.0f}')

Output:
population: RMSE increases by 162
median_income: RMSE increases by 58258
housing_median_age: RMSE increases by 6247

Scrambling median income adds $58,258 to the typical error; scrambling population adds only $162. Both methods agree: income matters most, then house age, and population barely matters.

Standardized weights are cheap and also show the direction (+ or −), but they are only valid for linear models. Permutation importance is model-agnostic, measured on test data, and in business units, so it is the better general-purpose tool, and it lets you compare a linear model and a tree on the same footing.

Both methods share one weakness: if two features carry nearly the same information, they split the credit between them, and each looks less important than it really is.

Questions

• Do you still need Grid Search to tune a linear regression model? Why?
• Do you still need to build a pipeline to do the necessary pre-processing?