Multiple Linear Regression

You have a house's square footage — that's one feature. But you also know the number of bedrooms, the lot size, the year built. Running separate simple regressions for each gives conflicting pictures: the sqft coefficient is 0.200 when sqft is the only predictor, but what happens when you add bedrooms? Houses with more bedrooms tend to be larger, so the sqft coefficient was partly capturing bedroom effects. You need to isolate each variable's true contribution while holding all others constant.

That "holding constant" — ceteris paribus — is what multiple regression does, and it's also why it's harder to interpret than running $p$ separate lines. This is not the same as fitting $p$ independent simple regressions; the coefficients change when you add or remove features because each one now represents a partial effect.

The Plan — Four Steps to Multiple Regression

We'll extend the simple line to multiple features, trace a prediction through the design matrix, find the optimal weights with OLS, then examine what happens when features are correlated.

Step 1 — From Simple to Multiple

Simple: $\hat{y}=w_0+w_1\times \text{sqft}$

Multiple: $\hat{y}=w_0+w_1\times \text{sqft}+w_2\times \text{bedrooms}$

Adding bedrooms introduces $w_2$, the partial effect of bedroom count while holding sqft fixed. This "ceteris paribus" interpretation — Latin for "all other things equal" — is what each coefficient in multiple regression represents. It's also why the coefficients change when you add or remove features.

Anchor dataset:

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

X = np.array([
    [650,  2],
    [850,  2],
    [1100, 3],
    [1400, 3],
    [1600, 4],
    [1900, 4]
])
y = np.array([180, 220, 280, 340, 370, 430])

The Design Matrix $X$

Append a column of ones for the intercept. The augmented matrix $X_{\text{aug}}$ and weight vector $\mathbf{w}$ are:

$X_{\text{aug}}=\left[\begin{array}{ccc}1& 650& 2\\ 1& 850& 2\\ 1& 1100& 3\\ 1& 1400& 3\\ 1& 1600& 4\\ 1& 1900& 4\end{array}\right],\mathbf{w}=\left[\begin{array}{c}{w}_0\\ w_1\\ w_2\end{array}\right]$

Predictions for all samples at once: $\hat{\mathbf{y}}=X_{\text{aug}}\mathbf{w}$.

Manual prediction trace with illustrative weights $w_0=30$, $w_1=0.17$, $w_2=15$:

SSE = $9.5^2+15.5^2+18^2+27^2+8^2+17^2=90.25+240.25+324+729+64+289=1736.5$

MSE = $1736.5/6=289.4$ — much worse than the simple regression MSE of 22.2 because these weights aren't optimal.

✓ Step 1 complete. Extending from one to two features changes the model from a line to a plane. The illustrative weights ($w_0=30,w_1=0.17,w_2=15$) are not optimal — SSE = 1,736.5.

Step 2 — The Design Matrix $X$

Predicting all 6 samples at once requires matrix operations. We augment $X$ with a column of ones for the intercept, then predictions are a single matrix multiplication:

$X_{\text{aug}}=\left[\begin{array}{ccc}1& 650& 2\\ 1& 850& 2\\ 1& 1100& 3\\ 1& 1400& 3\\ 1& 1600& 4\\ 1& 1900& 4\end{array}\right],\hat{\mathbf{y}}=X_{\text{aug}}\mathbf{w}$

✓ Step 2 complete. The design matrix bundles all feature columns plus an intercept column. Predictions become a single (6×3) @ (3×1) matrix multiply.

Step 3 — OLS Solution and Coefficient Interpretation

The optimal weights minimize SSE across all samples simultaneously:

${\mathbf{w}}^{\ast }=(X^{\top }X{)}^{-1}{X}^{\top }\mathbf{y}$

In practice, use np.linalg.solve instead of inverting directly:

X_aug = np.column_stack([np.ones(6), X])
w_ols = np.linalg.solve(X_aug.T @ X_aug, X_aug.T @ y)
print(f"w₀ = {w_ols[0]:.4f}")
print(f"w₁ = {w_ols[1]:.6f}")
print(f"w₂ = {w_ols[2]:.4f}")

w₀ = 58.3271
w₁ = 0.175000
w₂ = 9.6053

Coefficient Interpretation

• $w_0=58.33$: baseline price when sqft = 0 and bedrooms = 0. Mathematically required; not a meaningful real-world prediction (no house has zero sqft).
• $w_1=0.175$: holding bedrooms fixed, each additional square foot adds $175 to the predicted price.
• $w_2=9.61$: holding sqft fixed, each additional bedroom adds $9,610 to the predicted price.

Critical observation: $w_1$ changed from 0.200 (simple regression) to 0.175 (multiple regression). This is omitted variable bias. Larger houses tend to have more bedrooms — in simple regression, the sqft coefficient was absorbing bedroom effects. Adding bedrooms to the model isolated each variable's true independent contribution.

Within each bedroom group, the sqft slope is shallower (≈0.175) than the combined slope (0.200). The combined slope was inflated because larger sqft correlates with more bedrooms.

✓ Step 3 complete. OLS gives $w_0=58.33$, $w_1=0.175$, $w_2=9.61$. The sqft coefficient dropped from 0.200 (simple) to 0.175 (multiple) — omitted variable bias exposed.

Step 4 — sklearn Verification and Multicollinearity Warning

model = LinearRegression()

```python
model = LinearRegression()
model.fit(X, y)

print(f"w₀: {model.intercept_:.2f}")
print(f"w₁ (sqft):     {model.coef_[0]:.4f}")
print(f"w₂ (bedrooms): {model.coef_[1]:.4f}")

w₀: 58.30
w₁ (sqft):     0.1750
w₂ (bedrooms): 9.6000

new_house = np.array([[1200, 3]])
print(f"Predicted price: ${model.predict(new_house)[0]:.1f}k")

Predicted price: $297.3k

Manual check: $58.3+0.175\times 1200+9.6\times 3=58.3+210+28.8=297.1\text{k}$ ✓

Now check the correlation between our two features:

np.corrcoef(X[:, 0], X[:, 1])[0, 1]

0.9972

Correlation of 0.997 is nearly perfect multicollinearity. When $X^{\top }X$ is nearly singular, the inversion becomes numerically unstable. On a real dataset this would cause coefficient instability — tiny changes in the training data could swing $w_1$ and $w_2$ dramatically while leaving predictions nearly the same. This is the signal to apply Ridge regularization (post 14).

✓ Step 4 complete. sklearn matches the manual OLS. High feature correlation (0.997) threatens coefficient stability — a preview of regularization.

Related Concepts and Honest Limitations

A frustrating discovery with multiple regression is that the "partial effect" interpretation breaks silently under multicollinearity. The model still fits well (low training error), but the individual coefficients become statistically meaningless — they can even flip signs under resampling. Always check VIF (Variance Inflation Factor) when features are correlated.

Adding features always reduces or holds constant the training SSE — it can never increase it. This means you can always improve training error by adding noise features, which is why adjusted $R^2$ (next post) penalizes complexity. The OLS formula $w^{\ast }=(X^{\top }X{)}^{-1}{X}^{\top }y$ looks clean on paper, but when $p$ approaches $n$ (more features than samples), $X^{\top }X$ becomes singular and the solution doesn't exist without regularization.

Test Your Understanding

1. The simple regression coefficient for sqft was 0.200 and dropped to 0.175 in multiple regression. Intuitively, what happens to this coefficient if you add a third feature that is completely uncorrelated with sqft?

2. Compute the prediction for a house with sqft = 750 and 2 bedrooms using the optimal OLS weights ($w_0=58.3$, $w_1=0.175$, $w_2=9.6$). Compare to the simple regression prediction for sqft = 750. Which is higher and why?

3. If bedrooms and sqft had zero correlation in our dataset, would the multiple regression $w_1$ still differ from the simple regression $w_1$? Why?

4. In the design matrix $X_{\text{aug}}$, what happens if two columns are identical? What does $X^{\top }X$ look like and why does $(X^{\top }X{)}^{-1}$ fail?

5. The OLS solution minimizes SSE over all possible $\mathbf{w}$. Could you find weights with lower SSE on the training set by adding a random noise feature? Would that generalize to test data?