Logistic regression

Training: log loss and gradient descent

You have the model's shape. Training finds the weights. Linear regression had an exact recipe for its best weights. Logistic regression has none. So you use gradient descent, from earlier in this phase: start somewhere, then take many small steps downhill on the loss.

The loss: log loss

The loss is log loss, from the lesson on loss functions. It scores each probability against what really happened. Write the truth as $y$: 1 if the customer left, 0 if they stayed.

• If the customer left, the loss is $-\log p$.
• If the customer stayed, the loss is $-\log(1 - p)$.

Here log is the natural log. The log loss for a whole dataset is the average over all customers.

In plain words: the loss is small when you gave a high probability to what really happened. It's large when you were confident and wrong.

Say the model gives a customer 0.9. If they left, the loss is $-\log 0.9 \approx 0.105$. If they stayed, it's $-\log 0.1 \approx 2.303$. Being sure and wrong costs about 22 times more.

import math

p = 0.9                              # the model's churn probability
print(round(-math.log(p), 3))        # loss if the customer left
print(round(-math.log(1 - p), 3))    # loss if they stayed

It prints 0.105 and then 2.303.

The step: a familiar gradient

Each training step nudges every weight. The gradient says which way and how far. For logistic regression it has a tidy form:

Take the symbols one at a time:

• $L$ is the log loss.
• $\nabla_{\mathbf{w}} L$ is its gradient. The upside-down triangle is called nabla. The gradient holds one number per weight, saying which way is uphill.
• $n$ is the number of customers.
• $X$ is the table of features, one row per customer.
• $X^\top$ is that table flipped, from the Phase 2 lesson on transpose.
• $\mathbf{p}$ is the list of predicted probabilities.
• $\mathbf{y}$ is the list of true answers.

In plain words: for each customer, work out the miss: predicted probability minus truth. For each weight, multiply every customer's miss by their value of that feature, then average. That says which way is uphill for that weight, so you step the other way.

It has the same shape as linear regression's gradient: the miss, weighted by the features. (The softmax and cross-entropy gradient in Phase 2 was the version of this for many classes.)

Here's one step by hand. There are two customers and one feature, plus a column of ones for the bias. Customer A has feature 1 and left. Customer B has feature −1 and stayed. All weights start at 0.

• Both scores are 0, so both probabilities are 0.5.
• The misses: A gets 0.5 − 1 = −0.5. B gets 0.5 − 0 = 0.5.
• The feature's gradient: (−0.5 × 1 + 0.5 × −1) / 2 = −0.5.
• The bias's gradient: its "feature" is always 1, so (−0.5 + 0.5) / 2 = 0.
• Step with a learning rate of 0.5: the feature's weight becomes 0 − 0.5 × (−0.5) = 0.25.

The weight went up. Now a higher feature means a higher chance of leaving. That fits the data: the customer who left had the higher feature.

import numpy as np

X = np.array([[1.0, 1.0], [-1.0, 1.0]])   # feature, then a column of ones for the bias
y = np.array([1, 0])                      # A left, B stayed
w = np.zeros(2)

p = 1 / (1 + np.exp(-(X @ w)))            # both 0.5
gradient = X.T @ (p - y) / len(y)
w -= 0.5 * gradient                       # learning rate 0.5
print(gradient, w)

It prints the gradient, −0.5 and 0, then the new weights, 0.25 and 0. In NumPy, X.T is the flipped table and @ multiplies tables.

import math

y, h = 1, 0.000001          # a customer who left; a tiny nudge

def loss(z):
    p = 1 / (1 + math.exp(-z))
    return -(y * math.log(p) + (1 - y) * math.log(1 - p))

z = 0.5
slope = (loss(z + h) - loss(z)) / h
p = 1 / (1 + math.exp(-z))
print(round(slope, 4), round(p - y, 4))

Training on real customers

Now the full loop on the practice shop's 600 customers. It uses five features. Each one is standardized, from the data preparation lesson. That shifts and scales it so its average is 0 and its typical spread, the standard deviation, is 1. A column of ones carries the bias. The loop takes 2,000 steps with a learning rate of 0.5.

import numpy as np
import pandas as pd
import requests

df = pd.DataFrame(requests.get("https://api.gradient.test/v1/customers").json()["data"])
features = ["orders_last_year", "days_since_last_order", "support_tickets", "trail_club", "tenure_months"]
X_raw = df[features].to_numpy(dtype=float)
X = np.column_stack([(X_raw - X_raw.mean(axis=0)) / X_raw.std(axis=0), np.ones(len(df))])
y = df["churned"].to_numpy()

sigmoid = lambda z: 1 / (1 + np.exp(-z))
w = np.zeros(X.shape[1])
for step in range(2000):
    p = sigmoid(X @ w)
    w -= 0.5 * X.T @ (p - y) / len(y)

p = sigmoid(X @ w)
log_loss = -np.mean(y * np.log(p) + (1 - y) * np.log(1 - p))
print("log loss:", round(log_loss, 3), " accuracy at 0.5:", round(np.mean((p > 0.5) == y), 3))
for name, coef in zip(features + ["bias"], w):
    print(f"{name:<22} {coef:+.2f}")

It prints a log loss of 0.374. Then comes an accuracy of 0.838: with a cut-off of 0.5, the model labels about 84% of customers correctly.

Be careful with that number. Only about 19% of these customers churned. So guessing "stays" for everyone would already score about 81%. The evaluation lessons later in this phase show better ways to judge a yes/no model.

Below that, it prints one weight per feature, and the bias:

• orders_last_year, −0.82. More orders, less churn.
• days_since_last_order, +0.72. The longer since the last order, the more churn.
• support_tickets, +0.34. More complaints, more churn.
• trail_club, −0.38. Loyalty program members churn less.
• tenure_months, −0.33. Longer-standing customers churn less.
• bias, −1.93. The score of a customer who is average on every feature.

p, y, x = 0.5, 1, 2
w = 0.0
w -= 0.1 * (p - y) * x      # learning rate 0.1
print(w)

Reading the coefficients

The weights are also called coefficients. To read them, you need one idea from betting: odds.

Odds

Odds of "3 to 1" mean something is three times as likely to happen as not. The odds of an event are $p / (1 - p)$: the chance it happens, divided by the chance it doesn't.

In plain words: odds compare "yes" with "no". Odds of 3 mean "yes" is three times as likely as "no". Odds of 1 mean an even split.

With small numbers: a churn probability of 0.75 gives odds of 0.75 / 0.25 = 3. That's 3 to 1.

Log-odds: what the score really is

Here is the key fact. The score $z$ that the model adds up is the log of the odds:

In plain words: take the odds, then take their natural log. You get back the score the model added up. That's why the score is also called the log-odds.

Try it with a score of 1.1. The sigmoid gives a probability of 0.75. The odds are 3.0, and the log of 3 is about 1.1, the score you started with.

import math

z = 1.1
p = 1 / (1 + math.exp(-z))
odds = p / (1 - p)
print(round(p, 2), round(odds, 2), round(math.log(odds), 2))

It prints 0.75, 3.0 and 1.1.

What one coefficient does

Each coefficient adds to the log-odds. Adding to a log is the same as multiplying the thing inside it. You saw this in the lesson on exponentials: logs turn products into sums. So when a feature goes up by 1, its coefficient $w$ multiplies the odds by:

In plain words: raise $e$ to the coefficient. That's how many times the odds grow when the feature goes up by one unit, with the other features held fixed.

This multiplier is called the odds ratio. A positive coefficient gives an odds ratio above 1, so the odds grow. A negative one gives an odds ratio below 1, so the odds shrink. A coefficient of 0 gives exactly 1: no change.

Here's a worked example from the shop. The features were standardized, so "one unit" means one standard deviation. For days since the last order, that's about 94 days, roughly three months. Its coefficient is about +0.72.

• The odds ratio: $e^{0.72} \approx 2.05$. Three more quiet months roughly double the odds.
• Take a customer with a 20% churn chance. Their odds are 0.2 / 0.8 = 0.25, or 1 to 4.
• Three more quiet months: 0.25 × 2.05 ≈ 0.51.
• Turn odds back into a probability with odds / (1 + odds): 0.51 / 1.51 ≈ 0.34.

So the churn chance goes from 20% to about 34%. The odds doubled, but the probability didn't. A coefficient changes the odds by a fixed multiple. How much the probability moves depends on where it started.

import math

odds = 0.25                          # a 20% churn chance: 1 to 4
new_odds = odds * math.exp(0.72)     # three more quiet months
print(round(math.exp(0.72), 2), round(new_odds, 2), round(new_odds / (1 + new_odds), 2))

It prints 2.05, 0.51 and 0.34.

The negative coefficients work the same way in reverse. Orders have about −0.82, and $e^{-0.82} \approx 0.44$. So one standard deviation more orders cuts the odds of churning to less than half. Membership and tenure lower the odds too.

Every sign matches common sense. That's a useful sanity check on any model. Comparing the sizes of these coefficients is fair only because the features were standardized first. On raw units, a weight per day and a weight per order can't be compared.

p = 0.9
print(round(p / (1 - p), 1))

The same model in scikit-learn

In practice you'd use scikit-learn rather than write the loop yourself. The code below joins two steps with make_pipeline. StandardScaler standardizes the features, then LogisticRegression trains the model. train_test_split holds back a quarter of the customers for testing. Setting stratify keeps the same share of churners in both parts.

import numpy as np
import pandas as pd
import requests
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

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

model = make_pipeline(StandardScaler(), LogisticRegression()).fit(X_train, y_train)
probs = model.predict_proba(X_test)[:, 1]          # column 1 = P(churn)
print("first five churn probabilities:", probs[:5].round(2))
for threshold in [0.5, 0.3, 0.2]:
    flagged = probs >= threshold
    caught = (flagged & (y_test == 1)).sum()
    print(f"threshold {threshold}: flag {flagged.sum():>3} customers, catching {caught} of {(y_test == 1).sum()} churners")

predict_proba returns one column per class. Column 0 is the chance of staying. Column 1 is the chance of churning. Each row adds up to 1. The shortcut predict turns these into labels with a cut-off of 0.5.

The numbers won't match the from-scratch model exactly. This one trains on only three quarters of the customers. scikit-learn also adds a small penalty by default that keeps weights from growing large. That's called regularization, and it comes later in this phase.

The loop at the end tries three cut-offs. A lower cut-off flags more customers. It never catches fewer of the real churners, and usually catches more. The price is that more of the flagged customers would have stayed anyway.

Choosing the threshold

The cut-off is called the threshold. Above it, you act. 0.5 is a default, not a law. The right threshold depends on what each mistake costs.

Say a retention offer costs $5, and a lost customer costs $200. Sending the offer pays off when the expected saving beats the cost. The expected saving is the churn probability times $200. So the offer is worth it when $p \times 200 > 5$, which means $p > 5 / 200 = 0.025$.

In plain words: act when the probability times what you'd save is bigger than what acting costs. The break-even threshold is the cost divided by the value.

offer_cost = 5         # dollars per offer
customer_value = 200   # dollars lost when a customer leaves
print(offer_cost / customer_value)

It prints 0.025, far below 0.5. Flagging more customers is worth it, even though many of them would have stayed anyway.

This is a simplification. It assumes the offer always keeps the customer. If it works only some of the time, the threshold goes up. The evaluation lessons later in this phase (precision, recall and ROC curves) give you the tools to pick a threshold deliberately.

import numpy as np

months = np.array([0, 1, 2, 3, 4, 5, 6, 7], dtype=float)   # months since last order
churned = np.array([0, 0, 1, 0, 0, 1, 1, 1])                # 1 = the customer left

X = np.column_stack([months, np.ones_like(months)])         # feature, then a column of ones
sigmoid = lambda z: 1 / (1 + np.exp(-z))

w = np.zeros(2)
for step in range(5000):
    p = sigmoid(X @ w)
    w -= 0.1 * X.T @ (p - churned) / len(churned)

print("weight:", round(w[0], 2), " bias:", round(w[1], 2))
print("odds ratio per month:", round(np.exp(w[0]), 2))
print("probabilities:", sigmoid(X @ w).round(2))
print("boundary at month:", round(-w[1] / w[0], 2))