Multiple Linear Regression

Simple linear regression finds the best line through one feature. Multiple linear regression extends that to $p$ features simultaneously — and the key word is simultaneously. Each coefficient now measures the isolated effect of its feature while all other features are held constant. That "holding constant" is what makes multiple regression harder to interpret and more powerful than running $p$ separate simple regressions.

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.

OLS Solution for Multiple Features

${\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}")

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.

Predicting with sklearn

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}")

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

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

Multicollinearity Warning

If sqft and bedrooms are highly correlated, $X^{\top }X$ becomes nearly singular — the inversion is numerically unstable.

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

Correlation of 0.997 is nearly perfect multicollinearity. 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).

Simple vs Multiple Regression Comparison

Related Concepts and Honest Limitations

Multiple regression's "partial effect" interpretation is only valid when the other features are actually held constant — which requires them to be linearly independent. Near-multicollinearity breaks this assumption without making the model fail to fit. The model might achieve low training error while the individual coefficients are statistically meaningless. 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. Adjusted $R^2$ (next post) corrects for this by penalizing model complexity.

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?