Regression Metrics and Reading a Residual

Context

A residual is one listing's price minus the model's prediction for it, and a regression metric condenses the residuals of a whole set of listings into one number, such as the average size of a miss. A residual plot draws every residual against its prediction, so the pattern a single number averages away stays in view.

Every score a pricing model is judged by is built this way. A team deciding whether a model is good enough to show sellers reads its average miss in currency, a reviewer choosing between two models reads the same metric for both on the same held-out listings, and a modeler choosing the next feature reads the residuals split by neighborhood or by the size of the prediction.

Most metric mistakes are silent. A score computed on the training rows describes listings the model has already seen, a percentage error of 0.243 is reported as a quarter of a percent, and two R² values (the share of the price spread a model explains) from different listings are compared as if they measured the same thing. Residuals that average zero on the training rows look like a model that leans neither way, when every least squares fit with an intercept leaves them so. Lesson 3.1.3 separated the loss a fit minimizes from the number it reports, and this lesson computes that number and reads it.

Goal

• Compute a model's residuals on the training and test rows, and explain why the training residuals average zero while the test residuals need not.
• Compute MAE, MSE and RMSE with sklearn.metrics and by hand with NumPy, and explain why RMSE is never below MAE.
• Compute R² with r2_score and .score(), set it against a DummyRegressor baseline, and explain why two R² values from different rows cannot be compared.
• Draw a residual plot, check what it holds, and diagnose a fan and a missing feature from the residuals' spread and group means.
• Choose the metric that answers the question being asked, and report MAPE as the fraction it is.

Prerequisites

Optional, and the lesson stands without them. They are here if you want the idea behind R² drawn out, or every regression metric's formula in one place, before you start.

• YouTube - "R-squared, Clearly Explained!!!" by StatQuest with Josh Starmer. R² built up step by step from the variation around the mean and the variation around a fitted line.
• Official docs - scikit-learn User Guide, "Regression metrics", section 3.4.6 of the 1.9 documentation. The formula behind each regression metric beside the function that computes it, R², MAE, MSE and MAPE among them.

Table of Contents

• Residuals
• MAE, MSE and RMSE
• R² and the Baseline It Is Measured Against
• Reading a Residual Plot
• Choosing the Metric That Matches the Question

Theory

Nothing to install, scikit-learn has been in (bootcamp) since Lesson 3.1.1. Every Theory block below goes into one script, metrics.py, saved in ~/dev/bootcamp beside listings.csv and run with python3 metrics.py, each block added under the one before it.

Residuals

What. A residual is one row's target minus the model's prediction for that row, 2.3.4's residual, one per listing. A positive residual is a listing priced above its prediction, so the model guessed too low, and a negative one is a listing it guessed too high. Why. Every metric in this lesson summarizes the residuals, and the residuals keep what a summary drops: which listings the model misses, in which direction, and by how much. When. After every fit, on the test rows first, because the training rows are the ones the fit chose its parameters on.

The first block rebuilds Lesson 3.1.1's price-on-area model on Lesson 3.1.2's split.

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

listings = pd.read_csv("listings.csv")
data = listings.dropna(subset=["area_sqm"])
X = data[["area_sqm"]]
y = data["price"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
model = LinearRegression().fit(X_train, y_train)

Subtracting the array predict returns from y_test gives a Series that keeps each test listing's label, the labels 3.1.2 showed the split keeps.

train_pred = model.predict(X_train)
test_pred = model.predict(X_test)
train_res = y_train - train_pred
test_res = y_test - test_pred
print(len(train_res), len(test_res))  # 2300 575

# each test listing's row, its prediction and residual beside it
check = data.loc[X_test.index].assign(predicted=test_pred, residual=test_res)
cols = ["neighborhood", "area_sqm", "price", "predicted", "residual"]
print(check[cols].head(3).round({"predicted": 0, "residual": 0}))

     neighborhood  area_sqm   price  predicted  residual
765     Southgate     123.2  342000   461717.0 -119717.0
2542    Southgate     108.6  391000   413314.0  -22314.0
1346     Eastgate     107.5  528000   409667.0  118333.0

data.loc[X_test.index] takes each test listing's row by its label, and 2.6.1's assign adds the prediction and the residual. The second and third listings are 1.1 square meters apart, so the model, which sees only the area, predicts nearly the same price for both and misses them in opposite directions, 22,314 too high for the Southgate listing and 118,333 too low for the Eastgate one.

On the training rows the residuals sum to zero up to float error, and so does their mean. That is 2.3.4's column of ones, which fit_intercept=True adds, as Lesson 3.1.1 showed: least squares leaves the residuals at a right angle to every column it was given, so their dot product with the ones, which is their sum, is zero whatever the prices are. The check is on the mean, because a sum of 2,300 residuals, a third of them over 100,000 in size, can carry more rounding than np.isclose's fixed allowance of 1e-08 against zero, the limit 2.3.4 measured.

print(np.isclose(train_res.mean(), 0))  # True
print(round(test_res.sum()))            # -3530921
print(round(test_res.mean()))           # -6141

Nothing holds the test rows' sum at zero, because they had no part in the fit. Their mean of -6,141 says the model predicts new listings about 6,141 too high on average.

Without the column of ones, nothing holds even the training mean at zero:

origin = LinearRegression(fit_intercept=False).fit(X_train, y_train)
origin_res = y_train - origin.predict(X_train)
print(round(origin_res.mean()))  # 9031

Forced through zero price at zero area, the line leaves its training residuals averaging 9,031, so a zero training mean reports the intercept, as 2.3.4 found on five listings, and nothing about the fit.

DON'T:

• DON'T read the training residuals' zero mean as a good fit. Every LinearRegression with fit_intercept=True leaves one on its training rows, however badly it fits them.
• DON'T test the training sum with np.isclose(train_res.sum(), 0). Its allowance against zero is a fixed 1e-08, and a sum of 2,300 residuals this size can carry more rounding than that, so a fit whose residuals do sum to zero can print False.
• DON'T compute a residual as the prediction minus the price. Every sign flips, and a listing the model priced too low reads as one it priced too high.

MAE, MSE and RMSE

What. Three summaries of the residuals in sklearn.metrics. MAE, the mean absolute error, averages their sizes. MSE, the mean squared error, averages their squares. RMSE, the root mean squared error, is the square root of MSE, which brings it back to currency. Why. 575 residuals cannot be read one at a time, and each of the three condenses them for a different question. Lesson 3.1.2's score, the average size of a miss computed with NumPy, is MAE. When. Every time a regression model is scored, on the test rows, with MAE and RMSE side by side.

Each function takes the targets first and the predictions second.

from sklearn.metrics import (
    mean_absolute_error,
    mean_squared_error,
    root_mean_squared_error,
)

mae = mean_absolute_error(y_test, test_pred)
mse = mean_squared_error(y_test, test_pred)
rmse = root_mean_squared_error(y_test, test_pred)
# 87,952 13,631,340,909 116,753
print(f"{mae:,.0f} {mse:,.0f} {rmse:,.0f}")

By hand, each is one line of 2.2.3's NumPy on the residuals:

mae_hand = np.abs(test_res).mean()
mse_hand = (test_res ** 2).mean()
rmse_hand = np.sqrt(mse_hand)
# 87,952 13,631,340,909 116,753
print(f"{mae_hand:,.0f} {mse_hand:,.0f} {rmse_hand:,.0f}")

The average test listing is missed by 87,952. RMSE says 116,753 in the same currency, and MSE's 13,631,340,909 is in 3.1.3's squared currency, which RMSE's root undoes.

RMSE is never below MAE. Squaring before averaging weights every miss by its own size, so a large miss counts for more than its share, and the two are equal only when every miss has the same size. Two sets of four misses, both with an MAE of 100, show the two ends:

even = np.array([100.0, -100.0, 100.0, -100.0])
one_big = np.array([0.0, 0.0, 0.0, 400.0])
for misses in [even, one_big]:
    print(np.abs(misses).mean(), np.sqrt((misses ** 2).mean()))

100.0 100.0
100.0 200.0

One miss of 400 beside three perfect predictions doubles RMSE and leaves MAE where it was. On the listings, the 11 test listings missed by more than 300,000 play that part:

big = np.abs(test_res) > 300_000
print(big.sum())                            # 11
kept = (~big).to_numpy()                    # a plain array, picks by position
kept_mae = mean_absolute_error(y_test[kept], test_pred[kept])
kept_rmse = root_mean_squared_error(y_test[kept], test_pred[kept])
print(f"{kept_mae:,.0f} {kept_rmse:,.0f}")  # 82,127 104,357
# 28,802 22,230
print(f"{rmse - mae:,.0f} {kept_rmse - kept_mae:,.0f}")

Leaving out about 2 percent of the test listings moves RMSE from 116,753 to 104,357 and MAE only from 87,952 to 82,127, and the gap between them falls from 28,802 to 22,230. The gap grows with how unequal the misses are, and the largest misses weigh most in it, 11 listings carrying almost a quarter of it, which is why a report gives both.

Many tutorials ask mean_squared_error for the root with squared=False. In scikit-learn 1.9.1 the function has no such keyword, and the call raises a TypeError before computing anything:

try:
    mean_squared_error(y_test, test_pred, squared=False)
except TypeError as err:
    print(err)  # got an unexpected keyword argument 'squared'

DON'T:

• DON'T pass squared=False to mean_squared_error. scikit-learn 1.9.1 raises TypeError: got an unexpected keyword argument 'squared', and root_mean_squared_error returns the root.
• DON'T quote RMSE as the average miss. The average miss is MAE, 87,952, and RMSE's 116,753 carries the extra weight of the largest misses.

R² and the Baseline It Is Measured Against

What. R², the coefficient of determination, is one minus the model's sum of squared residuals over the sum of squared residuals left by predicting the evaluated rows' own mean price for every row. 1 is a perfect fit, 0 is no better than that mean, and below 0 is worse than it. Why. MAE and RMSE say how large the misses are, not whether a model that ignored every feature would miss by as much. R² sets the model against exactly that baseline (the simplest prediction a model has to beat). When. Beside a currency metric, never instead of one, and only beside another R² computed on the same rows.

from sklearn.metrics import r2_score

r2 = r2_score(y_test, test_pred)
sse = (test_res ** 2).sum()                  # the model's squared misses
sst = ((y_test - y_test.mean()) ** 2).sum()  # the test mean's squared misses
test_score = model.score(X_test, y_test)
train_score = model.score(X_train, y_train)
print(round(r2, 3))                          # 0.565
print(round(1 - sse / sst, 3))               # 0.565
print(round(test_score, 3))                  # 0.565
print(round(train_score, 3))                 # 0.546

r2_score and the formula agree, and a fitted regressor's .score(X, y), unless the regressor is one of Lesson 3.6.3's hyperparameter searches told to score by another metric, returns the R² of its predictions for X against y, the same 0.565. The training rows give 0.546, below the test score on this split and above it on other seeds, one split's draw.

scikit-learn's DummyRegressor (an estimator that ignores the features and predicts one number for every row, the mean of the training targets unless told otherwise) is that baseline as an estimator:

from sklearn.dummy import DummyRegressor

baseline = DummyRegressor().fit(X_train, y_train)
base_pred = baseline.predict(X_test)
base_r2 = baseline.score(X_test, y_test)
base_mae = mean_absolute_error(y_test, base_pred)
print(round(base_pred[0]))    # 380125
print(round(y_train.mean()))  # 380125
print(round(y_test.mean()))   # 381130
print(f"{base_r2:.6f}")       # -0.000032
print(f"{base_mae:,.0f}")     # 139,562

It predicts the training mean, 380,125, for every listing and misses a test listing by 139,562 on average, against the model's 87,952. Its test R² sits just below zero, and no constant can score above zero: R² measures against the evaluated rows' own mean, 381,130 here, and the mean of a set of prices is the constant with the smallest squared error on them, as 3.1.3's scan showed, so every other constant scores below zero.

R² is measured against the spread of the rows it is computed on, so R² values from different rows are measured against different things, and R² has no floor. The same predictions, scored on three neighborhoods' test listings alone:

for name in ["Lakeview", "Northend", "Hillcrest"]:
    rows = check[check["neighborhood"] == name]
    hood_r2 = r2_score(rows["price"], rows["predicted"])
    hood_mae = mean_absolute_error(rows["price"], rows["predicted"])
    spread = rows["price"].std()
    print(name, round(hood_r2, 3), f"{hood_mae:,.0f}", f"{spread:,.0f}")

Lakeview 0.47 117,885 203,062
Northend 0.108 121,490 149,286
Hillcrest -0.165 120,665 126,465

The model misses all three by about the same in currency, near 120,000, and their R² values run from 0.47 down to -0.165. Lakeview's prices spread widest, a standard deviation of 203,062, so the same miss is the smallest share of its spread. Hillcrest's spread least, and on its own listings the model does worse than Hillcrest's mean price would, a negative R².

DON'T:

• DON'T compare R² values computed on different rows. Lakeview's 0.47 and Hillcrest's -0.165 come from misses of nearly the same size.
• DON'T read the baseline's -0.000032 as a broken estimator. A constant fitted on the training rows scores below zero on any rows whose mean differs from the training mean.
• DON'T report model.score(X_train, y_train) as the model's R². Its 0.546 describes the rows the fit chose its parameters on.

Reading a Residual Plot

What. A residual plot draws one point per listing, its prediction on the x axis and its residual on the y axis, with a line at zero. A model whose misses are noise leaves an even band around that line, and any shape in the band is a miss the model makes the same way again and again. Why. MAE, RMSE and R² compress 575 residuals into one number each, and an even band, a fan, a curve and two separate clouds can all compress into the same number. The plot shows which one it is. When. After every fit, on the test rows, before any metric is reported as the model's.

The chart is 2.7.2's Figure and Axes, saved and checked by printing what it holds:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 5))
points = ax.scatter(check["predicted"], check["residual"], alpha=0.4)
ax.axhline(0, color="black", linestyle="--")
ax.set_xlabel("Predicted price")
ax.set_ylabel("Residual, price minus prediction")
ax.set_title("Price on area, test listings")
fig.savefig("residuals.png", dpi=150, bbox_inches="tight")
print(len(points.get_offsets()))  # 575
print(len(ax.lines))              # 1
plt.close(fig)

ax.axhline(0) draws a horizontal line at height 0 across the whole Axes, the twin of the axvline in 2.7.2's Assignment, and joins ax.lines. The 575 points are every test listing.

The cloud hugs the zero line at the left, where the predictions are small, and spreads farther above and below it as the prediction grows, a fan. The spread in each quarter of the predictions, cut with 2.6.1's pd.qcut, measures what the eye sees:

labels = ["q1", "q2", "q3", "q4"]
check["quarter"] = pd.qcut(check["predicted"], q=4, labels=labels)
by_quarter = check.groupby("quarter")["residual"].agg(["count", "mean", "std"])
print(by_quarter.round(0))

         count     mean       std
quarter                          
q1         144 -13472.0   55408.0
q2         144   4733.0  101261.0
q3         143  -7105.0  125068.0
q4         144  -8726.0  160070.0

Each quarter's mean stays within 13,472 of zero while the standard deviation nearly triples, from 55,408 to 160,070: the misses grow with the price, and no single currency figure describes both ends. The fix you can already name is 2.2.3's log transform, on whose scale a miss of the same proportion has the same size at every price, and Lesson 3.3.1 applies it to the price and the area together.

A pattern tied to a column the model never saw is measured rather than read off the plot, here with 2.6.2's groupby:

by_hood = check.groupby("neighborhood")["residual"].agg(["count", "mean"])
print(by_hood.sort_values("mean").round(0))

              count      mean
neighborhood                 
Hillcrest        73 -120177.0
Northend         67 -107643.0
Southgate        96  -54403.0
Downtown        125  -13549.0
Riverside        80   28789.0
Lakeview         67  109482.0
Eastgate         67  145253.0

The model sees only the area, so it prices two listings of the same area alike whatever their neighborhoods, and the residuals carry the difference: Eastgate's listings sit 145,253 above their predictions on average and Hillcrest's 120,177 below. Residual means far from zero across the groups of a column the model does not have mark a missing feature, here a category Lesson 3.2.2 turns into columns a model can take.

The x axis holds the prediction and not the price, for a reason 2.7.1's corr shows:

with_pred = check["residual"].corr(check["predicted"])
with_price = check["residual"].corr(check["price"])
print(round(with_pred, 3))   # -0.051
print(round(with_price, 3))  # 0.618

Against the prediction the residuals barely move. Against the price they climb, because the price is inside every residual: a listing priced high has a high residual unless its prediction rose just as far.

Four shapes are worth knowing on sight. An even band, the same width from end to end, is noise the features cannot explain, the healthy case. A fan that widens as the prediction grows, this model's shape, says the misses grow with the prediction. A curve, above zero at both ends and below in the middle or the reverse, says the model is straight where the data bends. Separate clouds, one above zero and another below, are groups the features cannot tell apart, the shape a missing category draws once its groups sit far enough apart.

DON'T:

• DON'T read the test mean of -6,141 as a model that leans neither way. The neighborhood means run from -120,177 to 145,253 and cancel in it.
• DON'T quote one currency miss for listings at every price. The spread grows from 55,408 in the lowest quarter of the predictions to 160,070 in the highest.
• DON'T put the price on the x axis. The price is inside every residual, so the cloud tilts upward along it, a correlation of 0.618 against -0.051 on the prediction axis.

Choosing the Metric That Matches the Question

What. Each metric answers one question. MAE: how far off is a price on average, in currency. RMSE: the same, with every miss weighted by its own size. MAPE, the mean absolute percentage error: how far off as a share of each listing's own price. R²: how much of the spread in prices the model explains beyond the evaluated rows' mean, and nothing about money. Why. Two models can rank in opposite orders under two metrics, so the metric decides which model wins, and only the question being asked can decide the metric. When. Before the models are scored, so the question fixes the metric and no one picks it afterwards by which number looks best.

from sklearn.metrics import mean_absolute_percentage_error

mape = mean_absolute_percentage_error(y_test, test_pred)
mape_hand = (np.abs(test_res) / y_test).mean()
print(round(mape, 3))       # 0.243
print(round(mape_hand, 3))  # 0.243
print(f"{mape:.1%}")        # 24.3%

mean_absolute_percentage_error divides each miss by its listing's price, averages the shares and returns a fraction: 0.243 means the average listing is missed by 24.3 percent of its price, which Module 1's :.1% format prints as a percentage.

Lesson 3.1.3 fitted the same feature under absolute error with QuantileRegressor(quantile=0.5, alpha=0.0), median regression. Both models, on the same test rows, by all four metrics:

from sklearn.linear_model import QuantileRegressor

median_model = QuantileRegressor(quantile=0.5, alpha=0.0)
median_model.fit(X_train, y_train)
for name, fitted in [("least squares", model), ("median", median_model)]:
    pred = fitted.predict(X_test)
    print(
        name,
        f"{mean_absolute_error(y_test, pred):,.0f}",
        f"{root_mean_squared_error(y_test, pred):,.0f}",
        round(mean_absolute_percentage_error(y_test, pred), 3),
        round(r2_score(y_test, pred), 3),
    )

least squares 87,952 116,753 0.243 0.565
median 85,031 117,789 0.22 0.558

Median regression misses the average listing by less, 85,031 against 87,952, and by a smaller share of its price, 0.22 against 0.243. Least squares has the lower RMSE, 116,753 against 117,789, and the higher R², 0.565 against 0.558. Here each model wins the two metrics built on the error its own fit minimized, so which model is better has no answer until someone says what a miss costs.

The question decides. A seller asking how far off an estimate is likely to be is answered by MAE, in currency. A team for whom one miss of 400,000 does more damage than four of 100,000 reads RMSE. A team comparing estimates for apartments near 150,000 and houses near 1,500,000 reads MAPE, because a miss of 50,000 is a different failure at each price. A modeler asking whether area explains price at all reads R². Each choice is made before the models are scored.

DON'T:

• DON'T read MAPE's 0.243 as a quarter of a percent. It is a fraction, 24.3 percent, and f"{mape:.1%}" prints it as one.
• DON'T choose the metric after scoring the models. Least squares wins RMSE and R² and median regression wins MAE and MAPE, so a metric picked afterwards can crown either.
• DON'T swap the arguments of r2_score or mean_absolute_percentage_error. Both take the targets first, and with the arguments swapped they return other numbers, 0.298 and 0.224, without an error.

Lab Exercises

Lab Exercise 1: Distance to the Center, in Currency

Work in distance_metrics.py in ~/dev/bootcamp. Load listings.csv, keep the listings that have an area, and split area_sqm and distance_km against price with test_size=0.2 and random_state=42. The split shuffles the same 2,875 rows with the same seed as the Theory, so its test rows are the Theory's 575 listings. Fit a LinearRegression.

Check that the training residuals average zero up to float error, and print the test residuals' sum and mean in whole units, with a comment saying which way the model leans on new listings.

Print MAE, MSE and RMSE from sklearn.metrics, the same three by hand from the residuals, and the gap between RMSE and MAE, and compare them in a comment with the area model's 87,952, 116,753 and 28,802.

Finally, count the test listings this model misses by more than 300,000, and print MAE, RMSE and their gap without them.

Lab Exercise 2: What Distance Explains, and What It Leaves

Work in distance_check.py in ~/dev/bootcamp, with the rows and split of Lab Exercise 1.

Fit price on distance_km alone and fit a DummyRegressor, then print each one's test R² to six decimal places and its MAE in whole units, with a comment on what distance explains on its own.

Fit Lab Exercise 1's two-feature model and print its test R² from r2_score and from .score(). Draw its residual plot the way the Theory did, save it as distance_residuals.png, and print how many points and lines it holds. Print the residual standard deviation in each quarter of the predictions and the residual mean per neighborhood, lowest first.

Finish with the model's MAPE as a percentage with one decimal, and a comment saying whether distance removed the fan and the neighborhood pattern.

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import (
    mean_absolute_error,
    mean_squared_error,
    root_mean_squared_error,
)
from sklearn.model_selection import train_test_split

data = pd.read_csv("listings.csv").dropna(subset=["area_sqm"])
X = data[["area_sqm", "distance_km"]]
y = data["price"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
model = LinearRegression().fit(X_train, y_train)

train_res = y_train - model.predict(X_train)
test_pred = model.predict(X_test)
test_res = y_test - test_pred
print(np.isclose(train_res.mean(), 0))  # True
print(round(test_res.sum()))            # -2963356
print(round(test_res.mean()))           # -5154
# it predicts new listings about 5,154 too high on average

mae = mean_absolute_error(y_test, test_pred)
mse = mean_squared_error(y_test, test_pred)
rmse = root_mean_squared_error(y_test, test_pred)
# 83,663 12,393,113,124 111,324
print(f"{mae:,.0f} {mse:,.0f} {rmse:,.0f}")

mae_hand = np.abs(test_res).mean()
mse_hand = (test_res ** 2).mean()
# 83,663 12,393,113,124 111,324
print(f"{mae_hand:,.0f} {mse_hand:,.0f} {np.sqrt(mse_hand):,.0f}")
print(f"{rmse - mae:,.0f}")  # 27,662
# area alone: 87,952, 116,753 and 28,802 on the same 575 listings,
# so distance lowers both misses and barely moves the gap

big = np.abs(test_res) > 300_000
print(big.sum())  # 11
kept = (~big).to_numpy()
kept_mae = mean_absolute_error(y_test[kept], test_pred[kept])
kept_rmse = root_mean_squared_error(y_test[kept], test_pred[kept])
gap = kept_rmse - kept_mae
# 77,811 98,267 20,457
print(f"{kept_mae:,.0f} {kept_rmse:,.0f} {gap:,.0f}")

import matplotlib.pyplot as plt
import pandas as pd
from sklearn.dummy import DummyRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import (
    mean_absolute_error,
    mean_absolute_percentage_error,
    r2_score,
)
from sklearn.model_selection import train_test_split

data = pd.read_csv("listings.csv").dropna(subset=["area_sqm"])
X = data[["area_sqm", "distance_km"]]
y = data["price"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

dist_model = LinearRegression().fit(X_train[["distance_km"]], y_train)
dist_pred = dist_model.predict(X_test[["distance_km"]])
baseline = DummyRegressor().fit(X_train, y_train)
base_pred = baseline.predict(X_test)
for name, pred in [("distance", dist_pred), ("baseline", base_pred)]:
    r2 = r2_score(y_test, pred)
    mae = mean_absolute_error(y_test, pred)
    print(name, f"{r2:.6f}", f"{mae:,.0f}")
# distance alone explains almost none of the spread, and its
# average miss is within a few hundred of the baseline's

distance 0.001401 139,441
baseline -0.000032 139,562

model = LinearRegression().fit(X_train, y_train)
test_pred = model.predict(X_test)
test_r2 = r2_score(y_test, test_pred)
score = model.score(X_test, y_test)
print(round(test_r2, 3))  # 0.605
print(round(score, 3))    # 0.605

check = data.loc[X_test.index].assign(
    predicted=test_pred, residual=y_test - test_pred
)
fig, ax = plt.subplots(figsize=(8, 5))
points = ax.scatter(check["predicted"], check["residual"], alpha=0.4)
ax.axhline(0, color="black", linestyle="--")
ax.set_xlabel("Predicted price")
ax.set_ylabel("Residual, price minus prediction")
ax.set_title("Price on area and distance, test listings")
fig.savefig("distance_residuals.png", dpi=150, bbox_inches="tight")
# 575 1
print(len(points.get_offsets()), len(ax.lines))
plt.close(fig)

labels = ["q1", "q2", "q3", "q4"]
check["quarter"] = pd.qcut(check["predicted"], q=4, labels=labels)
print(check.groupby("quarter")["residual"].std().round(0))
by_hood = check.groupby("neighborhood")["residual"].mean()
print(by_hood.sort_values().round(0))

quarter
q1     65467.0
q2     93514.0
q3    117934.0
q4    150862.0
Name: residual, dtype: float64
neighborhood
Hillcrest    -93202.0
Downtown     -63884.0
Northend     -62914.0
Southgate    -15780.0
Riverside      3157.0
Lakeview     118397.0
Eastgate     139864.0
Name: residual, dtype: float64

mape = mean_absolute_percentage_error(y_test, test_pred)
print(f"{mape:.1%}")  # 23.0%
# the spread still grows from 65,467 to 150,862, a fan, and the
# neighborhood means still run from -93,202 to 139,864: distance
# changed the misses without removing either pattern

Theory Assessment

Q1: Two Listings of the Same Size

The Theory's second and third test listings, 108.6 square meters in Southgate at 391,000 and 107.5 square meters in Eastgate at 528,000, have residuals of -22,314 and 118,333. Say what each sign says about its prediction, and why two listings of almost the same area are missed in opposite directions. Then explain why the 2,300 training residuals average zero up to float error, the 575 test residuals average -6,141, and the training residuals of LinearRegression(fit_intercept=False) average 9,031.

Hint: A residual is the price minus the prediction, and the model's only input is the area. For the rest, ask which column fit_intercept=True adds, what least squares guarantees about every column it was given, and on which rows it guarantees it.

Q2: Squares, Roots and a Constant

On the test rows MAE is 87,952 and RMSE 116,753. Explain why RMSE can never come out below MAE, what the misses would have to look like for the two to be equal, and why leaving out 11 listings moves RMSE from 116,753 to 104,357 but MAE only from 87,952 to 82,127. Then explain why DummyRegressor() scores -0.000032 on the test rows, and why no constant prediction can score above zero there.

Hint: Compare the misses [100, -100, 100, -100] and [0, 0, 0, 400]. For the second part, ask which constant has the smallest squared error on a set of prices, and whose mean R² measures against.

Q3: A Fan, Seven Neighborhoods and the Wrong Axis

In the area model's residual plot, the residuals' standard deviation grows from 55,408 in the lowest quarter of the predictions to 160,070 in the highest, while each quarter's mean stays within 13,472 of zero, and the neighborhood means run from -120,177 in Hillcrest to 145,253 in Eastgate. Say what each finding says about the model, which of the four shapes each one is, and what fix each points to. Then say what the same residuals plotted against the price would have shown, and why.

Hint: One finding is about how wide the cloud is at each prediction, the other about where groups of listings sit relative to zero. For the last part, write the residual out as a subtraction and ask what it does as the price grows.

Q4: A Report That Reads Every Number

A teammate's report on the area model makes four claims.

1. The model leans neither high nor low on new listings, because np.isclose(train_res.mean(), 0) prints True.
2. Its MAPE of 0.243 means it misses a listing's price by about a quarter of a percent.
3. It prices Lakeview listings far more accurately than Northend ones, R² 0.47 against 0.108.
4. RMSE sits 28,802 above MAE, so a few listings are missed by far more than the rest.

Say which claims hold. For each one that does not, give the number or the check that refutes it.

Hint: Every claim quotes a number this lesson printed. Ask, for each, which rows the number was computed on, what it is measured against, and what unit it is in.

gap = y_train.mean() - y_test.mean()
shortfall = len(y_test) * gap ** 2 / sst
print(round(gap))           # -1006
print(f"{-shortfall:.6f}")  # -0.000032

Assignment

Assignment: A Residual Report the Reviewer Can Check

Before the area model goes to the pricing team, a reviewer asks for two things, each with a printed check, in residual_report.py in ~/dev/bootcamp, built from the Theory's rows, split and model.

The first is one figure, saved once as residual_report.png, 14 by 5 inches, holding two titled residual plots side by side, the training rows on the left and the test rows on the right. The reviewer wants scikit-learn's own residual chart, not a hand-drawn one: PredictionErrorDisplay, imported from sklearn.metrics, has a method from_predictions(y_true, y_pred, kind="residual_vs_predicted", ax=ax) that draws the residuals against the predictions on the Axes passed as ax=, with a dashed line at zero, and returns a display object whose scatter_ attribute holds the drawn points, the same kind of object ax.scatter returned in the Theory. For each panel, print its title, how many rows it was given and how many points it drew, and make each panel draw every row it was given.

The second is a table of the test listings per neighborhood, built with 2.6.2's named aggregation and sorted by mean residual: the number of listings, the median price, the mean residual, the MAE and the MAPE, with the mean residual and the MAE in whole units and the MAPE to three places. Below it, answer in comments which neighborhoods the model prices too high and which too low, which neighborhood's sellers it misleads most in proportion to their price, and why that is not the neighborhood with the largest MAE.

Hints:

• Count each panel's points before you trust the figure, the way the Theory counted its scatter. If a panel drew fewer points than the rows it was given, help(PredictionErrorDisplay.from_predictions) names the parameter that decides how many it draws.
• Named aggregation reads columns that already exist, so put each test listing's absolute residual, and that residual's share of its price, on the test rows' frame first, then average each per neighborhood.

import matplotlib.pyplot as plt
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import PredictionErrorDisplay
from sklearn.model_selection import train_test_split

data = pd.read_csv("listings.csv").dropna(subset=["area_sqm"])
X = data[["area_sqm"]]
y = data["price"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
model = LinearRegression().fit(X_train, y_train)
train_pred = model.predict(X_train)
test_pred = model.predict(X_test)

fig, ax = plt.subplots()
shown = PredictionErrorDisplay.from_predictions(
    y_train, train_pred, kind="residual_vs_predicted", ax=ax
)
# 2300 1000
print(len(y_train), len(shown.scatter_.get_offsets()))
plt.close(fig)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
panels = [
    (y_train, train_pred, "training rows"),
    (y_test, test_pred, "test rows"),
]
for ax, (truth, pred, title) in zip(axes, panels):
    shown = PredictionErrorDisplay.from_predictions(
        truth, pred, kind="residual_vs_predicted", subsample=None, ax=ax
    )
    ax.set_title(title)
    drawn = len(shown.scatter_.get_offsets())
    print(title, len(truth), drawn)
fig.savefig("residual_report.png", dpi=150, bbox_inches="tight")
plt.close(fig)

training rows 2300 2300
test rows 575 575

check = data.loc[X_test.index].assign(residual=y_test - test_pred)
check["abs_residual"] = check["residual"].abs()
check["abs_share"] = check["abs_residual"] / check["price"]
table = check.groupby("neighborhood").agg(
    listings=("residual", "size"),
    median_price=("price", "median"),
    mean_residual=("residual", "mean"),
    mae=("abs_residual", "mean"),
    mape=("abs_share", "mean"),
).sort_values("mean_residual")
print(table.round({"mean_residual": 0, "mae": 0, "mape": 3}))

# too high, a mean residual below zero: Hillcrest, Northend, Southgate
# and Downtown. Too low: Riverside, Lakeview and Eastgate.
# Hillcrest's sellers are misled most in proportion, 0.412: misses of
# 120,665 on prices near 298,000. Eastgate has the largest MAE,
# 147,786, but on prices near 533,000 that is a share of 0.257.

              listings  median_price  mean_residual       mae   mape
neighborhood                                                        
Hillcrest           73      298000.0      -120177.0  120665.0  0.412
Northend            67      292000.0      -107643.0  121490.0  0.403
Southgate           96      338000.0       -54403.0   66476.0  0.221
Downtown           125      274000.0       -13549.0   39068.0  0.155
Riverside           80      358000.0        28789.0   56984.0  0.139
Lakeview            67      480000.0       109482.0  117885.0  0.203
Eastgate            67      533000.0       145253.0  147786.0  0.257

Additional Resources

Optional, and nothing later depends on either. Worth opening for four data sets a summary cannot tell apart, or for why a least squares fit's R² on its training rows never falls when a feature is added.

• Wikipedia - "Anscombe's quartet". Four small data sets with nearly identical means, variances, correlation, fitted line and R², which look nothing alike once drawn.
• Wikipedia - "Coefficient of determination", sections "Inflation of R²" and "Adjusted R²". Why a least squares fit's R² on the rows it was fitted to never falls when a feature is added, and the adjusted version that charges for each one.

Summary

A model's record on a set of listings is its residuals, and every metric in this lesson keeps one property of that record and throws the rest away. MAE keeps the typical size of a miss, RMSE that size with the large misses weighted up, MAPE their size against each price, and R² their size against what a mean would have left. None of them is the record, which is why the plot, the one view that keeps every residual, comes before any of them is reported.

What gives a number its meaning is the comparison behind it. A training mean of zero compared with nothing says only that an intercept was there, an R² set beside an R² from other rows compares the rows, and a metric picked after the models were scored compares nothing but the hopes of whoever picked it. The comparisons that hold are made on the same rows, against the same baseline, with the metric fixed by the question before the numbers arrive.

The residuals are also where the next step comes from. The fan asks for the price on a log scale, the neighborhood means ask for a feature the model was never given, and a test R² of 0.565 showed neither. From here on, every feature a model gains can be judged the way area was judged here, by the residuals it leaves on rows the fit never saw.

Cheat Sheet

Part A - Reference Table

Part B - Compact Code Block

from sklearn.metrics import (mean_absolute_error, r2_score,
                             root_mean_squared_error,
                             mean_absolute_percentage_error)

res = y_test - pred                        # price minus prediction
mae = mean_absolute_error(y_test, pred)    # targets first
rmse = root_mean_squared_error(y_test, pred)
r2 = r2_score(y_test, pred)                # vs the rows' own mean
mape = mean_absolute_percentage_error(y_test, pred)
print(f"{mape:.1%}")                       # a fraction, shown as %

fig, ax = plt.subplots(figsize=(8, 5))
points = ax.scatter(pred, res, alpha=0.4)  # prediction on x
ax.axhline(0, color="black", linestyle="--")
print(len(points.get_offsets()))           # every row drawn?