2.6 - Multiple linear regression: the five-variable mask model

In chapter 2, Martha extends the line of linear regression to a plane and, with more variables, to a hyperplane: $\hat{y} = w_1 x_1 + w_2 x_2 + \dots + w_n x_n + b$. Liam brings eighteen months of hospital records, Grace suggests the clinical variables that matter, and the team fits a model that predicts the masks needed each day from five variables: patients, temperature, emergencies, health alerts and weekends.

The other two notebooks of this chapter use a single variable. Here we implement the model with several, the way the book describes it. To be able to check the result, we generate the data from the equation the team obtained in the book, plus some noise, and we see whether each method recovers it:

1. Import the libraries.
2. From a line to a hyperplane: Martha's example.
3. Eighteen months of data.
4. Training and test months.
5. Gradient descent with Noah's update equations.
6. Scaling the features.
7. The exact solution: least squares and scikit-learn.
8. How good is the model?
9. When the data has errors: MSE or MAE?
10. The forecast for next week.

Step 1: Import the libraries

---

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

np.set_printoptions(precision=3, suppress=True)

Step 2: From a line to a hyperplane

---

With several variables, the prediction is a weighted sum: each variable $x_i$ is multiplied by its weight $w_i$ and the bias $b$ is added. Written with vectors, $\hat{y} = \mathbf{w} \cdot \mathbf{x} + b$, which is exactly what NumPy's @ operator computes. Martha's example in the book uses three variables: patients, temperature and the health alert level, with $\hat{y} = 0.6\,x_1 - 0.4\,x_2 + 20\,x_3 + 50$, for a day with 300 patients, 8 °C and no alert:

w = np.array([0.6, -0.4, 20.0])      # masks per patient, per degree, per alert level
b = 50.0                             # the daily baseline
x = np.array([300, 8, 0])            # 300 patients, 8 degrees, no alert

y_hat = w @ x + b
print(f"Predicted masks: {y_hat:.1f}")
print("Term by term:", w * x, "+", b)

Predicted masks: 226.8
Term by term: [180.   -3.2   0. ] + 50.0

226.8 masks: 180 for the 300 patients, −3.2 for the temperature, nothing for the alert, plus the baseline of 50. As Martha says, the decimal does not mean a fraction of a mask; it is an estimate, and in practice the centre would prepare about 227. A prediction with five variables is the same operation with longer vectors.

Step 3: Eighteen months of data

---

Liam's records cover eighteen months. We simulate them, one row per day, with the variables of the model the team fits in the book: the patients seen (fewer at weekends), the temperature (colder in winter), the emergencies (more in winter), three two-week periods of health alert, and whether the day is a weekend. The masks used are generated with the team's equation,

$$\text{masks} = 42.3 + 0.05\,\text{patients} - 0.37\,\text{temperature} + 3.2\,\text{emergencies} + 18.7\,\text{alert} + 4.1\,\text{weekend}$$

plus random noise with a standard deviation of 4 masks. The models will only see the data, never this equation, so we can check afterwards how close each one gets to it.

rng = np.random.default_rng(7)
dates = pd.date_range('2024-01-01', '2025-06-30', freq='D')
n = len(dates)
day_of_year = dates.dayofyear.values
winter = np.cos(2 * np.pi * (day_of_year - 20) / 365)          # 1 in mid-January, -1 in mid-July

weekend = (dates.dayofweek >= 5).astype(int)
temperature = 16 - 8 * winter + rng.normal(0, 2.5, n)
patients = np.where(weekend == 1, rng.normal(245, 18, n), rng.normal(295, 18, n)).round()
emergencies = rng.poisson(np.where(winter > 0.5, 13, 10))
alert = np.zeros(n, dtype=int)
for start in ['2024-01-15', '2024-12-20', '2025-02-03']:
    first = dates.get_loc(pd.Timestamp(start))
    alert[first:first + 14] = 1

book_w = np.array([0.05, -0.37, 3.2, 18.7, 4.1])
book_b = 42.3
X = np.c_[patients, temperature, emergencies, alert, weekend]
masks = X @ book_w + book_b + rng.normal(0, 4, n)

features = ['patients', 'temperature', 'emergencies', 'alert', 'weekend']
data = pd.DataFrame(X, index=dates, columns=features).assign(masks=masks)
print(data.round(1).head())
print(f"\n{n} days, {alert.sum()} of them with a health alert")

            patients  temperature  emergencies  alert  weekend  masks
2024-01-01     309.0          8.4          9.0    0.0      0.0   91.6
2024-01-02     301.0          9.1         12.0    0.0      0.0   83.4
2024-01-03     283.0          7.7         15.0    0.0      0.0  101.9
2024-01-04     332.0          6.1         16.0    0.0      0.0  109.6
2024-01-05     337.0          7.1         13.0    0.0      0.0   93.9

547 days, 42 of them with a health alert

fig, axes = plt.subplots(1, 3, figsize=(16, 4))
axes[0].plot(data.index, data['masks'], lw=0.8)
axes[0].set_title('Masks used each day')
axes[0].tick_params(axis='x', rotation=30)
axes[1].scatter(data['patients'], data['masks'], s=8, c=data['weekend'], cmap='coolwarm')
axes[1].set_xlabel('patients'); axes[1].set_title('Masks and patients (red: weekend)')
axes[2].scatter(data['temperature'], data['masks'], s=8, c=data['alert'], cmap='coolwarm')
axes[2].set_xlabel('temperature (°C)'); axes[2].set_title('Masks and temperature (red: alert)')
plt.tight_layout()
plt.show()

547 days, 42 of them with a health alert. The masks rise in both winters, when it is colder and there are more emergencies, and jump during the alerts (red in the right-hand plot). The number of patients, on its own, says little: its weight in the book's equation is small, and weekends, with fewer patients, still use as many masks. That is why the model needs all the variables at once.

Step 4: Training and test months

---

Hazel's rule in chapter 2 is to never evaluate a model on the data it was trained on. With data ordered in time, the honest split is by date: we train on the first fifteen months and keep the last three, which the model will never see, to test it. That is also how the hospital will use it: trained on the past, asked about the future.

train = dates < pd.Timestamp('2025-04-01')
X_train, y_train = X[train], masks[train]
X_test, y_test = X[~train], masks[~train]
m = len(y_train)
print(f"Training: {m} days ({dates[train][0].date()} to {dates[train][-1].date()})")
print(f"Test: {len(y_test)} days ({dates[~train][0].date()} to {dates[~train][-1].date()})")

Training: 456 days (2024-01-01 to 2025-03-31)
Test: 91 days (2025-04-01 to 2025-06-30)

Step 5: Gradient descent with Noah's update equations

---

Noah's notes give the update of each weight and of the bias for linear regression with the MSE cost:

$$w_j := w_j - \alpha \, \frac{2}{m} \sum_{i=1}^{m} (\hat{y}^{(i)} - y^{(i)})\, x_j^{(i)} \qquad b := b - \alpha \, \frac{2}{m} \sum_{i=1}^{m} (\hat{y}^{(i)} - y^{(i)})$$

With a matrix $X$ of $m$ rows and one column per variable, the five sums for the weights are a single product, $X^\top (\hat{y} - y)$. First we apply them to the variables as they come. Gradient descent diverges if the learning rate is too large, and with these variables the largest rate that is still stable can be computed: it depends on the biggest variable, the number of patients, which is in the hundreds.

def gradient_descent(X, y, learning_rate, iterations):
    m = len(y)
    w, b = np.zeros(X.shape[1]), 0.0
    for _ in range(iterations):
        error = X @ w + b - y                              # y_hat - y for every day
        w = w - learning_rate * (2 / m) * (X.T @ error)    # Noah's update, all the weights at once
        b = b - learning_rate * (2 / m) * error.sum()
    return w, b

def mse(X, y, w, b):
    return np.mean((X @ w + b - y) ** 2)

A = np.c_[X_train, np.ones(m)]
largest = np.linalg.eigvalsh((2 / m) * A.T @ A).max()
print(f"Largest stable learning rate with the raw variables: {2 / largest:.2e}")

w_raw, b_raw = gradient_descent(X_train, y_train, learning_rate=0.95 * 2 / largest, iterations=5000)
print("\nAfter 5000 iterations at 95% of that rate:")
print("  weights:", w_raw, " bias:", round(b_raw, 2))
print(f"  training MSE: {mse(X_train, y_train, w_raw, b_raw):.2f}")

Largest stable learning rate with the raw variables: 1.24e-05

After 5000 iterations at 95% of that rate:
  weights: [ 0.206 -0.432  3.34   0.201  0.376]  bias: 0.09
  training MSE: 76.17

The largest stable learning rate is only 1.24e-05, because the patients are in the hundreds while the alert and the weekend are 0 or 1. At that rate, after 5000 iterations, the model is still far from the data: the alert weight is 0.201 instead of about 18, the bias is 0.09 instead of about 44, and the training MSE is 76.17. Each step moves the weight of the patients a lot and the rest hardly at all, so the descent crawls along a long, narrow valley. It is the problem the from-scratch notebook of this chapter solves by scaling its single variable, now with five.

Step 6: Scaling the features

---

Noah's notes warn that gradient descent is sensitive to the scale of the variables, and recommend normalizing the data. We standardize each variable: subtract its mean and divide by its standard deviation, both computed on the training days only. Every variable then has mean 0 and standard deviation 1, the same learning rate suits all the weights, and the descent can take large steps. At the end we convert the weights and the bias back to the original units, so that each weight is again "masks per patient", "masks per degree" and so on:

$$w_j = \frac{w'_j}{\sigma_j} \qquad b = b' - \sum_j \frac{w'_j\,\mu_j}{\sigma_j}$$

mu, sigma = X_train.mean(axis=0), X_train.std(axis=0)
Z_train = (X_train - mu) / sigma

for iterations in [10, 50, 200, 1000]:
    w_s, b_s = gradient_descent(Z_train, y_train, learning_rate=0.1, iterations=iterations)
    print(f"{iterations:>5} iterations: training MSE {mse(Z_train, y_train, w_s, b_s):8.3f}")

w_gd = w_s / sigma                       # back to masks per patient, per degree, ...
b_gd = b_s - np.sum(w_s * mu / sigma)
print("\nWeights in original units:", w_gd, " bias:", round(b_gd, 2))

   10 iterations: training MSE  112.756
   50 iterations: training MSE   15.246
  200 iterations: training MSE   15.238
 1000 iterations: training MSE   15.238

Weights in original units: [ 0.044 -0.374  3.207 18.808  3.363]  bias: 44.3

With the variables standardized, the same function reaches a training MSE of 15.246 in 50 iterations and 15.238 from 200 on, close to the variance of the noise we added (16). Back in the original units, the model needs 0.044 masks per patient, −0.374 per degree, 3.207 per emergency, 18.808 more on an alert day and 3.363 more at weekends, over a base of 44.3.

Step 7: The exact solution: least squares and scikit-learn

---

For linear regression with the MSE cost there is also a direct solution, the least squares method that Martha mentions: the weights that minimize the squared error solve a system of linear equations, the normal equations $(A^\top A)\,\theta = A^\top y$, where $A$ is $X$ with an extra column of ones for the bias. NumPy solves it in one call, and scikit-learn's LinearRegression, which the team uses in the book, does the same inside fit:

theta = np.linalg.lstsq(A, y_train, rcond=None)[0]
model = LinearRegression().fit(X_train, y_train)

table = pd.DataFrame({'book equation': np.r_[book_w, book_b],
                      'gradient descent': np.r_[w_gd, b_gd],
                      'least squares': theta,
                      'scikit-learn': np.r_[model.coef_, model.intercept_]},
                     index=features + ['bias'])
print(table.round(3))

             book equation  gradient descent  least squares  scikit-learn
patients              0.05             0.044          0.044         0.044
temperature          -0.37            -0.374         -0.374        -0.374
emergencies           3.20             3.207          3.207         3.207
alert                18.70            18.808         18.808        18.808
weekend               4.10             3.363          3.363         3.363
bias                 42.30            44.304         44.304        44.304

The three methods agree to three decimals, and they are close to the book's equation: 0.044 against 0.05 per patient, −0.374 against −0.37 per degree, 3.207 against 3.2 per emergency, 18.808 against 18.7 for an alert and 44.304 against 42.3 for the bias. The largest difference is the weekend, 3.363 against 4.1: only one day in seven is a weekend, so its weight is estimated from fewer days and the noise shows more in it. With real data there is no equation to compare with, which is why the next step tests the model on days it has not seen.

Step 8: How good is the model?

---

We now use the three test months, which no method has seen. Besides the MSE and the MAE (the mean absolute error, in masks), we compute the coefficient of determination $R^2 = 1 - \frac{\sum (y - \hat{y})^2}{\sum (y - \bar{y})^2}$, the share of the day-to-day variation in masks that the model explains.

y_pred = model.predict(X_test)
print(f"Test MSE: {mean_squared_error(y_test, y_pred):.2f}   (the noise we added has variance 4² = 16)")
print(f"Test MAE: {mean_absolute_error(y_test, y_pred):.2f} masks")
print(f"Test R²:  {r2_score(y_test, y_pred):.3f}")

plt.figure(figsize=(12, 4))
plt.plot(dates[~train], y_test, label='masks used')
plt.plot(dates[~train], y_pred, label='prediction')
plt.title('The three test months')
plt.legend()
plt.show()

Test MSE: 13.34   (the noise we added has variance 4² = 16)
Test MAE: 2.75 masks
Test R²:  0.889

On the three test months the MSE is 13.34, the MAE 2.75 masks and R² 0.889: the model explains about 89% of the day-to-day variation, and what is left is mostly the noise we put in the data. The prediction follows the masks used closely.

Step 9: When the data has errors: MSE or MAE?

---

In the book, the team compares four cost functions and chooses the MSE, with a rule from Noah: if outliers appear, switch to the MAE. The squared error makes a large miss count much more (an error of 10 costs 100, as Martha points out), so a few wrong records can drag the whole model. We add five data-entry errors to the training days, masks typed with an extra zero, and fit the model twice: minimizing the MSE (LinearRegression) and minimizing the MAE (QuantileRegressor with the median, quantile=0.5, and no penalty).

y_typos = y_train.copy()
wrong_days = [40, 130, 222, 300, 390]
y_typos[wrong_days] = y_train[wrong_days] * 10              # e.g. 85 typed as 850
print("Five wrong records:", y_typos[wrong_days].round())

fits = {'MSE, clean data': LinearRegression().fit(X_train, y_train),
        'MSE, with the errors': LinearRegression().fit(X_train, y_typos),
        'MAE, with the errors': QuantileRegressor(quantile=0.5, alpha=0, solver='highs').fit(X_train, y_typos)}
rows = {name: np.r_[f.coef_, f.intercept_, mean_absolute_error(y_test, f.predict(X_test))] for name, f in fits.items()}
print(pd.DataFrame(rows, index=features + ['bias', 'test MAE']).round(3))

Five wrong records: [1146.  934.  966.  813. 1106.]
             MSE, clean data  MSE, with the errors  MAE, with the errors
patients               0.044                 0.220                 0.039
temperature           -0.374                -0.504                -0.390
emergencies            3.207                 4.601                 3.167
alert                 18.808                 4.373                18.741
weekend                3.363                36.153                 3.234
bias                  44.304               -17.253                46.452
test MAE               2.750                11.607                 2.755

Five wrong records out of 456 are enough to ruin the MSE fit: the alert weight falls from 18.808 to 4.373, the weekend weight jumps to 36.153 and the test MAE rises from 2.75 to 11.607 masks. The MAE fit, trained on the same wrong data, stays almost where the clean fit was (alert 18.741, test MAE 2.755). That is Noah's rule in numbers. Even so, the first remedy is to find and correct the wrong records: the MAE limits the damage, it does not tell us where the errors are.

Step 10: The forecast for next week

---

Finally, the week the team uses in the book to test its model: Monday with 295 patients, 14 °C and 12 expected emergencies; Tuesday with 310 patients, 15 °C and 10 emergencies; Wednesday with 285 patients, 16 °C and 11 emergencies; no health alert and no weekend.

week = pd.DataFrame([[295, 14, 12, 0, 0], [310, 15, 10, 0, 0], [285, 16, 11, 0, 0]],
                    index=['Monday', 'Tuesday', 'Wednesday'], columns=features)
week['book equation'] = week[features].values @ book_w + book_b
week['our model'] = model.predict(week[features].values)
print(week.round(2))

           patients  temperature  emergencies  alert  weekend  book equation  \
Monday          295           14           12      0        0          90.27   
Tuesday         310           15           10      0        0          84.25   
Wednesday       285           16           11      0        0          85.83   

           our model  
Monday         90.63  
Tuesday        84.51  
Wednesday      86.23  

Our model predicts 90.63, 84.51 and 86.23 masks for Monday, Tuesday and Wednesday. The book's equation gives exactly 90.27, 84.25 and 85.83, which the book rounds up to 91, 85 and 86. Both are inside Liam's historical ranges: 85 to 95 masks for a Monday like this one, 82 to 86 for the Tuesday and 80 to 90 for the Wednesday.

Conclusion

---

• With several variables, linear regression is still a weighted sum, $\hat{y} = \mathbf{w} \cdot \mathbf{x} + b$: a hyperplane instead of a line, and each weight reads as the effect of one variable when the others do not change.
• Noah's update equations work unchanged for any number of variables, but only if the variables are on similar scales. Standardizing them turns thousands of useless iterations into a few hundred useful ones.
• For linear regression, least squares gives the exact answer directly, and scikit-learn uses it. Gradient descent matters because it is the method that also works for the models of the following chapters.
• The cost function is a decision: the MSE is the usual choice, but a few wrong records can ruin it, and the MAE resists them.