Regression and Ranking Metrics

Regression metrics describe the sizes of prediction errors; choosing one for model selection favors some kinds of errors over others. Ranking metrics evaluate the order of items and, often, the positions where users encounter them. The choice depends on how the predictions will be used.

MAE and RMSE ask for different statistics

For errors \(e_i=y_i-\hat y_i\), mean absolute error is \(MAE=\frac1n\sum_i|e_i|\), and root mean squared error is \(RMSE=\sqrt{\frac1n\sum_i e_i^2}\). Errors of 1 and 3 give MAE 2 and RMSE \(\sqrt5\approx2.236\), both in the target’s units. Squaring gives the larger error more weight. The square root is increasing, so minimizing RMSE on one fixed evaluation set selects the same predictions as minimizing MSE. Averaging separate fold RMSEs is a different aggregation and need not preserve that equivalence.

The loss-functions article showed that the constant minimizing squared error is the mean and the constant minimizing absolute error is a median. For population prediction at a given input, the corresponding targets are the conditional mean and a conditional median, with finite conditional second moment for squared error and finite conditional absolute expectation for absolute error. A fitted model estimates these targets within the limits of its model family and training data.

import numpy as np
from scipy.optimize import minimize_scalar

rng = np.random.default_rng(0)
y = rng.lognormal(0.0, 1.0, 200_000)          # right-skewed, mean far above median

mae = minimize_scalar(lambda c: np.abs(y - c).mean()).x
rmse = minimize_scalar(lambda c: np.sqrt(((y - c) ** 2).mean())).x
print(f"MAE  minimizer {mae:.4f}   empirical median {np.median(y):.4f}")
print(f"RMSE minimizer {rmse:.4f}   empirical mean   {y.mean():.4f}")
# MAE  minimizer 0.9998   empirical median 0.9998
# RMSE minimizer 1.6506   empirical mean   1.6506

In this lognormal example, the RMSE-minimizing constant is about 65% higher than the MAE-minimizing constant. Large values pull the mean upward; the median depends on their order rather than how far above it they lie. An expected total calls for means, while a median describes the halfway point of outcomes. Inventory decisions additionally depend on the costs of shortages and excess stock, which can favor a different quantile.

Pinball loss at quantile \(\tau\) has a conditional \(\tau\)-quantile as its population optimum. At \(\tau=0.5\), it equals half the absolute loss; squared loss is a separate loss targeting the mean. Quantile regression can estimate lower and upper prediction limits without specifying a full parametric distribution. The resulting interval still needs coverage checks on new data.

Percentage errors are not symmetric

MAPE is mean absolute percentage error: average \(|y-\hat y|/|y|\) over observations. The table shows each observation’s contribution as a fraction, so 0.5 means 50%. It also shows one common symmetric MAPE (sMAPE) formula, \(2|y-\hat y|/(|y|+|\hat y|)\). Its symmetry refers to exchanging actual and predicted values, not to equal-sized errors above and below a fixed actual.

actual_pred = [(100, 50), (100, 150), (100, 200), (1, 2), (0.01, 0.02)]
print(f"{'actual':>8} {'pred':>8} {'abs err':>9} {'MAPE':>8} {'sMAPE':>8}")
for a, p in actual_pred:
    mape = abs(a - p) / abs(a)
    smape = abs(a - p) / ((abs(a) + abs(p)) / 2)
    print(f"{a:8.2f} {p:8.2f} {abs(a - p):9.2f} {mape:8.4f} {smape:8.4f}")
#   actual     pred   abs err     MAPE    sMAPE
#   100.00    50.00     50.00   0.5000   0.6667
#   100.00   150.00     50.00   0.5000   0.4000
#   100.00   200.00    100.00   1.0000   0.6667
#     1.00     2.00      1.00   1.0000   0.6667
#     0.01     0.02      0.01   1.0000   0.6667

At a fixed nonzero actual, MAPE gives equal penalties to equal absolute errors above and below it, as the first two rows show. For positive actuals and nonnegative predictions, under-prediction has MAPE at most 1.0, attained at a prediction of zero; over-prediction has no such bound. MAPE also weights each absolute error by the reciprocal of its actual value, giving small actuals more influence. This weighting can favor low forecasts. For actuals 1 and 10, a constant forecast of 1 gives MAPE \((0+0.9)/2=0.45\); forecasting 10 gives \((9+0)/2=4.5\). Each forecast has the same MAE, 4.5, but MAPE weights the error on the smaller actual more heavily. In the first two rows, sMAPE assigns different penalties to equal absolute errors: 0.6667 below the actual and 0.4000 above it. This particular comparison does not establish a universal upward bias for models optimized with sMAPE.

The last two rows show scale invariance: doubling 0.01 receives the same percentage error as doubling 1.00. That is useful when proportional error matters, but it may not match costs measured in absolute units. MAPE is undefined at a zero actual, and the displayed sMAPE formula is undefined when both actual and prediction are zero. With strictly positive actuals and predictions, MAE on their logarithms measures multiplicative discrepancies. For forecasting, MASE uses a common scale: divide forecast MAE by the mean absolute naive error computed on the training series. The naive forecast is usually the previous observation, or the observation one season earlier for seasonal data; a zero scaling denominator makes MASE undefined.

R-squared depends on the variation in the target

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_absolute_error

rng = np.random.default_rng(0)
n = 20_000
print(f"{'feature sd':>11} {'noise sd':>9} {'R^2':>8} {'MAE':>8} {'RMSE':>8}")
for sd in (1.0, 5.0, 50.0):
    X = rng.normal(0, sd, (n, 1))
    y = 2.0 * X[:, 0] + rng.normal(0, 3.0, n)      # noise sd is 3.0 in every row
    pred = LinearRegression().fit(X, y).predict(X)
    print(f"{sd:11.1f} {3.0:9.1f} {r2_score(y, pred):8.4f}"
          f" {mean_absolute_error(y, pred):8.4f}"
          f" {np.sqrt(((y - pred) ** 2).mean()):8.4f}")
#  feature sd  noise sd      R^2      MAE     RMSE
#         1.0       3.0   0.3049   2.4074   3.0229
#         5.0       3.0   0.9173   2.3870   2.9937
#        50.0       3.0   0.9991   2.4045   3.0179

Each row generates fresh data, fits a new linear model, and evaluates it on that same training sample. The noise distribution stays fixed while the spread of the feature changes. RMSE remains near the noise standard deviation of 3. For centered normal errors, MAE is \(\sigma\sqrt{2/\pi}\), about 2.394 here, consistent with the printed values. In \(R^2=1-\mathrm{SSE}/\mathrm{SST}\), SSE is the sum of squared prediction errors and SST is the sum of squared deviations from the sample mean. Spreading the feature increases SST substantially while SSE stays at roughly the same scale, so \(R^2\) rises from 0.3049 to 0.9991.

On a fixed evaluation set with a nonconstant target, \(R^2\) compares squared error with predicting that set’s mean. It can be negative when the model is worse than that reference. For a constant target, SST is zero, so the displayed ratio needs a separate convention. The evaluation-set mean is a reference used to define the metric, not a deployable forecast learned from training data. Across datasets, this reference changes with target variation, so a higher \(R^2\) alone does not establish smaller prediction errors. Report an error in the units of the target alongside it, and use held-out data when assessing performance on new observations.

Ranking metrics weight positions

For an ordered list, we now evaluate which items appear first. In search or recommendation, early positions often matter most. The example below represents one query with 20 candidate items: rel holds their true relevance, and the score array determines their order. Wrapping each array in a list gives scikit-learn one row per query. Its NDCG (normalized discounted cumulative gain) uses relevance itself as the gain, discounted at rank \(r\) by \(1/\log_2(r+1)\), and divides the total by the best possible total at the same cutoff. For two items with relevance 3 and 1, the ideal total is \(3+1/\log_2 3\approx3.631\). Reversing them gives \(1+3/\log_2 3\approx2.893\), or NDCG@2 of about 0.797. Some other implementations use \(2^{\mathrm{relevance}}-1\) as gain, so check that convention before comparing scores. Use nonnegative relevance judgments; the usual 0-to-1 interpretation also requires a positive ideal total. A query with no relevant items needs an explicit handling rule when averaging scores across queries.

import numpy as np
from scipy.stats import kendalltau
from sklearn.metrics import ndcg_score

n = 20
rel = np.zeros(n); rel[0] = 3; rel[1] = 2; rel[2] = 2; rel[3] = 1; rel[4] = 1
base = np.arange(n, 0, -1).astype(float)          # a perfect ranking

def swapped(i, j):
    s = base.copy(); s[i], s[j] = s[j], s[i]; return s

for label, s in (("perfect", base),
                 ("swap ranks 1-2", swapped(0, 1)),
                 ("swap ranks 19-20", swapped(18, 19))):
    print(f"{label:18s} Kendall tau {kendalltau(base, s).statistic:.4f}"
          f"  NDCG@5 {ndcg_score([rel], [s], k=5):.4f}")
# perfect            Kendall tau 1.0000  NDCG@5 1.0000
# swap ranks 1-2     Kendall tau 0.9895  NDCG@5 0.9393
# swap ranks 19-20   Kendall tau 0.9895  NDCG@5 1.0000

Here Kendall tau compares each score ordering with the base ordering. Both changes swap adjacent items, reversing one pair and producing 0.9895; exchanging more distant items would reverse more pairs. NDCG@5 falls from 1 to 0.9393 after the top swap. The bottom swap is outside the cutoff and also exchanges two items with zero relevance, so it would leave even full-list NDCG unchanged. This example illustrates different evaluation criteria without isolating position discounting alone. For an interface displaying five results, NDCG@5 measures the relevance ordering in those five positions relative to the ideal ordering of the supplied candidates. Relevant items missing from the candidate set are not recovered by this calculation, so candidate retrieval may need its own evaluation.

MRR averages the reciprocal rank of the first relevant item across queries and suits tasks where finding one useful answer is the goal. Untruncated MRR has no cutoff. Recall@k divides the number of relevant items retrieved in the top \(k\) by the total number of relevant items, so retrieving only half gives 0.5. NDCG@k can reach 1 when the top \(k\) is ideal even if more relevant items remain below it. MAP@k averages truncated AP across queries; implementations differ on whether its normalization uses all relevant items or at most \(k\), so state the convention. For cutoff-based metrics, choose \(k\) to reflect the number of results users can inspect.

MAE against RMSE on real data

The California housing data records median house values for 20,640 census block groups, in units of $100,000. Each row describes a block group, not an individual home, so the dollar errors below concern predicted block-group medians. The recorded target is capped, which matters when interpreting errors near its upper limit. The first dataset fetch requires internet access; later calls can use the local cache.

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error

X, y = fetch_california_housing(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0)
pred = HistGradientBoostingRegressor(random_state=0).fit(Xtr, ytr).predict(Xte)

err = yte - pred
mae = mean_absolute_error(yte, pred)
rmse = mean_squared_error(yte, pred) ** 0.5
print(f"MAE  {mae:.4f}   RMSE {rmse:.4f}   ratio {rmse / mae:.3f}")
cap = y.max()
print(f"recorded cap {cap:.5f}: {(yte == cap).mean():.4f} of test rows equal the cap")

big = np.abs(err) > 3 * np.std(err)
print(f"{big.mean() * 100:.2f}% of rows are >3sd errors and contribute "
      f"{(err[big] ** 2).sum() / (err ** 2).sum() * 100:.1f}% of the squared error")
# MAE  0.3175   RMSE 0.4772   ratio 1.503
# recorded cap 5.00001: 0.0464 of test rows equal the cap
# 1.84% of rows are >3sd errors and contribute 34.4% of the squared error

The MAE of 0.3175 represents about $31,750, and RMSE of 0.4772 about $47,720. Their ratio is 1.503. Centered normal errors have a population ratio of about 1.253, but the ratio alone does not identify a tail distribution. Inspecting the individual squared errors shows how concentrated their contribution is.

The code flags rows whose absolute error exceeds three times the residual standard deviation. These are 1.84% of the test rows and account for 34.4% of total test squared error. This describes the evaluated errors; it does not measure those rows’ contribution to the training objective. Whether to emphasize such large misses depends on their cost in the application.

The middle line counts 287 test targets, or 4.64%, exactly equal to the recorded maximum of 5.00001. These are capped measurements: their uncapped values are not available. If the goal is to predict the recorded target, evaluating against that target is appropriate. If the goal is the uncapped price, errors against capped labels cannot establish accuracy above the ceiling. The target definition determines what these scores can tell us.

Exercises

1. Which metric picks which model. Fit models minimizing squared and absolute error on skewed data and evaluate both under both metrics. OLS means ordinary least squares; LAD means least absolute deviations. Here both fit a straight line, and QuantileRegressor(quantile=0.5, alpha=0) supplies the unregularized LAD fit.

You should get: lower test RMSE for OLS and lower test MAE for LAD in the supplied simulation. Explain why these observed rankings are not guarantees for other test samples.

import numpy as np
from sklearn.linear_model import LinearRegression, QuantileRegressor
from sklearn.metrics import mean_absolute_error

rng = np.random.default_rng(0)
n = 4000
X = rng.uniform(0, 10, (n, 1))
# heavy right tail whose size grows with x
y = 2 * X[:, 0] + rng.lognormal(0.0, 1.2, n) * (1 + X[:, 0] / 5)

Xtr, Xte, ytr, yte = X[:3000], X[3000:], y[:3000], y[3000:]
ols = LinearRegression().fit(Xtr, ytr)
lad = QuantileRegressor(quantile=0.5, alpha=0.0, solver="highs").fit(Xtr, ytr)

for name, m in (("squared error (OLS)", ols), ("absolute error (LAD)", lad)):
    p = m.predict(Xte)
    print(f"{name:22s} RMSE {np.sqrt(((yte - p) ** 2).mean()):7.4f}"
          f"  MAE {mean_absolute_error(yte, p):7.4f}"
          f"  median residual {np.median(yte - p):+7.4f}")
# squared error (OLS)    RMSE  6.3193  MAE  3.6638  median residual -1.8627
# absolute error (LAD)   RMSE  6.6869  MAE  3.1630  median residual +0.0106

2. What NDCG@k cannot see. Move five relevant items out of the top ten and compare the original and modified rankings at cutoffs 10 and 50.

You should get: a larger decrease in NDCG@10 than in NDCG@50 for this example, and an explanation of how the cutoff changes the credit given to demoted items.

import numpy as np
from sklearn.metrics import ndcg_score

rel = np.zeros(50); rel[:10] = 3                  # ten relevant items

good = np.arange(50, 0, -1).astype(float)         # all ten ranked at the top
mixed = good.copy()
mixed[5:10] = 0.0                                 # five of them pushed to the bottom
mixed[45:50] = 45.0                               # five irrelevant items promoted

for label, s in (("all relevant in top 10", good), ("5 pushed to the bottom", mixed)):
    print(f"{label:26s} NDCG@10 {ndcg_score([rel], [s], k=10):.4f}"
          f"  NDCG@50 {ndcg_score([rel], [s]):.4f}")
# all relevant in top 10     NDCG@10 1.0000  NDCG@50 1.0000
# 5 pushed to the bottom     NDCG@10 0.6489  NDCG@50 0.8450

3. A high R-squared that does not beat a baseline. Build a dataset where \(R^{2}\) is above 0.9 and the model is worse than a trivial baseline on the question being asked.

You should get: a high \(R^{2}\) alongside a naive baseline that beats the model on the metric that matters.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_absolute_error

rng = np.random.default_rng(0)
n = 3000
t = np.arange(n)
level = 100 + np.cumsum(rng.normal(0, 1, n))       # a random walk
y = level + rng.normal(0, 2, n)

# "predict today's value from yesterday's" — one-step prediction using an observed lag
X = y[:-1].reshape(-1, 1)
target = y[1:]
Xtr, Xte = X[:2400], X[2400:]
ytr, yte = target[:2400], target[2400:]

pred = LinearRegression().fit(Xtr, ytr).predict(Xte)
naive = Xte[:, 0]                                   # predict yesterday's value

print(f"model  R^2 {r2_score(yte, pred):7.4f}  MAE {mean_absolute_error(yte, pred):7.4f}")
print(f"naive  R^2 {r2_score(yte, naive):7.4f}  MAE {mean_absolute_error(yte, naive):7.4f}")
ratio = mean_absolute_error(yte, pred) / mean_absolute_error(yte, naive)
print(f"MAE relative to the test-set naive baseline {ratio:.4f}")
# model  R^2  0.9392  MAE  2.4236
# naive  R^2  0.9411  MAE  2.3863
# MAE relative to the test-set naive baseline 1.0156