Judging a regressor: error and R²

What this chapter covers

The last chapter fit a line and quoted an R² of 0.585 like it meant something. This chapter is where we earn that number, and the five others you'll reach for alongside it. A regression model hands you a column of predictions; the truth is another column next to it; every metric here is a different way to squeeze those two columns down to one number that says how wrong the model is. The catch, and the whole reason this is its own chapter, is that those numbers are not interchangeable — they live in different units, they reward different things, and one of them will lie to you the moment an outlier shows up.

We build all six from scratch in NumPy: MSE, its square root RMSE, MAE, R², adjusted R², and MAPE. Then we check every one against sklearn.metrics and they match to the last decimal. The running question the whole way through is the one nobody asks often enough: what unit is this number in? MSE is in the target's units squared, which is why you can't read it. RMSE and MAE are back in the target's own units, so you can. R² is unitless on purpose. MAPE is a percent. Keep that question in your head and half of these metrics stop being interchangeable and start being tools with different jobs.

The predictions we judge are the ones from the linear-regression chapter: the same California housing model, the same 80/20 split, its guesses on the 4,128 held-out blocks. The target is median house value in units of $100,000, so every dollar figure below traces back to that unit.

A bit of history

Squared error is the oldest of these by two centuries, and it arrived as an optimization criterion long before anyone thought of it as a score. Adrien-Marie Legendre published least squares in 1805; Carl Friedrich Gauss published his own version in 1809 and tied it to the normal distribution, arguing that minimizing the sum of squared residuals is the right thing to do when your noise is Gaussian. That sum of squares, divided by n, is the MSE. So MSE didn't start life as a report-card number — it started as the thing you minimize, and it's still the best of the bunch at that job. RMSE and MAE are just the two obvious ways to make an error you can actually interpret: take the root, or drop the squaring entirely and use absolute value, an idea that predates Gauss and goes back to Roger Boscovich in the 1750s.

R² is much younger and comes from genetics, not astronomy. Sewall Wright, working out path analysis to untangle how much of the variation in a trait traced to which cause, introduced the coefficient of determination in 1921 — the fraction of variance one thing accounts for in another. It landed in the same decade that Ronald Fisher was building the analysis of variance, and the two ideas are the same idea: split the total variation into the part your model explains and the part it doesn't, and report the ratio. That's why R² is unitless and why it's the number people quote across completely different problems — it was designed from the start as a proportion, not a measurement.

The intuition

Every metric on this page starts from the same picture: the prediction and the truth, plotted against each other. If the model were perfect every point would sit on the 45-degree line, where predicted equals actual. It isn't, so the points scatter around that line, and the vertical distance from each point to it is that prediction's error — its residual. Every metric is a different rule for collapsing that cloud of vertical gaps into a single summary.

Here are the model's predictions on the held-out California blocks, actual on the x-axis and predicted on the y, with the y = x line where a perfect model would put everything. Both axes are in units of $100,000.

Two things jump out, and both are things a metric will have to summarize. The cloud is tilted along the line, which is good — high predictions mostly go with high actuals. But look at the right edge: a vertical wall of points at actual = 5.0 where the model predicts all over the place, because the target is capped at $500,000 and the model, not knowing about the cap, keeps guessing higher. Those are big residuals. A metric that squares them will fixate on them; a metric that doesn't will shrug. Same cloud, different verdicts, and that difference is the whole chapter.

The math

Write $y_i$ for the true value of block $i$, $\hat{y}_i$ for the model's prediction, $\bar{y}$ for the mean of the true values, and $n$ for the number of blocks. Every metric is a function of the residuals $y_i - \hat{y}_i$.

Mean squared error averages the squared residuals. Squaring is doing two jobs: it kills the sign so misses can't cancel, and it makes a big miss count far more than a small one.

Because it squares, its unit is the target's unit squared — dollars-squared, an area. Root mean squared error undoes that, dropping you back into the target's own units:

Mean absolute error takes the other road: don't square at all, just take the absolute value. Same units as RMSE, but every residual counts in proportion to its size, no more.

R² is a ratio of two sums of squares. The residual sum of squares is your model's total squared error; the total sum of squares is the error you'd make predicting the mean $\bar{y}$ for every block:

Both sums are in the target's units squared, so the ratio is unitless — that's the point of it. R² of 1 is perfect, 0 is no better than guessing the mean, and negative means you'd have been better off with the mean. Adjusted R² is the same number with a penalty for the number of features $p$, so adding a useless column can't inflate it:

Mean absolute percentage error divides each residual by the truth before averaging, which is what makes it scale-free — and what makes it explode if any $y_i$ is near zero:

What each one rewards, and what it punishes

The split that matters most is squared versus absolute. MSE and RMSE square the residual, so they care about big misses out of all proportion — a residual twice as large counts four times as much. That's exactly what you want when you're training a model, because it makes the fit work hardest where it's most wrong, and it's exactly what you don't want when you're reporting, because one outlier can dominate the whole number. MAE weights linearly, so it reports the typical miss and stays calm when a few points go wild. Neither is more correct; they answer different questions. "How bad is my worst behavior" is RMSE. "How wrong am I on a normal day" is MAE.

R² and MAPE are the two that let you compare across problems, and both pay for it. R² is scale-free because it's measured against the variance of the target itself, which is a strength — it tells you how much better than the dumb mean-guess you're doing — and a trap, because a high R² on a target that barely varies is easy and a low R² on a wildly variable one might be great. MAPE is scale-free because it's a percentage, which makes it the one metric you can say out loud in a meeting, but it's biased: it blows up near zero targets and it penalizes over-prediction less than under-prediction, so a model that systematically guesses low can score suspiciously well. The honest move is to never report just one.

The predictions we're judging

The numbers below all come from one file: the linear model's predictions on the 4,128 held-out California blocks, saved next to their true values in test_predictions.csv. The model was fit on the other 16,512 blocks — same OLS fit, same seed-0 split as the linear-regression chapter — and it uses eight features, which is the $p$ that adjusted R² will charge us for. The target is median house value in units of $100,000, and it ranges from 0.15 to the capped 5.0, which matters twice over: there are no zeros, so MAPE is well-defined, and there's that ceiling, which is where the fat residuals come from.

Loading them is the only I/O in the from-scratch file — every metric takes it from here:

def mse(y_true, y_pred):
    """Mean squared error. UNIT: target units, SQUARED.

    Average of the squared residuals. If the target is in $100k, this number is
    in ($100k)^2 — an area, not a price, which is exactly why you can't read it
    off as "the model is off by X." Squaring makes every miss positive and makes
    a big miss count far more than a small one, so a single outlier moves this a
    lot. Great to optimize, hard to interpret.
    """
    resid = y_true - y_pred
    return float(np.mean(resid ** 2))

That's MSE, and it's worth saying what its number will be before we even run it: whatever comes out is in $100k-squared, an area, and you should not try to picture it. It's a great thing to minimize and a terrible thing to interpret, which is the recurring theme.

Build them, one metric at a time

RMSE is MSE with the interpretation bolted back on. One square root and the number is in the target's units again — a figure you can hold next to a house price.

def rmse(y_true, y_pred):
    """Root mean squared error. UNIT: target units.

    The square root of the MSE, which undoes the squaring and lands you back in
    the target's own units — dollars, not dollars-squared. That's the whole
    reason it exists: it's the MSE made readable. Still outlier-sensitive (the
    squaring happened before the root), but now it's a number you can compare to
    the scale of what you're predicting.
    """
    return float(np.sqrt(mse(y_true, y_pred)))

MAE is the other summary of the residuals, and the robust one. Absolute value instead of square, so the units are the target's from the start and no single point can lever the total:

def mae(y_true, y_pred):
    """Mean absolute error. UNIT: target units.

    Average of the absolute residuals. Same units as RMSE — the target's units —
    but it weights every miss linearly instead of squaring it, so one wild point
    contributes its size and no more. This is the robust one: the typical error,
    unimpressed by outliers.
    """
    return float(np.mean(np.abs(y_true - y_pred)))

R² compares your squared error against the squared error of the mean-guess. Note the two sums the math defined — ss_res is the model's, ss_tot is the mean's — and that the ratio comes out unitless because both are in the same squared units:

def r2(y_true, y_pred):
    """Coefficient of determination. UNIT: none (a fraction, <= 1).

    R^2 = 1 - SS_res / SS_tot, where
      SS_res = sum of squared residuals (our model's squared error), and
      SS_tot = sum of squared deviations from the mean of y (the error you'd get
               predicting the mean for everything).
    So it's the fraction of the target's variance the model explains. 1.0 is
    perfect, 0.0 is no better than always guessing the mean, negative is worse
    than the mean. Unitless by construction — a ratio of two things in the same
    (squared) units — which is what lets you compare it across datasets.
    """
    ss_res = float(np.sum((y_true - y_pred) ** 2))
    ss_tot = float(np.sum((y_true - np.mean(y_true)) ** 2))
    return 1.0 - ss_res / ss_tot

Adjusted R² wraps that with the feature penalty. It needs to know how many predictors went into the model, because the whole point is to make each one pay rent:

def adjusted_r2(y_true, y_pred, n_features):
    """Adjusted R^2. UNIT: none.

    Plain R^2 never goes down when you add a feature, even a useless one, so it
    rewards fatter models. Adjusted R^2 charges rent for each feature: it
    discounts R^2 by the number of predictors `p` relative to the sample size
    `n`, and only rises when a new feature earns more than it costs.

        1 - (1 - R^2) * (n - 1) / (n - p - 1)
    """
    n = len(y_true)
    r2_val = r2(y_true, y_pred)
    return 1.0 - (1.0 - r2_val) * (n - 1) / (n - n_features - 1)

And MAPE, the percentage. The division by y_true is what makes it scale-free and also what makes it fragile — it's the one line in this file that can divide by zero, and the only reason it doesn't here is that house values are never zero:

def mape(y_true, y_pred):
    """Mean absolute percentage error. UNIT: percent.

    Average of |residual / true value|, times 100. Because it divides by the
    truth, it's scale-free — a $5k miss on a $50k house (10%) counts the same as
    a $50k miss on a $500k house (10%). That makes it easy to explain to a
    non-technical audience, but it has two teeth: it blows up when any true value
    is near zero, and it punishes over-prediction less than under-prediction.
    Undefined if any `y_true` is exactly zero (California's target never is).
    """
    return float(np.mean(np.abs((y_true - y_pred) / y_true)) * 100.0)

Run all six on the test-set predictions and here is what the model actually scores, every bar in its own unit because they don't share one. Read each with the label — a bar of MSE and a bar of MAE next to each other is a category error, and the chart is drawn to make you feel that:

The four numbers to carry forward, each with its unit: RMSE is 0.727 $100k — call it a typical error of about $72,700. MAE is 0.529 $100k, roughly $52,900, noticeably smaller than RMSE, and that gap is itself information: RMSE exceeds MAE exactly when the errors are uneven, when a minority of big misses (the capped blocks) are pulling the squared metric up. MSE is 0.528 $100k-squared, a number I genuinely can't interpret and won't try to. R² is 0.585, unitless — the model explains about 59% of the variance in house value, and adjusted R² barely moves it to 0.584 because eight features against 4,128 test points is a rounding error of a penalty. MAPE is 32.06%, meaning the model is off by about a third of the true value on a typical block, which sounds worse than the R² does and is a fair reminder that "explains 59% of variance" and "off by a third" are both true at once.

Watch it work

This is the demonstration that makes the squared-versus-absolute split stop being abstract. Twenty toy points, a fixed fit line drawn through them, and one point — the orange one — that we drag straight up, frame by frame, so its residual grows from nothing to twelve units. The predictions never change; only that one true value moves. Underneath, three bars: MSE, RMSE, and MAE, recomputed over all twenty points at every frame, sharing one axis so you can watch them diverge.

Press play. The dashed orange line is the growing residual of the point we're dragging.

Watch what the bars do. At the start, with the point sitting near the line, the three metrics are close together — MSE 0.17 units², RMSE 0.42 units, MAE 0.35 units, all in the same rough neighborhood. Then drag the one point out to a residual of twelve units and they come apart violently. MSE runs to 7.72 units² — it goes up by a factor of forty-five off the strength of a single point, because that residual gets squared before it's averaged in. RMSE climbs to 2.78 units, a factor of about seven, tracking the square of the outlier but softened by the root. And MAE crawls from 0.35 to 0.95 units, less than tripling, because to a linear average one point twelve units off is just one point twelve units off. One outlier, three completely different stories.

This is the whole practical lesson of the chapter in one animation. If your data has outliers you can't clean — a capped target, a fat-tailed distribution, a few data-entry disasters — MSE and RMSE will report your model as dramatically worse than it is on the bulk of the data, because they're answering "how bad is the worst of this," and the worst is an outlier. MAE will keep telling you about the typical block. Which one you want depends entirely on whether those big misses are the thing you care about or the thing you want to ignore, and that's a judgment call the metric can't make for you — but you can't even make it if you only ever look at one number.

The full implementation

The whole file, no library, top to bottom — six metrics and one loader, each in its stated unit. This is the code every number above and every bar in the animation actually ran:

"""Regression metrics, built from scratch in pure NumPy.

Every function here answers the same question — "how wrong is this regressor?" —
and every one answers it in a different unit. That unit is the whole point of
the chapter, so each docstring states it plainly.

The inputs are always the same: `y_true`, the real targets, and `y_pred`, the
model's predictions, both 1-D arrays of the same length. None of these functions
know or care where the predictions came from; they judge a column of numbers
against another column of numbers. We feed them the test-set predictions of the
linear model from the linear-regression chapter (California housing, target in
units of $100,000).

Pure NumPy — no ML library in this file. The `# region:` markers are what the
book's include directives pull in, one metric at a time.
"""

import numpy as np
import pandas as pd


# region: mse
def mse(y_true, y_pred):
    """Mean squared error. UNIT: target units, SQUARED.

    Average of the squared residuals. If the target is in $100k, this number is
    in ($100k)^2 — an area, not a price, which is exactly why you can't read it
    off as "the model is off by X." Squaring makes every miss positive and makes
    a big miss count far more than a small one, so a single outlier moves this a
    lot. Great to optimize, hard to interpret.
    """
    resid = y_true - y_pred
    return float(np.mean(resid ** 2))
# endregion


# region: rmse
def rmse(y_true, y_pred):
    """Root mean squared error. UNIT: target units.

    The square root of the MSE, which undoes the squaring and lands you back in
    the target's own units — dollars, not dollars-squared. That's the whole
    reason it exists: it's the MSE made readable. Still outlier-sensitive (the
    squaring happened before the root), but now it's a number you can compare to
    the scale of what you're predicting.
    """
    return float(np.sqrt(mse(y_true, y_pred)))
# endregion


# region: mae
def mae(y_true, y_pred):
    """Mean absolute error. UNIT: target units.

    Average of the absolute residuals. Same units as RMSE — the target's units —
    but it weights every miss linearly instead of squaring it, so one wild point
    contributes its size and no more. This is the robust one: the typical error,
    unimpressed by outliers.
    """
    return float(np.mean(np.abs(y_true - y_pred)))
# endregion


# region: r2
def r2(y_true, y_pred):
    """Coefficient of determination. UNIT: none (a fraction, <= 1).

    R^2 = 1 - SS_res / SS_tot, where
      SS_res = sum of squared residuals (our model's squared error), and
      SS_tot = sum of squared deviations from the mean of y (the error you'd get
               predicting the mean for everything).
    So it's the fraction of the target's variance the model explains. 1.0 is
    perfect, 0.0 is no better than always guessing the mean, negative is worse
    than the mean. Unitless by construction — a ratio of two things in the same
    (squared) units — which is what lets you compare it across datasets.
    """
    ss_res = float(np.sum((y_true - y_pred) ** 2))
    ss_tot = float(np.sum((y_true - np.mean(y_true)) ** 2))
    return 1.0 - ss_res / ss_tot
# endregion


# region: adjusted_r2
def adjusted_r2(y_true, y_pred, n_features):
    """Adjusted R^2. UNIT: none.

    Plain R^2 never goes down when you add a feature, even a useless one, so it
    rewards fatter models. Adjusted R^2 charges rent for each feature: it
    discounts R^2 by the number of predictors `p` relative to the sample size
    `n`, and only rises when a new feature earns more than it costs.

        1 - (1 - R^2) * (n - 1) / (n - p - 1)
    """
    n = len(y_true)
    r2_val = r2(y_true, y_pred)
    return 1.0 - (1.0 - r2_val) * (n - 1) / (n - n_features - 1)
# endregion


# region: mape
def mape(y_true, y_pred):
    """Mean absolute percentage error. UNIT: percent.

    Average of |residual / true value|, times 100. Because it divides by the
    truth, it's scale-free — a $5k miss on a $50k house (10%) counts the same as
    a $50k miss on a $500k house (10%). That makes it easy to explain to a
    non-technical audience, but it has two teeth: it blows up when any true value
    is near zero, and it punishes over-prediction less than under-prediction.
    Undefined if any `y_true` is exactly zero (California's target never is).
    """
    return float(np.mean(np.abs((y_true - y_pred) / y_true)) * 100.0)
# endregion


def load_predictions(path="../data/test_predictions.csv"):
    """Load the committed (y_true, y_pred) pair for the real test set.

    These are the linear model's predictions on the held-out California housing
    blocks, in units of $100,000. Returns (y_true, y_pred) as float arrays.
    """
    df = pd.read_csv(path)
    return df["y_true"].to_numpy(float), df["y_pred"].to_numpy(float)

The library version

You would never hand-roll these in practice; sklearn.metrics has all of them, and reaching for the library is the right call. The only reason to build them once is so you know exactly what the library is computing, because a couple of its choices will bite you if you don't. root_mean_squared_error is a separate function from mean_squared_error — sklearn used to fold it into a squared=False flag and then deprecated that, so on any recent version you call the dedicated one. And mean_absolute_percentage_error returns a fraction, not a percent: it gives you 0.32, and if you print it as "0.32% error" you're off by two orders of magnitude. We multiply by 100 to line it up with ours.

def sklearn_metrics(y_true, y_pred):
    """Every metric sklearn ships, as a dict, in the same units as impl.py.

    mape is scaled to a percent (sklearn returns a fraction) so it matches our
    `mape`. rmse uses sklearn's dedicated function rather than sqrt-ing the mse.
    """
    return {
        "mse": float(mean_squared_error(y_true, y_pred)),          # target units^2
        "rmse": float(root_mean_squared_error(y_true, y_pred)),    # target units
        "mae": float(mean_absolute_error(y_true, y_pred)),         # target units
        "r2": float(r2_score(y_true, y_pred)),                     # unitless
        "mape": float(mean_absolute_percentage_error(y_true, y_pred) * 100.0),  # percent
    }

There's no adjusted-R² helper in sklearn — that one stays in our file, because it needs to know the feature count, which is a property of your model, not of the two prediction columns. Everything else is a one-liner, and it's what you should ship.

Scratch versus library

Same predictions, our six formulas against sklearn's five, and they land on the same numbers. Each metric is drawn in its own row with its own scale — the two bars per metric are identical lengths, which is the only honest way to show "these agree" when the metrics themselves live in different units:

The agreement is exact — the largest difference between any of our numbers and sklearn's is zero to the precision the assertion checks, which is what you'd hope for from formulas this short. MSE 0.528 $100k-squared, RMSE 0.727 $100k, MAE 0.529 $100k, R² 0.585 unitless, MAPE 32.06%, matched line for line. The point of building them from scratch was never that sklearn might be wrong; it was to make the units impossible to lose track of. When you call mean_squared_error and get 0.528, you now know without thinking that it's an area, and that the number you actually want to quote is its root.

Takeaways

Report RMSE, and report it in the target's units. It's the one number that's both in the units you care about and sensitive to the misses that hurt, so "the model is off by about $72,700 on a typical block" is a sentence a stakeholder understands and a sentence that's true. RMSE is my default headline for a regressor, and the units are half of why — a metric you can't attach a unit to is a metric you can't defend in a room.

Keep MAE next to it, and read the gap. When RMSE sits well above MAE — 0.727 against 0.529 here — that spread is telling you the errors are uneven, that a minority of big misses is inflating the squared metric. If those big misses are outliers you don't care about, quote MAE and move on; if they're the capped $500,000 houses that your business does care about, that gap is the most important thing on the page. The outlier animation is the whole argument: one point moved MSE by 45x and MAE by less than 3x, and only by watching both did you learn which kind of wrong your model is.

Use R² for the scale-free question — how much better than guessing the mean — and never as your only number, because 0.585 and "off by a third of the value" are the same model. Reach for MAPE when a percentage is what the audience needs and you've checked there are no near-zero targets to blow it up. And MSE: minimize it, don't report it. It's the best training objective in this chapter and the worst report-card number, because it's the one figure here you cannot say what unit it's in without stopping to think — dollars squared — and a number you have to translate before you can believe it is a number that will eventually fool you. That's the running lesson: always know what unit you're holding.