Linear regression

The exact answer

You don't have to try lines one by one. For linear regression there's an exact recipe that finds the best weights in one go. NumPy's np.linalg.lstsq runs it for you (the name is short for "least squares"). On the same three points it finds the best line, $\hat{y} = 1.5x + 0.667$:

import numpy as np

x = np.array([1.0, 2.0, 3.0])
y = np.array([2.0, 4.0, 5.0])
A = np.column_stack([x, np.ones_like(x)])
(w, b), *_ = np.linalg.lstsq(A, y, rcond=None)
print(round(w, 3), round(b, 3))
print(round(np.sum((y - (w * x + b)) ** 2), 3))

It prints 1.5 and 0.667, then the total of 0.167. That's smaller than line B's 0.5. No other straight line gets a smaller total for these three points.

The np.column_stack line builds a table with two columns: the feature, and a column of ones. lstsq finds one weight per column. The weight on the column of ones is multiplied by 1 every time, so it acts as the intercept.

import numpy as np
X = np.array([[1.0, 1.0], [2.0, 1.0], [3.0, 1.0]])
y = np.array([2.0, 4.0, 5.0])
print((np.linalg.inv(X.T @ X) @ X.T @ y).round(3))

Real data: orders and lifetime value

Now the same recipe on the practice shop's 600 customers. A customer's lifetime value is everything they've spent at the shop so far, in dollars. The code predicts it from one feature: how many orders the customer placed last year. Then it draws the customers as dots and the line through them.

import numpy as np
import pandas as pd
import requests
import matplotlib.pyplot as plt

df = pd.DataFrame(requests.get("https://api.gradient.test/v1/customers").json()["data"])
x = df["orders_last_year"].to_numpy(dtype=float)
y = df["lifetime_value"].to_numpy()

A = np.column_stack([x, np.ones_like(x)])          # a column of ones gives the intercept
(w, b), *_ = np.linalg.lstsq(A, y, rcond=None)
print(f"lifetime value ≈ {w:.1f} × orders + {b:.1f}")

plt.scatter(x, y, s=10, alpha=0.5)
xs = np.linspace(0, x.max(), 50)
plt.plot(xs, w * xs + b, color="C1")
plt.xlabel("orders last year")
plt.ylabel("lifetime value ($)")

It prints lifetime value ≈ 77.0 × orders + 302.5. Each extra order goes with about $77 more lifetime value. A customer with no orders last year is predicted at about $302.

Several features, with scikit-learn

scikit-learn is Python's standard machine-learning library. All its models work the same way:

• fit(X, y) learns the weights from data;
• predict(X) makes predictions;
• score(X, y) measures how good they are.

For linear regression, the learned weights end up in coef_, short for coefficients (another word for weights). The intercept ends up in intercept_.

The code below uses four features. train_test_split sets aside a quarter of the customers as a test set. The model never sees them while it learns, so they give an honest check. Setting random_state=0 makes the shuffle the same on every run.

import pandas as pd
import requests
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

df = pd.DataFrame(requests.get("https://api.gradient.test/v1/customers").json()["data"])
features = ["orders_last_year", "avg_order_value", "tenure_months", "trail_club"]
X_train, X_test, y_train, y_test = train_test_split(df[features], df["lifetime_value"], test_size=0.25, random_state=0)

model = LinearRegression().fit(X_train, y_train)
for name, coef in zip(features, model.coef_):
    print(f"{name:<18} {coef:8.2f}")
print(f"intercept          {model.intercept_:8.2f}")
print("R² on held-out customers:", round(model.score(X_test, y_test), 3))

Here's how to read the weights. Each one assumes the other three features stay fixed.

• orders_last_year, about 127. Each extra order adds roughly $127 of predicted lifetime value.
• avg_order_value, about 13.5. Each extra dollar on the average order adds about $13. That extra dollar is paid on every order, which is why it's worth more than $1.
• tenure_months, about 21. Each extra month as a customer adds about $21.
• trail_club, about −7. This is 1 for members of the shop's loyalty program and 0 for everyone else. Its weight comes out small and slightly negative. Once you know a customer's orders and order value, membership adds little extra information.

Notice that an order was worth about $77 in the one-feature model, and about $127 here. A weight always depends on which other features are in the model, because it's measured with them held fixed.

The intercept is about −1,496. That's the prediction for a customer with every feature at 0, and no real customer looks like that: nobody has a $0 average order. So it's only the starting point of the sum. The last line, R², scores the model on the test customers. It comes out around 0.82. The next section explains what that number means.

R²: better than guessing the average?

Before you trust a model, compare it with the laziest guess there is: predict the average every time. The shop's customers average about $776 of lifetime value. So the lazy guess says $776 for everyone. A useful model should miss by much less than the lazy guess does.

R² (said "R squared") measures how much less:

The top of the fraction is the model's sum of squared errors, from earlier. The bottom is the same total for the lazy guess. In it, $\bar{y}$, said "y-bar", is the average of the true values.

In plain words: R² is the share of the lazy guess's squared error that the model gets rid of. An R² of 0.8 means the model removes 80% of it.

A worked example. The true values are 2, 4 and 6, so the average is 4. The model predicts 3, 4 and 5.

• The model misses by −1, 0 and 1. Square and add: 1 + 0 + 1 = 2.
• The lazy guess of 4 misses by −2, 0 and 2. Square and add: 4 + 0 + 4 = 8.
• R² = 1 − 2/8 = 0.75. The model removes three quarters of the lazy guess's error.

y = [2, 4, 6]            # what really happened
y_hat = [3, 4, 5]        # the model's predictions
mean = sum(y) / len(y)   # 4.0

model_error = sum((a - p) ** 2 for a, p in zip(y, y_hat))
baseline_error = sum((a - mean) ** 2 for a in y)
print(model_error, baseline_error, 1 - model_error / baseline_error)

It prints 2, 8.0 and 0.75. Here's how to read any R²:

• 1 is perfect: the model never misses.
• 0 means the model is no better than guessing the average.
• Below 0 can happen on test data. It means the model does worse than the average would.

The customer model scored about 0.82 on its test customers. It removes about 82% of the lazy guess's squared error.

y = [1, 2, 3]
y_hat = [1, 2, 4]
mean = sum(y) / len(y)
top = sum((a - p) ** 2 for a, p in zip(y, y_hat))
bottom = sum((a - mean) ** 2 for a in y)
print(1 - top / bottom)

Look at the residuals

R² is a single number, and a single number can hide a pattern. So also look at the residuals: the misses $y - \hat{y}$ from earlier. Plot each residual against its prediction.

• If the dots form a shapeless cloud around zero, the straight line caught the pattern. What's left is noise.
• If they form a curve or a funnel, the model missed something.

Think of a tailor checking a suit. Small, random gaps here and there are fine. If the sleeves are too short on every customer, that pattern is the real problem.

For the shop's customers, a likely culprit is that lifetime value behaves like a product: roughly orders × order value. A linear model only adds things up, so it can't make a product on its own. Adding that product as an extra feature, or taking logs, could fit better.

import numpy as np

spend = np.array([1.0, 2.0, 3.0, 4.0, 5.0])      # weekly ad spend, in $100s
signups = np.array([2.0, 5.0, 6.0, 8.0, 11.0])   # signups that week

A = np.column_stack([spend, np.ones_like(spend)])
(w, b), *_ = np.linalg.lstsq(A, signups, rcond=None)
predicted = w * spend + b
residuals = signups - predicted
r2 = 1 - np.sum(residuals ** 2) / np.sum((signups - signups.mean()) ** 2)

print(f"signups ≈ {w:.2f} × spend + {b:.2f}")
print("residuals:", residuals.round(2))
print("R²:", round(r2, 3))
print("prediction for $600:", round(w * 6 + b, 2))

Limits to keep in mind

• It only draws straight lines. It can't bend unless you give it bent features, like squares, logs or products. The other option is a more flexible model.
• Outliers pull hard. Misses are squared, so one far-off point can drag the whole line toward it. The next lesson looks at this.
• Don't extrapolate. That means predicting far outside the data you learned from. The shop's customers placed 0 to 18 orders last year. The line says nothing reliable about one with 200.
• Weights aren't causes. "Members are worth more" doesn't mean membership makes them worth more.
• Compare weights only after scaling. A weight per dollar and a weight per order use different units. You can't compare them until the features are put on the same scale (standardized, from the data preparation lesson).