Linear Regression

In a real job you would pull this from BigQuery or a Postgres database. Here we simulate it with numpy so you can run it instantly with no setup.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

np.random.seed(42)          # same random numbers every run — reproducibility

n_orders = 1000
distance_km = np.random.uniform(0.5, 8.0, n_orders)    # orders range 0.5–8 km
noise = np.random.normal(0, 4, n_orders)                # ±4 min real-world variability
delivery_time_min = 8.6 + 7.3 * distance_km + noise    # true relationship + noise

df = pd.DataFrame({
    'distance_km': distance_km,
    'delivery_time_min': delivery_time_min
})

print(df.describe())        # always run describe() first — know your data before touching it

This is non-negotiable. A scatter plot in 30 seconds tells you whether Linear Regression is the right tool, before you write another line of code.

plt.figure(figsize=(8, 5))
plt.scatter(df['distance_km'], df['delivery_time_min'],
            alpha=0.3, s=15, color='#378ADD')
plt.xlabel('Distance (km)')
plt.ylabel('Delivery Time (min)')
plt.title('DoorDash Orders — Distance vs Delivery Time')
plt.tight_layout()
plt.show()

# What to look for:
#   Clear upward trend  → Linear Regression is appropriate
#   No trend / random   → LR will learn nothing useful
#   Curve (parabola)    → Wrong tool — try polynomial features or tree model
#   Tight cluster       → Great signal, model will be accurate

This is the single most important step beginners skip. If you evaluate on the same data you trained on, you are asking a student to grade their own exam using the answer sheet they already memorised. The score is meaningless.

X = df[['distance_km']]      # double brackets → 2D array (n_samples, n_features)
y = df['delivery_time_min']  # single brackets → 1D array (n_samples,)

X_train, X_test, y_train, y_test = train_test_split(
    X, y,
    test_size=0.2,       # 20% for testing = 200 orders
    random_state=42      # same split every run
)

print(f"Training set: {len(X_train)} orders")   # 800
print(f"Test set:     {len(X_test)} orders")    # 200

# The test set is the answer sheet.
# Treat it like one: don't look at it, don't tune on it, don't touch it
# until you are done building and ready for a final honest evaluation.

Three lines. That is all sklearn needs. The complexity is hidden inside fit().

model = LinearRegression()      # create the model object — no learning has happened yet

model.fit(X_train, y_train)    # THIS is where learning happens:
                                # - computes optimal slope and intercept via OLS
                                # - for 800 rows, this takes milliseconds
                                # - result: model.coef_ and model.intercept_ are set

print(f"Slope (coef):     {model.coef_[0]:.2f}")    # expect ~7.3 (min per km)
print(f"Intercept:        {model.intercept_:.2f}")  # expect ~8.6 (baseline minutes)

# These numbers ARE the model.
# coef_ = [7.28]    ← distance adds 7.28 min per km
# intercept_ = 8.71 ← base time before rider moves

A number without context is useless. Always compare your model to the dumbest possible baseline: predict the mean for every order. If your model cannot beat that, it has learned nothing.

mae = mean_absolute_error(y_test, y_pred)
print(f"MAE:  {mae:.2f} min")      # ~3.1 min on average — a PM can understand this

mse = mean_squared_error(y_test, y_pred)
rmse = np.sqrt(mse)
print(f"RMSE: {rmse:.2f} min")     # penalises large errors more heavily than MAE

r2 = r2_score(y_test, y_pred)
print(f"R²:   {r2:.3f}")           # ~0.78 — distance explains 78% of delivery time variation

# --- Compare to baseline ---
baseline_pred = np.full(len(y_test), y_train.mean())   # always predict training mean
baseline_mae  = mean_absolute_error(y_test, baseline_pred)
print(f"\nBaseline MAE (always predict mean): {baseline_mae:.2f} min")
print(f"Model MAE:                           {mae:.2f} min")
print(f"Improvement over baseline:           {((baseline_mae - mae) / baseline_mae * 100):.1f}%")

Numbers tell you how wrong you are. Charts tell you where and why. The residual plot is the most important diagnostic chart in Linear Regression.

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Left: scatter + fitted line
ax1.scatter(X_train['distance_km'], y_train, alpha=0.2, s=12, color='#378ADD', label='Training data')
x_line = np.linspace(0.5, 8.0, 100).reshape(-1, 1)
ax1.plot(x_line, model.predict(x_line), color='#ff6b6b', linewidth=2, label='Fitted line')
ax1.set_xlabel('Distance (km)')
ax1.set_ylabel('Delivery Time (min)')
ax1.set_title('Fitted Line')
ax1.legend()

# Right: residual plot — actual minus predicted vs predicted
residuals = y_test - y_pred
ax2.scatter(y_pred, residuals, alpha=0.3, s=12, color='#7F77DD')
ax2.axhline(0, color='#ff6b6b', linewidth=1.5, linestyle='--')
ax2.set_xlabel('Predicted Time (min)')
ax2.set_ylabel('Residual (actual − predicted)')
ax2.set_title('Residuals')

# Interpreting the residual plot:
#   GOOD → random cloud centred on 0, no pattern
#   BAD  → fan shape (wider spread at higher predictions) → heteroscedasticity
#   BAD  → curve → the relationship is not actually linear

plt.tight_layout()
plt.show()

# Add more features to the dataset
df['traffic_score']    = np.random.uniform(1, 10, n_orders)    # 1=clear 10=gridlock
df['restaurant_prep']  = np.random.uniform(5, 25, n_orders)   # minutes to prepare

# Re-split with all features
X_multi = df[['distance_km', 'traffic_score', 'restaurant_prep']]
y = df['delivery_time_min']
X_tr, X_te, y_tr, y_te = train_test_split(X_multi, y, test_size=0.2, random_state=42)

# Training is identical — sklearn handles multiple features automatically
model_multi = LinearRegression()
model_multi.fit(X_tr, y_tr)

# Inspect coefficients — one per feature
for feature, coef in zip(X_multi.columns, model_multi.coef_):
    print(f"{feature:22s}: {coef:.3f}")
# distance_km           :  7.312  ← 7.3 min per km (as expected)
# traffic_score         :  0.891  ← each traffic point adds ~0.9 min
# restaurant_prep       :  0.998  ← 1 extra prep minute = 1 extra delivery minute

y_pred_multi = model_multi.predict(X_te)
print(f"\nSimple LR R²:   {r2_score(y_test, y_pred):.3f}")        # ~0.78
print(f"Multiple LR R²: {r2_score(y_te, y_pred_multi):.3f}")      # ~0.92

Two more assumptions get silently violated far more often than the ones on this page's assumptions table: homoscedasticity (residual variance stays roughly constant across the range of predictions) and normality of the residuals. Homoscedasticity is checked with the same residual-versus-predicted plot from Step 7 above — look specifically for a fan or cone shape, spread widening as predictions grow, which shows up constantly in real-world data where absolute error naturally scales with the size of the thing being predicted, like larger orders having proportionally larger errors. Normality of residuals is checked with a Q-Q plot comparing residual quantiles to a theoretical normal distribution, and it mainly affects whether confidence intervals and p-values on the coefficients can be trusted, not the accuracy of the point predictions themselves. Both drift in gradually and are easy to miss unless you deliberately plot for them.

R squared only measures how much of this particular sample's variance in the target the fitted line accounts for — it says nothing about whether the model's assumptions hold, whether it will generalise to new orders, or why the relationship exists. The Multiple Linear Regression example earlier on this page reaches an R squared around 0.92 by adding traffic_score and restaurant_prep; that number looks identical whether traffic genuinely slows deliveries or whether traffic simply happens to be measured at the same times as some other unmodelled cause of delay. A high R squared earned on training data specifically can also just mean the model memorised patterns that will not hold on the held-out test set, which is exactly why this page insists on splitting the data before ever looking at R squared at all.

Multicollinearity mainly damages your ability to interpret individual coefficients, not the accuracy of the combined prediction. If distance_km and an estimated_travel_time feature are highly correlated, the model can split credit between them almost arbitrarily — one fit might give most of the weight to distance and a little to travel time, another equally valid fit might reverse that split — while the combined contribution to any given prediction stays nearly the same either way. That instability is a real problem if the goal is explaining which feature drives delivery time, but if the only goal is minimising prediction error on data that resembles the training distribution, a collinear model can predict just as well as one with a redundant feature removed. The real danger shows up later: if the correlation between those features shifts in production away from what training data showed, a model that leaned on that correlation can degrade sharply.

Ordinary Least Squares fits parameters by minimising squared error on whatever training data it is given, and handing it one more column — even a genuinely useless one, like a random number generator's output — gives the optimiser strictly more freedom to reduce that training error further. R squared on the training set is mathematically guaranteed to only increase or stay flat as features are added; it can never go down. That is exactly why adjusted R squared exists: it subtracts a penalty that grows with the number of predictors relative to the sample size, so a feature that is not pulling its weight can push adjusted R squared down even while plain R squared ticks up. Comparing models with different feature counts using plain R squared instead of the adjusted version, or instead of test set performance, is one of the most common ways beginners convince themselves a bigger model is a better one.

A fitted coefficient captures the association between a feature and the target after accounting for whatever other features happen to be in the model, on this particular dataset — it is a correlational quantity, not a causal one, no matter how large or stable it looks. If a restaurant_popularity feature came out with a large positive coefficient on delivery time, the honest read is only that popular restaurants and long delivery times moved together in this data; a very plausible alternative story is that popular restaurants get overwhelmed with orders and understaff their kitchen during exactly those busy hours — a confounder, not an effect of popularity itself. Distance is the rare feature on this page where correlation and causation likely do line up, since physically travelling further mechanically takes more time, but that confidence comes from reasoning about how delivery works in the real world, not from anything the regression itself proved. Establishing causation in general needs a controlled experiment or a dedicated causal-inference method, never a coefficient by itself.

I would check each one with a specific diagnostic rather than assuming from a scatter plot. Linearity and homoscedasticity both show up in a residuals-versus-predicted-values plot: a random cloud centred on zero with roughly constant spread supports both, a curved pattern means the true relationship is not linear, and a fan or cone shape means variance is not constant. Normality of residuals is checked with a Q-Q plot comparing residual quantiles to a theoretical normal distribution, and it mainly affects whether confidence intervals and p-values on the coefficients can be trusted, not the point predictions themselves. Independence of errors is checked by plotting residuals against time, order, or any variable outside the model — a pattern there means a systematic effect is missing from the model. Multicollinearity is checked with a correlation matrix or, more precisely, Variance Inflation Factor. I would emphasise that violated assumptions in ordinary least squares mostly threaten how much you can trust the model's inference — its confidence intervals, its p-values, how interpretable individual coefficients are — while the point predictions themselves can often still be reasonably useful even under a moderate violation.

R squared alone does not answer that, and I would ask a few things before trusting it. First, is 0.95 measured on a held-out test set or on the training data — a strong score on training data alone says almost nothing about generalisation, and the gap between training and test R squared is often more informative than either number alone. Second, how does it compare to a naive baseline, such as always predicting the mean — a workflow like the DoorDash example on this page always computes that baseline MAE first specifically so a strong-looking metric can be judged against a floor. Third, a high R squared does not validate the model's assumptions or imply the relationships found are causal; it is entirely possible to reach a high R squared with a leaked feature that is only available because the outcome already happened, or with a relationship that will not hold once the underlying data distribution shifts. I would want to see the baseline comparison, the train-versus-test gap, and a residual plot before calling 0.95 good rather than lucky or leaky.

R squared can only increase or stay the same as you add features to an Ordinary Least Squares model, because the optimiser is minimising training error and an extra column — even a useless one — gives it strictly more freedom to reduce that error further; there is no mechanism by which adding a feature can make training R squared go down. Adjusted R squared fixes this by including a penalty term that grows with the number of predictors relative to the number of observations, so a feature that does not meaningfully reduce error can cause adjusted R squared to fall even while plain R squared rises slightly. In practice I use adjusted R squared, or better, held-out test performance, whenever I am comparing two models with different numbers of features, since plain R squared is always biased toward the larger model regardless of whether the extra features are genuinely useful.

The core issue is that the model can no longer cleanly separate how much each of the two correlated features individually contributes, because many different splits of coefficient weight between them produce almost the same combined prediction — that makes individual coefficients unstable and their standard errors unreliable, so a statement like "feature A matters twice as much as feature B" cannot be trusted when the two are highly collinear. Whether it matters in production depends on the goal: if the job is explaining which factor drives the outcome, this is a real problem, and I would drop one of the redundant features, combine them, or use a regularised model like Ridge regression, which handles correlated features more gracefully by shrinking coefficients instead of letting them swing to extreme, unstable values. If the job is purely predictive accuracy on data that resembles the training distribution, collinearity by itself often does not hurt performance much, since the model still captures the combined effect correctly — the real danger is that if the correlation between those features breaks down in future data in a way it never did during training, a model that leaned on that correlation can degrade sharply.

A coefficient from a fitted regression is an association conditional on whatever other features happen to be in the model, on this particular dataset — by itself it is not proof of causation, no matter how large, stable, or statistically significant it looks. The honest process is to ask whether a plausible confounder could produce the same pattern without a causal link, and in general that is a real risk: a feature like restaurant popularity correlating with delivery time could easily be driven by understaffed kitchens during that restaurant's busiest hours rather than popularity itself causing delay. Distance is actually one of the more defensible cases for a near-causal read, since travelling a longer physical distance mechanically takes more time regardless of any other variable, and that mechanical story exists independently of the regression. But that confidence comes from domain reasoning about how delivery physically works, not from the coefficient itself — to establish causation rigorously in the general case, you would want a randomised experiment or a dedicated causal inference method, not just a large coefficient in an observational dataset.