Regression Metrics — MAE, RMSE, R²

A basketball commentator says "this team needs 12 points a quarter to win." The team scores 10, 11, 13, 9, 12, 8 — never exactly 12. MAE asks: how far off was each quarter on average? Answer: about 1.5 points. RMSE asks the same but doubles down on the 8-point quarter (4 under) — that squared miss from target hurts more than two smaller misses. MAPE asks: what percentage of the target was each miss?

Choose MAE when all errors cost equally — late by 5 minutes is 5× worse than late by 1 minute, nothing more. Choose RMSE when catastrophic errors cost disproportionately — one 45-minute delay is far worse than nine 5-minute delays.

import numpy as np
from sklearn.metrics import (mean_absolute_error,
                              mean_squared_error,
                              mean_absolute_percentage_error,
                              r2_score)

# ── Compute all four metrics from scratch ─────────────────────────────
y_true = np.array([41, 42, 29, 97, 35], dtype=float)
y_pred = np.array([32, 45, 28, 52, 38], dtype=float)

errors     = y_true - y_pred
abs_errors = np.abs(errors)

# MAE
mae_manual = abs_errors.mean()
mae_sklearn = mean_absolute_error(y_true, y_pred)

# MSE and RMSE
mse_manual  = (errors ** 2).mean()
rmse_manual = np.sqrt(mse_manual)
rmse_sklearn = np.sqrt(mean_squared_error(y_true, y_pred))

# MAPE
mape_manual  = (abs_errors / np.abs(y_true)).mean() * 100
mape_sklearn = mean_absolute_percentage_error(y_true, y_pred) * 100

# R²
ss_res = (errors ** 2).sum()
ss_tot = ((y_true - y_true.mean()) ** 2).sum()
r2_manual  = 1 - ss_res / ss_tot
r2_sklearn = r2_score(y_true, y_pred)

print("Manual vs sklearn verification:")
print(f"  MAE:   {mae_manual:.4f}  sklearn: {mae_sklearn:.4f}  match: {np.isclose(mae_manual, mae_sklearn)}")
print(f"  RMSE:  {rmse_manual:.4f}  sklearn: {rmse_sklearn:.4f}  match: {np.isclose(rmse_manual, rmse_sklearn)}")
print(f"  MAPE:  {mape_manual:.4f}%  sklearn: {mape_sklearn:.4f}%")
print(f"  R²:    {r2_manual:.4f}  sklearn: {r2_sklearn:.4f}  match: {np.isclose(r2_manual, r2_sklearn)}")

print(f"
Per-sample contribution to each metric:")
print(f"{'Error':>8} {'|error|':>8} {'error²':>8} {'% error':>8}")
print("─" * 38)
for e, ae, y in zip(errors, abs_errors, y_true):
    pct = abs(e) / abs(y) * 100
    print(f"  {e:>6.0f}  {ae:>8.0f}  {ae**2:>8.0f}  {pct:>7.1f}%")
print(f"  {'Mean':>6}  {mae_manual:>8.1f}  {mse_manual:>8.0f}  {mape_manual:>7.1f}%")
print(f"  {'':>6}  {'↑MAE':>8}  {'√→RMSE':>8}")

import numpy as np
from sklearn.metrics import r2_score
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 2000
distance = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic  = np.random.randint(1, 11, n).astype(float)
prep     = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
delivery = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
            + np.random.normal(0, 4, n)).clip(10, 120)

X = np.column_stack([distance, traffic, prep])
y = delivery
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

# ── R² for different models ───────────────────────────────────────────
models = {
    'Always predict mean (baseline)': None,
    'Linear Regression':              LinearRegression(),
    'Gradient Boosting':              GradientBoostingRegressor(
                                          n_estimators=200, random_state=42),
    'Shuffled labels (broken)':       'shuffled',
}

print(f"R² comparison — DoorDash delivery time:")
print(f"{'Model':<35} {'R²':>8} {'MAE':>8} {'RMSE':>8}")
print("─" * 64)

for name, model in models.items():
    if model is None:
        y_pred = np.full_like(y_te, y_tr.mean())
    elif model == 'shuffled':
        y_pred = np.random.permutation(y_te)   # random predictions
    else:
        model.fit(X_tr, y_tr)
        y_pred = model.predict(X_te)

    r2   = r2_score(y_te, y_pred)
    mae  = np.mean(np.abs(y_te - y_pred))
    rmse = np.sqrt(np.mean((y_te - y_pred)**2))
    flag = '← always predict mean' if model is None else            '← worse than baseline!' if r2 < 0 else ''
    print(f"  {name:<33}  {r2:>8.4f}  {mae:>8.4f}  {rmse:>8.4f}  {flag}")

# ── Adjusted R² — penalises adding useless features ───────────────────
# Standard R² always increases when you add more features (even noise)
# Adjusted R² penalises extra features that do not improve the model
def adjusted_r2(r2, n, p):
    """r2: R², n: samples, p: number of features"""
    return 1 - (1 - r2) * (n - 1) / (n - p - 1)

for n_noise in [0, 5, 10, 20]:
    X_noise = np.hstack([X_tr, np.random.randn(len(X_tr), n_noise)])
    X_te_n  = np.hstack([X_te, np.random.randn(len(X_te), n_noise)])
    lr = LinearRegression().fit(X_noise, y_tr)
    y_pred_n = lr.predict(X_te_n)
    r2_val  = r2_score(y_te, y_pred_n)
    adj_r2  = adjusted_r2(r2_val, len(X_te_n), X_noise.shape[1])
    n_feat  = 3 + n_noise
    print(f"  {n_feat} features ({n_noise} noise): R²={r2_val:.4f}  Adj R²={adj_r2:.4f}"
          f"  {'← noise inflated R²' if n_noise > 0 else ''}")

import numpy as np
from sklearn.metrics import (mean_absolute_error, mean_squared_error,
                              mean_absolute_percentage_error, r2_score)
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_validate, KFold
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 3000
distance = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic  = np.random.randint(1, 11, n).astype(float)
prep     = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
delivery = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
            + np.random.normal(0, 4, n)).clip(10, 120)
X = np.column_stack([distance, traffic, prep])
y = delivery

pipeline = Pipeline([
    ('sc', StandardScaler()),
    ('m',  GradientBoostingRegressor(n_estimators=200, learning_rate=0.1,
                                      max_depth=3, random_state=42)),
])

# ── Always compare against a baseline ─────────────────────────────────
from sklearn.model_selection import train_test_split
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)
pipeline.fit(X_tr, y_tr)
y_pred = pipeline.predict(X_te)

# Baselines
y_mean   = np.full_like(y_te, y_tr.mean())
y_median = np.full_like(y_te, np.median(y_tr))

print("Model vs baselines — DoorDash delivery time:")
print(f"{'Model':<25} {'MAE':>8} {'RMSE':>8} {'MAPE%':>8} {'R²':>8}")
print("─" * 62)

for label, pred in [
    ('Always predict mean',   y_mean),
    ('Always predict median', y_median),
    ('Gradient Boosting',     y_pred),
]:
    mae  = mean_absolute_error(y_te, pred)
    rmse = np.sqrt(mean_squared_error(y_te, pred))
    mape = mean_absolute_percentage_error(y_te, pred) * 100
    r2   = r2_score(y_te, pred)
    print(f"  {label:<23}  {mae:>8.4f}  {rmse:>8.4f}  {mape:>8.2f}  {r2:>8.4f}")

# ── Cross-validate multiple metrics at once ────────────────────────────
def rmse_scorer(est, X, y):
    return -np.sqrt(mean_squared_error(y, est.predict(X)))

from sklearn.metrics import make_scorer
cv_results = cross_validate(
    pipeline, X, y,
    cv=KFold(5, shuffle=True, random_state=42),
    scoring={
        'mae':  make_scorer(mean_absolute_error, greater_is_better=False),
        'rmse': make_scorer(rmse_scorer, greater_is_better=False),
        'r2':   'r2',
        'mape': make_scorer(mean_absolute_percentage_error,
                             greater_is_better=False),
    },
    return_train_score=False,
)

print(f"
5-fold CV across all metrics:")
print(f"  MAE:   {-cv_results['test_mae'].mean():.4f} ± {cv_results['test_mae'].std():.4f} min")
print(f"  RMSE:  {-cv_results['test_rmse'].mean():.4f} ± {cv_results['test_rmse'].std():.4f} min")
print(f"  MAPE:  {-cv_results['test_mape'].mean()*100:.2f}% ± {cv_results['test_mape'].std()*100:.2f}%")
print(f"  R²:    {cv_results['test_r2'].mean():.4f} ± {cv_results['test_r2'].std():.4f}")

# ── RMSE vs MAE ratio reveals outlier presence ─────────────────────────
rmse_val = -cv_results['test_rmse'].mean()
mae_val  = -cv_results['test_mae'].mean()
ratio    = rmse_val / mae_val
print(f"
RMSE / MAE ratio: {ratio:.2f}")
print("  Ratio ≈ 1.0: errors are uniform, no major outliers")
print("  Ratio > 2.0: large outlier errors present — investigate them")

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 3000
distance = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic  = np.random.randint(1, 11, n).astype(float)
prep     = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
delivery = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
            + np.random.normal(0, 4, n)).clip(10, 120)

X = np.column_stack([distance, traffic, prep])
y = delivery
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
pipeline = Pipeline([
    ('sc', StandardScaler()),
    ('m',  GradientBoostingRegressor(n_estimators=200, random_state=42)),
])
pipeline.fit(X_tr, y_tr)
y_pred = pipeline.predict(X_te)

residuals = y_te - y_pred   # positive = underestimate, negative = overestimate

print(f"Overall MAE: {mean_absolute_error(y_te, y_pred):.4f} min")

# ── Are residuals biased? ──────────────────────────────────────────────
mean_residual = residuals.mean()
print(f"
Mean residual: {mean_residual:+.4f} min")
print(f"  {'← model systematically underestimates' if mean_residual > 0.5 else '← model systematically overestimates' if mean_residual < -0.5 else '← no systematic bias'}")

# ── Error by delivery time bucket ─────────────────────────────────────
print(f"
MAE by delivery time bucket:")
for low, high in [(10,20),(20,30),(30,45),(45,60),(60,120)]:
    mask = (y_te >= low) & (y_te < high)
    if mask.sum() > 10:
        mae_bucket  = mean_absolute_error(y_te[mask], y_pred[mask])
        bias_bucket = residuals[mask].mean()
        print(f"  {low:2d}–{high:3d} min  n={mask.sum():4d}  "
              f"MAE={mae_bucket:.2f}  bias={bias_bucket:+.2f}")

# ── Error by distance bucket ──────────────────────────────────────────
dist_te = X_te[:, 0]   # distance feature in test set
print(f"
MAE by delivery distance:")
for low, high in [(0,2),(2,4),(4,7),(7,10),(10,15)]:
    mask = (dist_te >= low) & (dist_te < high)
    if mask.sum() > 10:
        mae_d = mean_absolute_error(y_te[mask], y_pred[mask])
        print(f"  {low:.0f}–{high:.0f} km   n={mask.sum():4d}  MAE={mae_d:.2f} min")

# ── Largest errors — what went wrong? ─────────────────────────────────
worst_idx   = np.argsort(np.abs(residuals))[-5:][::-1]
print(f"
Top 5 worst predictions:")
print(f"  {'Actual':>8} {'Predicted':>10} {'Error':>8} {'Dist':>6} {'Traffic':>8} {'Prep':>6}")
print("  " + "─" * 52)
for idx in worst_idx:
    print(f"  {y_te[idx]:>8.1f} {y_pred[idx]:>10.1f} {residuals[idx]:>+8.1f} "
          f"{X_te[idx,0]:>6.1f}  {X_te[idx,1]:>8.0f}  {X_te[idx,2]:>6.0f}")

Once a team settles on what "good" means, that definition frequently gets built directly into training, not just left for evaluation afterward. A delivery-time model whose real business metric is a percentile hit-rate is sometimes trained with quantile loss aimed directly at the ninetieth percentile, rather than the mean-squared-error loss that comes with sklearn's default regressor — training the model to be precisely accurate at the percentile that actually gates the SLA, instead of hoping that minimising average error happens to also fix the tail. A pricing model that cares about large errors on high-value items but does not want a handful of outliers to dominate training entirely might use a Huber loss, which behaves like squared error for small residuals and like absolute error for large ones — a training-time compromise that mirrors the same MAE-versus-RMSE tradeoff this module covers for evaluation.

import numpy as np

# ── Delivery-time model: an SLA hit-rate, not just MAE ─────────────────
# The business promise is "arrives within the estimate," which is a
# threshold-based hit-rate, not an average.
def delivery_sla_report(y_true, y_pred, sla_minutes=10):
    abs_error   = np.abs(y_true - y_pred)
    within_sla  = (abs_error <= sla_minutes).mean()
    p90_error   = np.percentile(abs_error, 90)
    return {
        'mae': abs_error.mean(),
        'within_sla_pct': within_sla * 100,
        'p90_error_min': p90_error,
    }

# ── Price model: revenue-weighted error, not plain MAPE ────────────────
# A cheap item's percentage error should not carry the same weight
# as an expensive item's percentage error in the business's eyes.
def price_wape_report(y_true, y_pred):
    abs_error  = np.abs(y_true - y_pred)
    wape       = abs_error.sum() / y_true.sum() * 100    # weighted by actual value
    plain_mape = (abs_error / y_true).mean() * 100        # unweighted, for comparison
    return {'wape_pct': wape, 'plain_mape_pct': plain_mape}

np.random.seed(42)
delivery_actual = np.random.normal(32, 8, 5000).clip(10, 90)
delivery_pred   = delivery_actual + np.random.normal(0, 4, 5000)
print("Delivery-time SLA report:")
for k, v in delivery_sla_report(delivery_actual, delivery_pred).items():
    print(f"  {k:<18}: {v:.2f}")

price_actual = np.concatenate([
    np.random.uniform(10, 200, 4000),        # everyday items
    np.random.uniform(5000, 100_000, 200),    # high-value items
])
price_pred = price_actual * np.random.normal(1.0, 0.06, len(price_actual))
print("\nPrice model report:")
for k, v in price_wape_report(price_actual, price_pred).items():
    print(f"  {k:<18}: {v:.2f}")
print("\nNotice: WAPE and plain MAPE diverge whenever error rates differ")
print("between cheap, high-volume items and rare, high-value ones.")

RMSE's quadratic penalty on large errors is a modelling choice, not a universal improvement — it only makes sense when large errors genuinely cost more than proportionally, like DoorDash's 45-minute delay triggering a refund. If your business cost is truly linear in the size of the error — being off by 10 minutes is exactly twice as bad as being off by 5, no more — MAE matches that cost structure honestly, and it is also far more robust to a handful of outliers dominating the reported number. Defaulting to RMSE because it's the convention, without checking whether your error costs are actually quadratic, means optimising and reporting against the wrong objective.

R² is relative to the variance of the target, not to any absolute error tolerance. A model with R²=0.92 on delivery times ranging from 10 to 120 minutes is explaining 92% of a large variance — but the remaining 8% unexplained can still translate into an MAE of 8 minutes, which may be operationally unacceptable even though the R² number looks excellent. This module's own guidance is explicit about this: report MAE or RMSE in the target's real units alongside R², because R² alone tells you nothing about whether the typical error is 30 seconds or 30 minutes.

Being in the same units does not mean being on the same scale. An RMSE of 5 minutes for delivery-time predictions is not comparable to an RMSE of 5 dollars for price predictions, and even within the same problem, an RMSE of 5 minutes on a dataset averaging 30-minute deliveries is a much larger relative error than an RMSE of 5 minutes on a dataset averaging 90-minute deliveries. RMSE and MAE are scale-dependent by construction — comparing them meaningfully across datasets or targets requires normalising first, whether that's dividing by the mean, reporting MAPE instead, or using a normalised RMSE (RMSE / range or RMSE / mean).

An aggregate MAE or RMSE is an average across every prediction, and averages hide systematic bias by construction. A model can post an excellent overall MAE of 4.2 minutes while being consistently 15 minutes late specifically for long-distance orders, or systematically biased for one customer segment — the aggregate number simply blends the good predictions with the bad ones. This is exactly why this module's residual analysis section exists: checking the mean residual, MAE by distance bucket, and MAE by delivery-time bucket separately is the only way to catch a model that looks fine in aggregate but is quietly failing a subgroup that never shows up in the headline metric.

MAPE divides by the actual value, so it breaks down — sometimes to infinity — whenever true values are zero or very close to zero, which is common in demand forecasting for low-volume products. It is also asymmetric in a way that's easy to miss: a prediction of twice the actual value produces a 100% error, but a prediction of half the actual value is capped at a 50% error, so MAPE structurally punishes overestimates more harshly than underestimates of the same relative size. It is genuinely useful for comparing errors across targets of very different scales, but "scale-independent" is not the same as "safe to use everywhere" — it should be avoided whenever the target can be zero or near zero.

The decision should follow the actual cost structure of your errors, not convention. Choose MAE when every unit of error costs proportionally the same — being off by 10 is exactly twice as bad as being off by 5, nothing more — and when you want a metric that isn't dominated by a handful of outliers. Choose RMSE when large errors are disproportionately costly in the real world, because squaring the error before averaging means a 10-minute miss contributes 4× more than a 5-minute miss, not 2×. A good answer also mentions checking the RMSE/MAE ratio in practice: a ratio near 1.0 means errors are fairly uniform and either metric tells a similar story; a ratio above 2.0 signals a few large outliers are inflating RMSE and worth investigating directly.

It depends on three things I'd check before answering. First, compared to what baseline — R²=0.91 sounds strong, but if a naive model (always predict the mean) already gets R²=0.85 on this target because the target itself is easy to predict, the real lift from the model is much smaller than 0.91 suggests. Second, what does that translate to in absolute error — R²=0.91 on a target with huge variance can still leave a large MAE in real units, which matters more operationally than the R² number itself. Third, is it measured on a held-out set with the same distribution as production — R² computed on training data or on a leaked split is not trustworthy at all. I wouldn't call any single R² value "good" without that context.

Negative R² means the model's squared errors are larger than they would be if you'd just predicted the mean of the target for every example — the model is actively worse than the simplest possible baseline. That's a strong signal something is broken, not just underperforming. I'd check, in order: whether train and test come from the same distribution (compare y_train.mean() and y_test.mean() — a big gap points to a bad split); whether the target was transformed during training (e.g. predicting log(y)) but evaluated without un-transforming the predictions back to y's scale; and whether the features used at inference time actually match what the model was trained on. Negative R² is rarely a subtle modelling issue — it's almost always a pipeline bug.

MAPE divides each error by the actual value, so for SKUs with true demand near zero — a product that sells 1 or 2 units a day — a small absolute error like being off by 3 units produces a triple-digit or even undefined percentage error. Those low-volume SKUs then dominate the averaged MAPE even though their absolute business impact is tiny, while high-volume SKUs where the forecast actually matters most get comparatively little weight in the metric. A better choice here is often a weighted MAE (weighted by revenue or volume) or WAPE (weighted absolute percentage error, which divides the sum of absolute errors by the sum of actuals rather than averaging per-item ratios), because both avoid the near-zero-denominator blowup that plain MAPE is vulnerable to.

I'd start with residual analysis rather than trusting the aggregate number. First check whether the mean residual is near zero overall — if it's shifted, the model has a global bias, not just a segment-specific one. Then break MAE down by the segment in question (and a few related cuts — by prediction range, by a key input feature) to see if that segment's error is meaningfully higher than the rest, and check the sign of its mean residual to see whether the model over- or under-predicts for that group specifically. An aggregate metric is a weighted average across every prediction, so a model can look excellent overall while being consistently wrong for a segment that's simply outnumbered by the rest of the data — the fix is always to disaggregate before concluding the model is fine.