How to Predict Numbers with Machine Learning in Python: Regression for Beginners (Food Delivery ML Series, Part 2)

The short answer

Regression is machine learning that predicts a number instead of a category. In this guide you will predict how many minutes a food delivery will take, getting within 4 minutes on average, which turns out to be as good as it is possible to get with this data.

Quick answer: To build a regression model in Python: (1) load your data with pandas, (2) choose the number you want to predict as the target, (3) split into training and test sets, (4) train LinearRegression() with .fit(), (5) measure errors with MAE, RMSE and R² against a simple average baseline, (6) improve with better features or tree-based models, and (7) predict new values and save the model with joblib.

What you will build

In Part 1 you built a model that answers a yes/no question: will this order be late? That is useful, but a customer staring at the app wants something more specific: "When will my food arrive?"

In this part you take the same 2,000 orders and predict the exact delivery time in minutes. Along the way you will:

• learn the difference between classification (categories) and regression (numbers);
• understand exactly what a model learns: features, weights, the bias, parameters and hyperparameters;
• measure a model's error in minutes, a unit everyone understands, and learn what MAE, RMSE, R² and the loss function really mean;
• find a model's blind spots by looking at the pattern of its errors;
• discover why a model can't get past a certain accuracy, however clever it is;
• teach a simple model about rain and distance working together (called interaction features), and watch it rediscover our hidden recipe almost exactly;
• turn predicted minutes back into an on-time/late warning and compare it with Part 1.

How this part connects to the series

This is Part 2 of the Delivery ML series. It reuses Part 1's dataset code line for line, with one change: instead of throwing away the delivery time and keeping only "late or not", we keep the minutes. If you completed Part 1, about half of this guide will feel familiar. That is deliberate: you learn only one new idea at a time.

Prerequisites

You need basic Python and the ideas from Part 1: features, labels, train/test split and overfitting. No new software, no maths beyond averages.

From earlier in the series

Knowledge

Software

The same setup as Part 1. If your environment from Part 1 still works, skip Step 1.

Hardware, time and cost

• Computer: any laptop, or a browser with Google Colab. Every model here trains in a few seconds.
• Time: 45–60 minutes reading every explanation; about 15 minutes if you only run the code.
• Cost: free. No downloads, GPU or paid API.

Theory break: classification vs regression

Both are supervised learning: the model learns from examples with known answers. The only difference is the kind of answer. Classification picks a category; regression predicts a number.

The quickest test: can you sensibly average two answers? The average of 40 and 50 minutes is 45 minutes, which makes sense, so that's regression. The average of "late" and "on time" means nothing, so that's classification.

Why predict the number at all?

A number carries more information than a category:

• Different decisions need different cut-offs. The app promises 45 minutes, but a premium service might promise 30. One minutes model serves both; a late/on-time model must be retrained for each promise.
• "How late" matters. 46 minutes and 80 minutes are both "late", but only one deserves an apology voucher.
• You can always turn a number into a category (above 45 means late), but you can't turn "late" back into minutes. You'll do exactly this in Step 7.

A beginner's trap: "logistic regression"

Part 1 used logistic regression, and this part uses linear regression. Despite the shared word, logistic regression is a classification model (it outputs a probability of a category). Linear regression is the real regression model. Blame history, not yourself.

Step 1: Set up your environment

Reuse your Part 1 environment, or open a fresh Google Colab notebook. Nothing new needs installing.

Coming from Part 1 on your own computer? Open a terminal in your project folder and reactivate the virtual environment:

cd first-ml-model

# Windows
venv\Scripts\activate

# Mac / Linux
source venv/bin/activate

Starting fresh? Follow Step 1 of Part 1, or install everything in one line inside an active virtual environment:

pip install scikit-learn pandas numpy joblib jupyter

Using Colab? Create a new notebook; everything is preinstalled.

Check your scikit-learn version

This guide uses root_mean_squared_error, added in scikit-learn 1.4:

import sklearn
print(sklearn.__version__)

If you see a version below 1.4, upgrade with pip install -U scikit-learn.

Tip: keep this tutorial in a new notebook or file, such as predict_delivery_time.py, next to your Part 1 file. Later parts of the series build on both.

Step 2: Recreate the orders and explore delivery times

We rebuild Part 1's 2,000 orders but keep the actual delivery time. The average order takes 40 minutes, and rain adds about 10 minutes on average.

The same hidden recipe, one change

As a reminder, each order's delivery time is kitchen prep time, plus travel time (about 3 minutes per km, 50% slower in rain, 30% slower at peak hour), minus 0.1 minutes per month of rider experience, plus random luck of typically ±5 minutes.

In Part 1 we turned that time into a yes/no label (is_late). Now the time itself is our target:

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)   # same seed as Part 1: the same 2,000 orders
n = 2000

orders = pd.DataFrame({
    "distance_km": rng.uniform(0.5, 12, n).round(1),
    "prep_time_min": rng.integers(5, 36, n),
    "is_raining": (rng.random(n) < 0.2).astype(int),
    "is_peak_hour": (rng.random(n) < 0.35).astype(int),
    "rider_experience_months": rng.integers(0, 49, n),
})

travel = orders.distance_km * 3 * (1 + 0.5 * orders.is_raining) * (1 + 0.3 * orders.is_peak_hour)
total = orders.prep_time_min + travel - 0.1 * orders.rider_experience_months + rng.normal(0, 5, n)
orders["delivery_min"] = total.clip(lower=5).round(1)   # the new target, in minutes

What does **.clip(lower=5)** do? Random luck occasionally makes a very short order come out at zero or even negative minutes, which is impossible. clip sets a floor of 5 minutes. It affects only 17 of the 2,000 orders. Checking that your target values are physically possible is a habit worth building on every regression project.

The dataset code, line by line

If any line above looked mysterious, here's what each one does. You'll reuse this pattern in every part of the series.

The 0/1 trick in the travel line is worth noticing: multiplying by a yes/no column switches an effect on or off. You'll use the same idea to build interaction features in Step 6.

Look at the first orders

orders.head()

Output:

   distance_km  prep_time_min  is_raining  is_peak_hour  rider_experience_months  delivery_min
0          9.4             14           0             1                       19          42.5
1          5.5             31           0             0                       36          45.9
2         10.4              6           0             0                        4          37.4
3          8.5             18           0             1                       44          48.6
4          1.6             29           1             1                       25          38.2

These are the same orders as Part 1. Order 1 took 45.9 minutes, which is why Part 1 labelled it late. Order 0 took 42.5, so it was on time. Now the model sees how late or early each order was, not just which side of 45 it fell on.

Summarise the target

print(orders["delivery_min"].describe().round(1))
count    2000.0
mean       40.4
std        16.5
min         5.0
25%        28.5
50%        39.6
75%        51.4
max        97.2

Remember the 16.5-minute spread. A model that just guesses the average will be off by roughly that much. Our job is to shrink it.

Which facts move delivery time?

print(orders.groupby("is_raining")["delivery_min"].mean().round(1))
print(orders.groupby("is_peak_hour")["delivery_min"].mean().round(1))

Rain adds about 10 minutes and rush hour about 5.5, on average. Keep those numbers in mind; the model will find something very similar in Step 4.

Step 3: Split the data and set a baseline

Hold back 400 orders for testing, then measure the laziest possible "model": always guess the average. It is off by 12.5 minutes per order. Every real model must beat that.

Split into training and test sets

from sklearn.model_selection import train_test_split

X = orders.drop(columns="delivery_min")   # the 5 features
y = orders["delivery_min"]                  # the target, in minutes

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

This is the same as Part 1 with one difference: no **stratify=y**. Stratifying keeps the share of each category equal in both sets. A continuous number like 42.5 minutes has no categories to balance, so a plain random split is the norm for regression.

The baseline: always guess the average

scikit-learn has a ready-made "dummy" model that ignores the features entirely and predicts the training average for every order:

from sklearn.dummy import DummyRegressor
from sklearn.metrics import mean_absolute_error

baseline = DummyRegressor(strategy="mean").fit(X_train, y_train)
print(f"Baseline MAE: {mean_absolute_error(y_test, baseline.predict(X_test)):.2f} minutes")

Output:

Baseline MAE: 12.54 minutes

The baseline predicts 40.3 minutes (the training average) for every single order, and is wrong by 12.54 minutes on average. MAE (mean absolute error) is explained properly in Step 5; for now read it as "average minutes off".

In Part 1, the lazy baseline was "always say on time" (62% accuracy). The regression version is "always say the average". Same lesson: a model is only useful if it clearly beats the lazy guess.

Step 4: Train a linear regression and read its formula

Linear regression learns a simple recipe: a starting number plus a fixed number of minutes for each feature. Ours learns that each kilometre adds 3.6 minutes and rain adds 10.6, close to the truth.

Train it

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(X_train, y_train)

Same .fit() as Part 1; only the model class changed.

How linear regression thinks

Linear regression assumes the answer is a weighted sum:

\text{minutes} = b + w_1 \times \text{distance} + w_2 \times \text{prep time} + w_3 \times \text{rain} + w_4 \times \text{peak} + w_5 \times \text{experience}

• b is the intercept, the starting number.
• w₁ … w₅ are the coefficients, the minutes each feature adds per unit.

Training finds the values of b and w that make the predictions as close as possible to the real delivery times in the training data. Technically it minimises the sum of squared errors, which is why it's also called ordinary least squares. You don't need the maths; scikit-learn solves it instantly.

Read what it learned

print(f"Intercept: {model.intercept_:.2f}")
for feature, coef in zip(X.columns, model.coef_):
    print(f"{feature:25s} {coef:+.2f}")

Output:

Intercept: -4.74
distance_km               +3.63
prep_time_min             +1.02
is_raining                +10.56
is_peak_hour              +5.67
rider_experience_months   -0.09

In Part 1, logistic regression's weights were abstract "pushes" towards late. Here every coefficient is in minutes, so you can read it like a recipe:

Do a prediction by hand

For a 6 km order, 18 minutes of prep, raining, at peak hour, with a rider of 6 months:

−4.74 + (3.63 × 6) + (1.02 × 18) + 10.56 + 5.67 + (−0.09 × 6)
= −4.74 + 21.78 + 18.36 + 10.56 + 5.67 − 0.54
= 51.1 minutes

model.predict() does exactly this sum for every order. That transparency is why linear regression is still used everywhere, from pricing to insurance.

Spot the flaw

Look at rain again. The model says rain adds a fixed 10.6 minutes, whether the order is 1 km or 12 km away. But in reality, rain makes travel 50% slower: that's about 1.5 extra minutes on a 1 km trip and 18 extra minutes on a 12 km trip. A plain linear model can't express "rain hurts more on long trips". Remember this; we fix it in Step 6.

Theory break: features, weights, parameters and hyperparameters

Features are the facts you give the model. Weights and the bias are what the model learns; together they're its parameters. Hyperparameters are settings you choose before training. The quickest test: who decides the value, and when?

You've now seen all four in action, so this is the right moment to name them clearly.

In the formula from earlier, every piece has a name:

minutes = −4.74 + 3.63 × distance + 1.02 × prep time + 10.56 × rain + 5.67 × peak − 0.09 × experience
          ↑ bias   ↑ weight ↑ feature

The weights and bias stay the same for every order; only the features change from order to order.

The bias, in a little more detail

• Why it's needed: without it, the formula would be forced through zero, meaning a 0 km order with 0 minutes of prep must take exactly 0 minutes. The bias lets the model shift every prediction up or down to fit reality.
• How it's found: LinearRegression calculates all six parameters at once with a formula, in a single step. Many other models (logistic regression, neural networks) learn by repeating small adjustments, and their bias changes a little on every step until it settles.
• It's frozen after training. Every prediction uses −4.74. Retrain on new data and you'll get a slightly different value.
• You can switch it off with LinearRegression(fit_intercept=False). Notice the difference: fit_intercept is a hyperparameter you choose, while the intercept it controls is a parameter the model learns.

Watch out: "bias" has two other meanings in machine learning. High bias means a model is too simple to fit the pattern (underfitting; you're about to see an example in Step 5). Bias in the news usually means unfairness towards a group of people. Context tells you which is meant.

The rule of thumb in code

model = RandomForestRegressor(n_estimators=200)   # inside the brackets: hyperparameters (you choose)
model.fit(X_train, y_train)                         # X_train: features (from the data)
model.coef_, model.intercept_                       # ending in _ after .fit(): learned parameters

(The last line works for LinearRegression; a random forest's learned parameters are its tree rules instead.)

Analogy: a student preparing for an exam. Each exam question's details are the features. The knowledge the student builds up by practising is the parameters. The study plan chosen beforehand (how many hours, how many practice papers) is the hyperparameters. Part 5 of this series is all about choosing good hyperparameters.

Step 5: Evaluate with MAE, RMSE and R²

The model is off by 4.6 minutes on an average order, against 12.5 for the baseline, and explains 86% of the variation in delivery times. 92% of its predictions land within 10 minutes.

In regression there's no "right" or "wrong", only how far off. Every prediction has an error (also called a residual):

\text{error} = \text{actual minutes} - \text{predicted minutes}

A negative error means the model was too pessimistic (food arrived earlier than predicted); a positive error means the food arrived later than predicted.

See the errors order by order

y_pred = model.predict(X_test)

results = pd.DataFrame({"actual": y_test.values, "predicted": y_pred.round(1)})
results["error"] = (results["actual"] - results["predicted"]).round(1)
print(results.head())

Output:

   actual  predicted  error
0    40.4       44.1   -3.7
1    36.8       37.9   -1.1
2    53.8       66.7  -12.9
3    43.1       47.0   -3.9
4    24.8       32.0   -7.2

Most errors are a few minutes. Order 2 is a big miss: predicted 66.7, actually 53.8. Scoring metrics squeeze all 400 of these errors into a single number.

The three metrics you need

from sklearn.metrics import mean_absolute_error, root_mean_squared_error, r2_score

print(f"MAE:  {mean_absolute_error(y_test, y_pred):.2f} minutes")
print(f"RMSE: {root_mean_squared_error(y_test, y_pred):.2f} minutes")
print(f"R2:   {r2_score(y_test, y_pred):.3f}")

Output:

MAE:  4.56 minutes
RMSE: 5.81 minutes
R2:   0.858

On just the five errors shown above, MAE is 5.8 minutes but RMSE is 7.0: squaring lets the single 12.9-minute miss dominate. Across all 400 test orders, the same effect makes RMSE (5.81) larger than MAE (4.56).

Which one should you report? Use MAE for business people ("we're usually within about 5 minutes"). Watch RMSE when big misses are especially bad (a 30-minute miss is much worse than three 10-minute misses). Use R² to compare models on the same data, but never alone, because it hides the size of the errors in minutes.

A friendlier way to say it

within_5 = (results["error"].abs() <= 5).mean()
within_10 = (results["error"].abs() <= 10).mean()
print(f"Within 5 minutes: {within_5:.0%}, within 10 minutes: {within_10:.0%}")
Within 5 minutes: 64%, within 10 minutes: 92%

For a product manager, "92% of our estimates are within 10 minutes" is more meaningful than any formula. Translating metrics into statements like this is a real professional skill.

Theory break: how to judge a regression model

Every regression metric is built from the same thing: the errors, one per prediction. MAE averages their size, RMSE punishes big ones, R² compares the model with simply guessing the average, and the pattern of errors shows where the model is blind.

Step 5 gave you three numbers. This section shows where each comes from, what it's good for, and the one check most beginners skip. Everything uses this model's real test results.

1. It all starts with the residuals

For every test order:

\text{error (residual)} = \text{actual minutes} - \text{predicted minutes}

We have 400 residuals, one per test order. Every metric below squeezes those 400 numbers into one.

2. Mean error: is the model biased in one direction?

Average the residuals with their signs:

errors = y_test.values - y_pred
print(round(errors.mean(), 2))   # 0.15

Our mean error is 0.15 minutes: close to zero, so overall the model isn't consistently too optimistic or too pessimistic. But this is a weak check on its own, because big positive and negative errors cancel out. A model that's 10 minutes late half the time and 10 minutes early the other half also scores 0. That's why we need the next metrics.

3. MAE: the typical size of a mistake

\text{MAE} = \text{average of } |\text{error}| = 4.56 \text{ minutes}

Drop the minus signs, then average. Every minute of error counts the same, whether it's part of a 1-minute miss or a 20-minute miss.

• Best for: explaining to anyone. "Our estimates are usually within about 4.6 minutes."
• Same units as the target, so it's immediately meaningful.

4. MSE and RMSE: punishing the big misses

\text{MSE} = \text{average of } \text{error}^2 = 33.7 \qquad \text{RMSE} = \sqrt{\text{MSE}} = \sqrt{33.7} = 5.81 \text{ minutes}

Squaring makes big errors count much more: a 2-minute miss counts as 4, but a 20-minute miss counts as 400. RMSE then takes the square root to get back to minutes.

• MSE is what linear regression actually minimises during training. It's the model's loss function: the number .fit() tries to make as small as possible. Its units ("minutes squared") are awkward to report, so we usually report RMSE instead.
• RMSE is always at least as large as MAE. The bigger the gap between them, the more your error comes from a few large misses. Ours: 5.81 vs 4.56.
• Best for: problems where one big miss is much worse than several small ones. A customer kept waiting 25 extra minutes is far angrier than five customers each waiting 5.

The diagram above shows the same effect on just five errors.

5. R²: how much better than guessing the average?

R² compares your model's squared errors with those of the lazy baseline that always predicts the average:

R^2 = 1 - \frac{\text{model's total squared error}}{\text{baseline's total squared error}} = 1 - \frac{13{,}480}{95{,}201} = 0.858

Read it as: "the model removed 86% of the squared error that guessing the average would make." People often say it "explains 86% of the variation in delivery times".

Caution: R² has no units, so it can't tell you whether errors are acceptable. An R² of 0.9 could mean 4-minute errors on deliveries or 400-minute errors on cross-country freight. Always report MAE or RMSE in real units next to it. R² is most

useful for comparing models on the same data.

6. The check most beginners skip: look at the pattern of errors

A good overall score can hide a model that's systematically wrong for certain kinds of order. Group the residuals and average them with their signs:

results_check = X_test.copy()
results_check["error"] = y_test.values - y_pred

rainy = results_check[results_check["is_raining"] == 1]
print("Rain, long trips (> 8 km): ", round(rainy[rainy["distance_km"] > 8]["error"].mean(), 1))
print("Rain, short trips (<= 3 km):", round(rainy[rainy["distance_km"] <= 3]["error"].mean(), 1))
print("Dry weather:                ", round(results_check[results_check["is_raining"] == 0]["error"].mean(), 1))
Rain, long trips (> 8 km):  6.1
Rain, short trips (<= 3 km): -6.2
Dry weather:                 0.2

This is the flaw from Step 4 caught in the data. Adding a fixed 10.6 minutes for rain is too much for short trips and too little for long ones. The errors cancel out in the overall mean error (0.15), which is exactly why looking at groups matters. A model that's systematically wrong in a predictable direction like this has high bias (it's too simple for the pattern); compare it with the overfitting decision tree, which has high variance.

In Step 6, the interaction features fix it. Rerun this check with the improved model and the long-trip error falls from +6.1 to +0.9 minutes, and the short-trip error from −6.2 to −0.5. Whenever a model disappoints, group its errors by your features; the pattern usually tells you which feature to engineer next.

7. A metric to be careful with: percentage error (MAPE)

Businesses love "we're 13% off on average" (our model's mean absolute percentage error is 13.4%). But percentages explode for small values: a 3-minute miss on a 5-minute order is a 60% error, while the same 3 minutes on a 60-minute order is 5%. Use MAPE only when the target never gets close to zero, and report it alongside MAE.

8. Which metric should you use?

Quick self-check

1. Model A has MAE 5 and RMSE 5.5. Model B has MAE 5 and RMSE 9. What's different about their mistakes?
2. A model's mean error is exactly 0. Does that mean it's accurate?
3. Someone reports "R² = 0.95" and nothing else. What's the first thing you'd ask?

Answers: (1) Model B makes some much bigger misses: its RMSE is far above its MAE, so a few large errors dominate. Model A's errors are more even. (2) No. Large errors in both directions can cancel out; check MAE and the error pattern by group. (3) "What's the MAE or RMSE in real units?" R² alone can't tell you whether the errors are acceptable.

Step 6: Improve the model

Three extra columns that combine distance with rain and rush hour cut the error from 4.6 to 3.8 minutes, and the model rediscovers our hidden recipe almost to the decimal. No algorithm can do much better, because the remaining error is pure luck.

Fix the flaw with interaction features

In Step 4 the model couldn't express "rain hurts more on long trips". The fix is to give it a feature that represents exactly that: distance multiplied by rain. When it's dry, the new column is 0; when it rains, it equals the distance.

def add_interactions(df):
    df = df.copy()
    df["distance_x_rain"] = df["distance_km"] * df["is_raining"]
    df["distance_x_peak"] = df["distance_km"] * df["is_peak_hour"]
    df["distance_x_rain_x_peak"] = df["distance_km"] * df["is_raining"] * df["is_peak_hour"]
    return df

X_train_fe = add_interactions(X_train)
X_test_fe = add_interactions(X_test)

better = LinearRegression().fit(X_train_fe, y_train)
for feature, coef in zip(X_train_fe.columns, better.coef_):
    print(f"{feature:25s} {coef:+.2f}")

Output:

distance_km               +3.01
prep_time_min             +1.02
is_raining                -0.23
is_peak_hour              +0.37
rider_experience_months   -0.09
distance_x_rain           +1.62
distance_x_peak           +0.80
distance_x_rain_x_peak    +0.38

Now compare with the recipe we hid in Step 2. Travel takes 3 minutes per km, ×1.5 in rain, ×1.3 at peak. Multiplied out, that is 3 + 1.5 (rain) + 0.9 (peak) + 0.45 (both) minutes per km:

The flat rain and peak effects (−0.23, +0.37) have shrunk to almost nothing, because the model now knows rain and rush hour act through distance. This is the power of feature engineering: the same simple algorithm, given one well-chosen idea, describes reality far better.

Compare with other models

Tree-based models can discover interactions by themselves. Let's compare everything, reporting both training and test error so we can spot overfitting:

from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
from sklearn.tree import DecisionTreeRegressor

models = {
    "Baseline (average)": (DummyRegressor(), X_train, X_test),
    "Linear regression": (LinearRegression(), X_train, X_test),
    "Linear + interactions": (LinearRegression(), X_train_fe, X_test_fe),
    "Decision tree": (DecisionTreeRegressor(random_state=42), X_train, X_test),
    "Random forest": (RandomForestRegressor(n_estimators=200, random_state=42), X_train, X_test),
    "Gradient boosting": (GradientBoostingRegressor(random_state=42), X_train, X_test),
}

for name, (m, Xtr, Xte) in models.items():
    m.fit(Xtr, y_train)
    pred = m.predict(Xte)
    train_mae = mean_absolute_error(y_train, m.predict(Xtr))
    print(f"{name:23s} train MAE={train_mae:5.2f}  test MAE={mean_absolute_error(y_test, pred):5.2f}"
          f"  RMSE={root_mean_squared_error(y_test, pred):5.2f}  R2={r2_score(y_test, pred):.3f}")

Output:

Baseline (average)      train MAE=13.55  test MAE=12.54  RMSE=15.43  R2=-0.001
Linear regression       train MAE= 4.58  test MAE= 4.56  RMSE= 5.81  R2=0.858
Linear + interactions   train MAE= 4.11  test MAE= 3.82  RMSE= 4.84  R2=0.901
Decision tree           train MAE= 0.00  test MAE= 5.89  RMSE= 7.62  R2=0.756
Random forest           train MAE= 1.75  test MAE= 4.41  RMSE= 5.57  R2=0.869
Gradient boosting       train MAE= 3.76  test MAE= 4.05  RMSE= 5.10  R2=0.891

You've seen the single decision tree overfit before, in Part 1. It's the same lesson in minutes: a training error of zero is a warning sign, not an achievement.

Theory break: the floor you can't get through

Why not keep going until the error is zero? Because of the random luck in our recipe: a traffic light, a lift, a wrong turn. We added it with a typical size of ±5 minutes (rng.normal(0, 5, n)), and no feature describes it.

No model can predict pure luck. Even a model that knew the hidden recipe perfectly would still be off by about 4 minutes MAE (5 minutes RMSE). Statisticians call this the irreducible error. Our best model scores 3.82 MAE and 4.84 RMSE: it has reached the floor (slightly below it, by luck of which orders landed in the test set).

In real projects you never know where the floor is, but the lesson holds: when several very different models all stall at about the same error, you've probably hit the limits of your data. The way forward is then better data and new features, such as live traffic or the restaurant's current queue, not a fancier algorithm.

Step 7: Predict new orders, convert to late/on time, and save

The model predicts 52 minutes for a 6 km order in rain at dinner time, 40 without the rain, and 35 off-peak, each within half a minute of the hidden truth. Turned into a late/on-time warning, it is right 90.8% of the time.

Predict new orders

The same three "what if" orders as Part 1, but now we get minutes instead of a probability:

what_if = pd.DataFrame([
    {"distance_km": 6.0, "prep_time_min": 18, "is_raining": 1, "is_peak_hour": 1, "rider_experience_months": 6},
    {"distance_km": 6.0, "prep_time_min": 18, "is_raining": 0, "is_peak_hour": 1, "rider_experience_months": 6},
    {"distance_km": 6.0, "prep_time_min": 18, "is_raining": 0, "is_peak_hour": 0, "rider_experience_months": 6},
])
print(better.predict(add_interactions(what_if)).round(1))

Output:

[52.2 40.4 35.2]

Important: new orders must go through add_interactions() too. The model was trained on 8 columns, so it must always receive the same 8 columns. Forgetting this step is one of the most common bugs in real projects. In Part 3 you'll learn to bundle feature steps and the model into a single pipeline so this can't happen.

The two models tell a consistent story, but the minutes are far more useful: the app can now say "arriving in about 52 minutes" instead of just "likely late".

Turn minutes into a late/on-time warning

A number can always become a category. Anything over 45 minutes is late:

late_pred = better.predict(X_test_fe) > 45
late_true = y_test > 45
print(f"Late/on-time accuracy from minutes: {(late_pred == late_true).mean():.1%}")
Late/on-time accuracy from minutes: 90.8%

Part 1's dedicated classifier scored 92.0% on the same question. So which is better?

• For a single yes/no decision, a model trained directly on that question usually wins by a little.
• For flexibility, the minutes model wins: change the promise to 30 or 60 minutes, or show a countdown, without retraining anything.

Many real systems use both: one model for the time estimate shown to customers, and one tuned specifically to flag risky orders for the operations team.

Save the model

import joblib

joblib.dump(better, "delivery_time_model.joblib")

loaded = joblib.load("delivery_time_model.joblib")
print(loaded.predict(add_interactions(what_if)).round(1))
[52.2 40.4 35.2]

The saved file holds the model, but not the **add_interactions** function. Whoever loads the model needs that function too. Keep it in a shared file, such as features.py, alongside the model. Part 4 of this series, where the model becomes an API, depends on getting this right.

The complete script

Here is the whole of Part 2 in one file. Save it as **predict_delivery_time.py** and run **python predict_delivery_time.py**.

"""Delivery ML Series, Part 2: predict the exact delivery time (regression)."""
import joblib
import numpy as np
import pandas as pd
from sklearn.dummy import DummyRegressor
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, r2_score, root_mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeRegressor

# 1. Recreate the 2,000 orders from Part 1, keeping the delivery time in minutes
rng = np.random.default_rng(42)
n = 2000
orders = pd.DataFrame({
    "distance_km": rng.uniform(0.5, 12, n).round(1),
    "prep_time_min": rng.integers(5, 36, n),
    "is_raining": (rng.random(n) < 0.2).astype(int),
    "is_peak_hour": (rng.random(n) < 0.35).astype(int),
    "rider_experience_months": rng.integers(0, 49, n),
})
travel = orders.distance_km * 3 * (1 + 0.5 * orders.is_raining) * (1 + 0.3 * orders.is_peak_hour)
total = orders.prep_time_min + travel - 0.1 * orders.rider_experience_months + rng.normal(0, 5, n)
orders["delivery_min"] = total.clip(lower=5).round(1)

# 2. Explore
print(orders.head())
print(orders["delivery_min"].describe().round(1))

# 3. Split and set a baseline
X = orders.drop(columns="delivery_min")
y = orders["delivery_min"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

baseline = DummyRegressor(strategy="mean").fit(X_train, y_train)
print(f"Baseline MAE: {mean_absolute_error(y_test, baseline.predict(X_test)):.2f} minutes")

# 4. Train a linear regression
model = LinearRegression()
model.fit(X_train, y_train)
for feature, coef in zip(X.columns, model.coef_):
    print(f"{feature:25s} {coef:+.2f}")

# 5. Evaluate
y_pred = model.predict(X_test)
print(f"MAE:  {mean_absolute_error(y_test, y_pred):.2f} minutes")
print(f"RMSE: {root_mean_squared_error(y_test, y_pred):.2f} minutes")
print(f"R2:   {r2_score(y_test, y_pred):.3f}")
results = pd.DataFrame({"actual": y_test.values, "predicted": y_pred.round(1)})
results["error"] = (results["actual"] - results["predicted"]).round(1)
print(results.head())
within_5 = (results["error"].abs() <= 5).mean()
within_10 = (results["error"].abs() <= 10).mean()
print(f"Within 5 minutes: {within_5:.0%}, within 10 minutes: {within_10:.0%}")


# 6. Improve with interaction features, then compare models
def add_interactions(df):
    df = df.copy()
    df["distance_x_rain"] = df["distance_km"] * df["is_raining"]
    df["distance_x_peak"] = df["distance_km"] * df["is_peak_hour"]
    df["distance_x_rain_x_peak"] = df["distance_km"] * df["is_raining"] * df["is_peak_hour"]
    return df


X_train_fe = add_interactions(X_train)
X_test_fe = add_interactions(X_test)

models = {
    "Baseline (average)": (DummyRegressor(), X_train, X_test),
    "Linear regression": (LinearRegression(), X_train, X_test),
    "Linear + interactions": (LinearRegression(), X_train_fe, X_test_fe),
    "Decision tree": (DecisionTreeRegressor(random_state=42), X_train, X_test),
    "Random forest": (RandomForestRegressor(n_estimators=200, random_state=42), X_train, X_test),
    "Gradient boosting": (GradientBoostingRegressor(random_state=42), X_train, X_test),
}
for name, (m, Xtr, Xte) in models.items():
    m.fit(Xtr, y_train)
    pred = m.predict(Xte)
    train_mae = mean_absolute_error(y_train, m.predict(Xtr))
    print(f"{name:23s} train MAE={train_mae:5.2f}  test MAE={mean_absolute_error(y_test, pred):5.2f}"
          f"  RMSE={root_mean_squared_error(y_test, pred):5.2f}  R2={r2_score(y_test, pred):.3f}")

best = models["Linear + interactions"][0]
for feature, coef in zip(X_train_fe.columns, best.coef_):
    print(f"{feature:25s} {coef:+.2f}")

# 7. Predict new orders, convert to late/on time, save
what_if = pd.DataFrame([
    {"distance_km": 6.0, "prep_time_min": 18, "is_raining": 1, "is_peak_hour": 1, "rider_experience_months": 6},
    {"distance_km": 6.0, "prep_time_min": 18, "is_raining": 0, "is_peak_hour": 1, "rider_experience_months": 6},
    {"distance_km": 6.0, "prep_time_min": 18, "is_raining": 0, "is_peak_hour": 0, "rider_experience_months": 6},
])
print(best.predict(add_interactions(what_if)).round(1))

late_pred = best.predict(X_test_fe) > 45
late_true = y_test > 45
print(f"Late/on-time accuracy from minutes: {(late_pred == late_true).mean():.1%}")

joblib.dump(best, "delivery_time_model.joblib")
loaded = joblib.load("delivery_time_model.joblib")
print(loaded.predict(add_interactions(what_if)).round(1))

Expected output (key lines):

Baseline MAE: 12.54 minutes
MAE:  4.56 minutes
RMSE: 5.81 minutes
R2:   0.858
Within 5 minutes: 64%, within 10 minutes: 92%
Linear + interactions   train MAE= 4.11  test MAE= 3.82  RMSE= 4.84  R2=0.901
[52.2 40.4 35.2]
Late/on-time accuracy from minutes: 90.8%

Common mistakes and how to fix them

Regression habits worth building now: check the target's range is physically possible, always compare with an "average" baseline, report errors in real units, and look at individual errors, not just the summary score.

Practice project: predict a home's exact electricity use

Back to the homes from Part 1's practice project. This time, instead of "over budget or not", predict the exact number of units (kWh) each home will use. You'll hit the same flaw this guide found, and fix it the same way.

Get the data

Use the same make_homes() function from Part 1's practice project (copy it again if you're starting fresh). The target is now monthly_kwh itself:

homes = make_homes()
X = homes.drop(columns="monthly_kwh")
y = homes["monthly_kwh"]

If you added an over_budget column in Part 1, drop that too: it's calculated from the target.

Your tasks

• 1. Explore (Step 2). Use describe() on monthly_kwh. What's the average and the spread (std)? What does a baseline that always guesses the average need to beat?
• 2. Split and baseline (Step 3). Split 80/20 with random_state=42 (no stratify for regression). What MAE does DummyRegressor() score?
• 3. Train linear regression (Step 4). Report MAE and R². Print the coefficients. They're in kWh per unit, so read them like a recipe: how many units does one extra hour of AC a day add? What do solar panels save?
• 4. Spot the flaw (Step 4). Look at the hidden recipe in make_homes(). AC use is multiplied by home size: an hour of AC in a 200 m² home uses far more than in a 40 m² flat. Can plain linear regression express that?
• 5. Fix it with an interaction feature (Step 6). Add a column ac_x_size = ac_hours_per_day × home_size_sqm and retrain. How much does MAE drop?
• 6. Challenge: the solar discount (Step 6). Solar cuts the whole bill by 35%, not by a fixed number of units. Add interaction columns that multiply has_solar by each of home_size_sqm, residents, ac_x_size, appliance_age_years and works_from_home. What MAE do you get now?
• 7. Find the floor (theory break). The recipe's random luck has a typical size of 30 kWh (rng.normal(0, 30, n)). Roughly what's the best MAE any model could reach? How close did you get?
• 8. Compare with gradient boosting (Step 6). Train GradientBoostingRegressor(random_state=42) on the original features (no interactions). How does it compare with your engineered linear models?

Check your results

Think about it: gradient boosting matched one hand-made feature for free, but your domain knowledge ("solar is a percentage discount") beat it. When would you choose each approach in a real project?

Share your notebook and answer in the comments, and keep it: Part 3's practice project hands you a messy version of this data.

Frequently asked questions

What is regression in machine learning?

Regression is supervised learning that predicts a number, such as a delivery time, a price or a temperature. It learns from examples whose correct numbers are known, then estimates the number for new examples.

What is the difference between classification and regression?

Classification predicts a category (late or on time); regression predicts a number (47 minutes). A quick test: if averaging two answers makes sense, it's regression.

How do I do linear regression in Python?

With scikit-learn: create LinearRegression(), call .fit(X_train, y_train), then .predict(X_test). Read the learned weights from .coef_ and the starting value from .intercept_.

What is a good MAE or R²?

It depends on the problem. Always compare with a baseline that predicts the average. Here the baseline's MAE is 12.5 minutes and our best model's is 3.8. An R² of 0.9 means the model explains 90% of the variation.

What is the difference between MAE and RMSE?

Both measure average error in the target's units. MAE treats every minute of error equally; RMSE squares errors first, so big misses count more. RMSE is always at least as large as MAE.

Why can't my model's error go to zero?

Because some variation is random luck that no feature captures. This irreducible error sets a floor. When very different models stall at the same error, you need better data, not a better algorithm.

What are interaction features?

New columns made by multiplying features together, such as distance × rain. They let simple models capture effects that depend on two things at once, like rain slowing long trips more than short ones.

Glossary

What to learn next: the Delivery ML series

This guide is Part 2 of an 8-part series. Every part uses the same food delivery orders and starts where the previous part ended.

Next up: Part 3, Clean Messy, Real-World Delivery Data. Keep predict_delivery_time.py and your add_interactions function: Part 3 builds them into a pipeline.

Conclusion

You've built a regression model that predicts delivery time to within about 4 minutes, which is the limit of what this data allows, and you've seen a simple model rediscover a hidden recipe from examples alone.

The big ideas carry over to every regression problem, from house prices to insurance claim amounts: beat an average baseline, report errors in real units, look at individual mistakes, give simple models the features they need, and recognise when you've hit the limit of your data. Try changing the luck in the recipe from 5 to 10 minutes, rerun everything, and see how the floor moves.