3.5 - Noah's Notes - Logistic regression from scratch¶
Four months after the supply-prediction project, Minermont Hospital faces a new question: will a patient attend their appointment? Across the hospital about 18% of appointments are missed, and in some specialties up to 27%. Physiotherapy has the highest rate, so the pilot starts there.
This time we do not predict a quantity, like the masks of chapter 2, but a yes/no answer: a binary classification problem. Logistic regression keeps the linear score of linear regression, $z = w \cdot x + b$, and passes it through the sigmoid function, which turns any number into a probability between 0 and 1:
$$\hat{p} = \sigma(z) = \frac{1}{1 + e^{-z}}$$
If $\hat{p}$ is above a threshold (usually 0.5), we predict that the patient will not attend.
In this notebook we build it from scratch with NumPy, just as we did with linear regression:
- Import the libraries.
- Create a synthetic dataset of Physiotherapy appointments.
- Prepare the data: training and test sets, and feature scaling.
- Define the sigmoid function.
- Define the cost function (cross-entropy, or log loss).
- Train the model with gradient descent.
- Evaluate it: confusion matrix, sensitivity, specificity, PPV and NPV.
- Choose the decision threshold.
- Interpret the coefficients.
- Predict for new patients.
Step 1: Import the libraries¶
numpy for the calculations, pandas to handle the data as a table and matplotlib for the plots. We will only use scikit-learn for one utility: splitting the data into training and test sets. The model itself is ours.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split # only to split the data
np.set_printoptions(precision=3, suppress=True)
print("Libraries loaded.")
Libraries loaded.
Step 2: A synthetic dataset of appointments¶
While Martha and Ethan gather the hospital's real attendance records, I am creating a synthetic dataset, as we did in chapter 2. Each row is one Physiotherapy appointment, described by the variables the team considered relevant:
| Column | Meaning |
|---|---|
age |
Age of the patient (18–85 years) |
days_since_last_visit |
Days since the patient's previous appointment |
distance_km |
Distance from home to the hospital, in km (Minermont serves a rural area) |
hour |
Time of the appointment (8 to 19 h) |
first_visit |
1 if it is a first visit, 0 if it is a follow-up |
previous_no_shows |
Number of appointments the patient has missed before |
no_show |
Target: 1 if the patient did not attend, 0 if they did |
The positive class (1) is the no-show, because it is the event we want to detect. Behind the data there is a hidden rule that we only use to generate the appointments; the model will have to discover the pattern from the examples.
def generate_appointments(n, seed):
"""Synthetic Physiotherapy appointments at Minermont: one row per appointment."""
rng = np.random.default_rng(seed)
age = rng.integers(18, 86, n) # years
days = rng.integers(1, 181, n) # days since the previous appointment
distance = np.round(rng.gamma(2.0, 6.0, n).clip(0.5, 60), 1) # km to the hospital (a rural area)
hour = rng.integers(8, 20, n) # appointment time, from 8:00 to 19:00
first = rng.binomial(1, 0.35, n) # 1 = first visit, 0 = follow-up
previous = rng.poisson(0.6, n).clip(0, 6) # no-shows in the patient's history
late = np.where(hour >= 17, hour - 16, 0) # late-afternoon slots are missed more
# The hidden rule behind the data. We only use it to generate the appointments:
# the model never sees it and has to discover the pattern from the examples.
z = (-1.55 - 0.069 * (age - 50) + 0.026 * (days - 90) + 0.11 * (distance - 12)
+ 1.5 * (previous - 0.6) + 0.56 * (late - 0.75) + 0.63 * (first - 0.35))
no_show = rng.binomial(1, 1 / (1 + np.exp(-z))) # 1 = did not attend
return pd.DataFrame({
'age': age, 'days_since_last_visit': days, 'distance_km': distance, 'hour': hour,
'first_visit': first, 'previous_no_shows': previous, 'no_show': no_show,
})
df = generate_appointments(1200, seed=42)
print(f"Appointments: {len(df)}")
print(f"No-show rate: {df['no_show'].mean():.1%}")
df.head()
Appointments: 1200 No-show rate: 26.2%
| age | days_since_last_visit | distance_km | hour | first_visit | previous_no_shows | no_show | |
|---|---|---|---|---|---|---|---|
| 0 | 24 | 159 | 4.1 | 11 | 0 | 2 | 1 |
| 1 | 70 | 86 | 43.7 | 11 | 1 | 0 | 0 |
| 2 | 62 | 178 | 7.1 | 18 | 0 | 0 | 0 |
| 3 | 47 | 138 | 3.5 | 8 | 0 | 0 | 0 |
| 4 | 47 | 118 | 9.1 | 8 | 0 | 3 | 1 |
Before modelling, a quick look at the data. How does the no-show rate change with some of the variables?
fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)
groups = [
(pd.cut(df['age'], [17, 35, 50, 65, 85], labels=['18-35', '36-50', '51-65', '66-85']), 'Age'),
(pd.cut(df['distance_km'], [0, 5, 15, 30, 60], labels=['< 5', '5-15', '15-30', '> 30']), 'Distance to the hospital (km)'),
(df['previous_no_shows'].clip(upper=3).map({0: '0', 1: '1', 2: '2', 3: '3+'}), 'Previous no-shows'),
]
for ax, (group, title) in zip(axes, groups):
df.groupby(group, observed=True)['no_show'].mean().plot.bar(ax=ax, color='tab:blue')
ax.set_title(title)
ax.set_xlabel('')
ax.tick_params(axis='x', rotation=0)
axes[0].set_ylabel('No-show rate')
plt.tight_layout()
plt.show()
Younger patients, patients who live far away and patients who have missed appointments before are absent more often. No single variable decides on its own: the model will have to combine them.
Step 3: Prepare the data¶
We split the appointments into a training set (75%), which the model learns from, and a test set (25%), which we keep aside to evaluate it on appointments it has never seen. stratify=y keeps the same proportion of no-shows in both sets.
Then we scale the features: we subtract the mean and divide by the standard deviation, so that every variable is measured in "standard deviations". Without this, a variable in km and another in days would have very different sizes and gradient descent would converge slowly. The mean and standard deviation are computed only on the training set and applied to both sets: the test set must not influence training in any way.
feature_names = ['age', 'days_since_last_visit', 'distance_km', 'hour', 'first_visit', 'previous_no_shows']
X = df[feature_names].to_numpy(dtype=float)
y = df['no_show'].to_numpy()
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, random_state=42, stratify=y)
# Scaling with the statistics of the TRAINING set only
mean = X_train.mean(axis=0)
std = X_train.std(axis=0)
X_train_s = (X_train - mean) / std
X_test_s = (X_test - mean) / std
print(f"Training set: {len(y_train)} appointments ({y_train.mean():.1%} no-shows)")
print(f"Test set: {len(y_test)} appointments ({y_test.mean():.1%} no-shows)")
Training set: 900 appointments (26.2% no-shows) Test set: 300 appointments (26.3% no-shows)
Step 4: The sigmoid function¶
The sigmoid $\sigma(z) = \frac{1}{1 + e^{-z}}$ turns the linear score $z$ into a probability: large positive scores give probabilities close to 1, large negative scores give probabilities close to 0, and $z = 0$ gives exactly 0.5, the default threshold.
def sigmoid(z):
"""Maps any real number (or array) to a value between 0 and 1."""
return 1 / (1 + np.exp(-z))
for z in [-4, -2, 0, 2, 4]:
print(f"sigmoid({z:>2}) = {sigmoid(z):.3f}")
z = np.linspace(-8, 8, 200)
plt.figure(figsize=(8, 4))
plt.plot(z, sigmoid(z), color='tab:orange', linewidth=3)
plt.axhline(0.5, color='gray', linestyle='--', label='Threshold 0.5')
plt.xlabel('Linear score z = w·x + b')
plt.ylabel('Probability of no-show')
plt.title('The sigmoid function')
plt.legend()
plt.grid(True)
plt.show()
sigmoid(-4) = 0.018 sigmoid(-2) = 0.119 sigmoid( 0) = 0.500 sigmoid( 2) = 0.881 sigmoid( 4) = 0.982
Step 5: The cost function: cross-entropy (log loss)¶
To measure how wrong the model is we cannot use the mean squared error of chapter 2: with a sigmoid it produces a cost surface with many valleys. Logistic regression uses the cross-entropy, or log loss:
$$J(w, b) = -\frac{1}{m} \sum_{i=1}^{m} \Big[ y_i \log(\hat{p}_i) + (1 - y_i) \log(1 - \hat{p}_i) \Big]$$
For a patient who did not attend ($y = 1$) only the first term counts, $-\log(\hat{p})$: it is almost 0 if the model gave a probability close to 1, and it grows without limit as that probability approaches 0. Confident mistakes are punished hard.
def compute_cost(X, y, w, b):
"""Cross-entropy (log loss) of the model with weights w and bias b."""
p = sigmoid(X @ w + b)
p = np.clip(p, 1e-15, 1 - 1e-15) # avoids log(0)
return -np.mean(y * np.log(p) + (1 - y) * np.log(1 - p))
# The cost of a single patient who did NOT attend (y = 1):
for p in [0.9, 0.5, 0.1]:
print(f"The model predicted p = {p}: cost = {-np.log(p):.3f}")
The model predicted p = 0.9: cost = 0.105 The model predicted p = 0.5: cost = 0.693 The model predicted p = 0.1: cost = 2.303
Step 6: Training with gradient descent¶
The gradients of the cross-entropy have exactly the same form as in linear regression; the only difference is that the prediction now goes through the sigmoid:
$$\frac{\partial J}{\partial w} = \frac{1}{m} X^\top (\hat{p} - y) \qquad \frac{\partial J}{\partial b} = \frac{1}{m} \sum_{i=1}^{m} (\hat{p}_i - y_i)$$
and the update rule is the one we already know: $w := w - \alpha \frac{\partial J}{\partial w}$ and $b := b - \alpha \frac{\partial J}{\partial b}$, where $\alpha$ is the learning rate. Because the features are scaled, a learning rate of 0.1 converges quickly.
def train_logistic(X, y, learning_rate=0.1, num_iterations=3000):
"""Batch gradient descent on the cross-entropy. Returns w, b and the cost of every iteration."""
m, n = X.shape
w = np.zeros(n)
b = 0.0
costs = []
for i in range(num_iterations):
p = sigmoid(X @ w + b) # a) predict with the current parameters
dw = (1 / m) * X.T @ (p - y) # b) gradients
db = (1 / m) * np.sum(p - y)
w = w - learning_rate * dw # c) one step downhill
b = b - learning_rate * db
costs.append(compute_cost(X, y, w, b))
if i % 500 == 0:
print(f"Iteration {i:4d}: cost = {costs[-1]:.4f}")
return w, b, costs
w, b, costs = train_logistic(X_train_s, y_train)
print(f"\nFinal cost: {costs[-1]:.4f}")
print("Learned weights (w):", w)
print(f"Learned bias (b): {b:.3f}")
plt.figure(figsize=(8, 4))
plt.plot(costs)
plt.xlabel('Iteration')
plt.ylabel('Cost (log loss)')
plt.title('Learning curve')
plt.grid(True)
plt.show()
Iteration 0: cost = 0.6807 Iteration 500: cost = 0.3523 Iteration 1000: cost = 0.3511 Iteration 1500: cost = 0.3511 Iteration 2000: cost = 0.3511 Iteration 2500: cost = 0.3511 Final cost: 0.3511 Learned weights (w): [-1.269 1.169 1.02 0.365 0.286 1.026] Learned bias (b): -1.808
The cost falls quickly and then flattens: the model has converged.
Step 7: Evaluating the model¶
With a threshold we turn each probability into a prediction, and comparing the predictions with what really happened gives four possible outcomes, the confusion matrix:
| Predicted: will not attend | Predicted: will attend | |
|---|---|---|
| Did not attend | True positive (TP) | False negative (FN) |
| Attended | False positive (FP) | True negative (TN) |
From it we compute the metrics used in medicine:
- Accuracy: $\frac{TP + TN}{\text{total}}$, the share of correct predictions.
- Sensitivity (recall): $\frac{TP}{TP + FN}$, the share of no-shows that the model detects.
- Specificity: $\frac{TN}{TN + FP}$, the share of attendances that the model recognises.
- Positive predictive value (PPV): $\frac{TP}{TP + FP}$, when the model says "will not attend", how often it is right.
- Negative predictive value (NPV): $\frac{TN}{TN + FN}$, when the model says "will attend", how often it is right.
A useful reference: a "model" that always predicts that every patient will attend already has an accuracy equal to the attendance rate, without detecting a single no-show. Our model has to beat that.
def predict_proba(X, w, b):
return sigmoid(X @ w + b)
def evaluate(y_true, p, threshold=0.5):
"""Confusion matrix and clinical metrics for a given threshold."""
y_pred = (p >= threshold).astype(int)
tp = np.sum((y_pred == 1) & (y_true == 1))
fn = np.sum((y_pred == 0) & (y_true == 1))
fp = np.sum((y_pred == 1) & (y_true == 0))
tn = np.sum((y_pred == 0) & (y_true == 0))
return {
'threshold': threshold, 'TP': tp, 'FN': fn, 'FP': fp, 'TN': tn,
'accuracy': (tp + tn) / len(y_true),
'sensitivity': tp / (tp + fn),
'specificity': tn / (tn + fp),
'PPV': tp / (tp + fp),
'NPV': tn / (tn + fn),
}
p_test = predict_proba(X_test_s, w, b)
r = evaluate(y_test, p_test, threshold=0.5)
print("Confusion matrix (test set, threshold 0.5):")
print(f" Predicted no-show Predicted attends")
print(f" Did not attend TP = {r['TP']:<14} FN = {r['FN']}")
print(f" Attended FP = {r['FP']:<14} TN = {r['TN']}")
print()
for k in ['accuracy', 'sensitivity', 'specificity', 'PPV', 'NPV']:
print(f"{k:>12}: {r[k]:.3f}")
print(f"\nBaseline (always predict 'attends'): accuracy = {1 - y_test.mean():.3f}")
Confusion matrix (test set, threshold 0.5):
Predicted no-show Predicted attends
Did not attend TP = 49 FN = 30
Attended FP = 22 TN = 199
accuracy: 0.827
sensitivity: 0.620
specificity: 0.900
PPV: 0.690
NPV: 0.869
Baseline (always predict 'attends'): accuracy = 0.737
The model beats the baseline, and it is very reliable when it says that a patient will attend (high NPV). But with the default threshold it only detects a part of the no-shows (the sensitivity): many absences go unnoticed.
Step 8: Choosing the threshold¶
The 0.5 threshold is only a convention. Lowering it makes the model flag more patients as likely no-shows: sensitivity goes up (fewer absences slip through), but so do the false alarms, so specificity and PPV go down. Which balance is better is not a mathematical question but a clinical and organisational one: here a false alarm costs a reminder phone call, while a missed no-show costs an empty slot on a waiting list.
table = pd.DataFrame([evaluate(y_test, p_test, t) for t in [0.5, 0.4, 0.3, 0.2]])
table.round(3)
| threshold | TP | FN | FP | TN | accuracy | sensitivity | specificity | PPV | NPV | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.5 | 49 | 30 | 22 | 199 | 0.827 | 0.620 | 0.900 | 0.690 | 0.869 |
| 1 | 0.4 | 54 | 25 | 32 | 189 | 0.810 | 0.684 | 0.855 | 0.628 | 0.883 |
| 2 | 0.3 | 61 | 18 | 46 | 175 | 0.787 | 0.772 | 0.792 | 0.570 | 0.907 |
| 3 | 0.2 | 68 | 11 | 64 | 157 | 0.750 | 0.861 | 0.710 | 0.515 | 0.935 |
With a threshold of 0.3 the model detects many more no-shows, at the cost of more reminder calls to patients who would have come anyway.
Step 9: Interpreting the coefficients¶
Because the features are scaled, the weights are comparable with each other: each one says how much the log-odds of a no-show change when that variable increases by one standard deviation, with the others fixed. Its exponential, $e^{w}$, is the odds ratio: a value of 2 doubles the odds of a no-show, a value of 0.5 halves them.
coefficients = pd.DataFrame({'weight': w, 'odds ratio': np.exp(w)}, index=feature_names)
coefficients = coefficients.reindex(coefficients['weight'].abs().sort_values(ascending=False).index)
print(coefficients.round(3))
plt.figure(figsize=(8, 4))
plt.barh(coefficients.index[::-1], coefficients['weight'][::-1],
color=['tab:red' if v > 0 else 'tab:blue' for v in coefficients['weight'][::-1]])
plt.axvline(0, color='black', linewidth=1)
plt.xlabel('Weight (change in log-odds per standard deviation)')
plt.title('What makes a no-show more likely (red) or less likely (blue)')
plt.tight_layout()
plt.show()
weight odds ratio age -1.269 0.281 days_since_last_visit 1.169 3.217 previous_no_shows 1.026 2.791 distance_km 1.020 2.774 hour 0.365 1.440 first_visit 0.286 1.332
The model has found the pattern the team suspected, and the same three variables the book highlights come out on top. Age is the strongest signal: each standard deviation of age (about 20 years) multiplies the odds of a no-show by about 0.28, so older patients miss far fewer appointments. The days since the last visit (×3.2), the previous no-shows (×2.8) and the distance to the hospital (×2.8) raise the risk. The time of the appointment (×1.4) and the type of visit (×1.3) matter less. These are associations learned from the data, not causes: they tell the team whom to call, not why a patient misses an appointment.
Step 10: Predicting for new patients¶
Finally, we use the trained model for two new appointments. The new data must be scaled with the same mean and standard deviation as the training set.
# age, days since last visit, distance (km), hour, first visit, previous no-shows
new_patients = np.array([
[35, 100, 18.0, 18, 1, 0], # young, long gap, lives far away, late slot, first visit
[33, 80, 14.0, 17, 1, 0], # similar, but a little closer to the usual profile
[62, 30, 6.0, 9, 0, 0], # older, recent follow-up, lives nearby, morning slot
])
p_new = predict_proba((new_patients - mean) / std, w, b)
for patient, p in zip(new_patients, p_new):
flag_05 = "flagged" if p >= 0.5 else "not flagged"
flag_03 = "flagged" if p >= 0.3 else "not flagged"
print(f"Patient aged {patient[0]:.0f}: P(no-show) = {p:5.1%} | threshold 0.5: {flag_05:<11} | threshold 0.3: {flag_03}")
Patient aged 35: P(no-show) = 57.5% | threshold 0.5: flagged | threshold 0.3: flagged Patient aged 33: P(no-show) = 36.0% | threshold 0.5: not flagged | threshold 0.3: flagged Patient aged 62: P(no-show) = 0.3% | threshold 0.5: not flagged | threshold 0.3: not flagged
The second patient shows why the threshold matters: with 0.5 the model would not flag them, with 0.3 the team would call to confirm the appointment.
Conclusion¶
We have built a logistic regression from scratch: the sigmoid turns a linear score into a probability, the cross-entropy measures the error, and gradient descent adjusts the weights exactly as it did in linear regression.
We have also seen that a classifier is not evaluated with a single number. Accuracy hides the fact that, with the default threshold, many no-shows go undetected; sensitivity, specificity, PPV and NPV tell the whole story, and the threshold lets the team choose the balance that suits the hospital.
In the next notebook, Logistic regression with scikit-learn, we will train the same model with a library and add two tools from this chapter: regularization, to avoid overfitting, and k-fold cross-validation, to choose its strength.