Linear Regression

Linear regression is the model every data scientist should understand deeply before reaching for anything more complex — not because it is simple, but because every failure mode it has (violated assumptions, multicollinearity, nonlinearity) reappears in subtler forms in every model that follows. The practitioner who can diagnose a linear regression confidently is already thinking in the right way about all regression.

Learning Objectives

• Write the linear regression equation and explain what each term means geometrically
• Understand the intuition behind Ordinary Least Squares without needing the full matrix derivation
• State the four key assumptions and know how to check each one with plots
• Fit a model with scikit-learn and interpret its coefficients correctly
• Recognise when a linear model is the right choice and when it is not

---

The Equation and What It Means

The linear regression model says:

y = β₀ + β₁x₁ + β₂x₂ + ... + βₙxₙ + ε

Each piece:

Symbol	Name	Meaning
y	Target	The value you are predicting
β₀	Intercept	The predicted value when all features equal zero
β₁...βₙ	Coefficients	How much y changes per one-unit increase in each feature, holding all others constant
xᵢ	Feature values	The observed input data
ε	Error term	The irreducible noise — variation in y not explained by the features

The phrase "holding all others constant" is the most important and most misunderstood part of reading coefficients. When you see that bedrooms has a coefficient of 15,000, it does not mean adding a bedroom unconditionally adds $15,000 — it means adding a bedroom while keeping square footage, location, and age fixed adds $15,000. This distinction matters enormously when features are correlated.

---

The Geometry: Why Least Squares

Linear regression draws the hyperplane that minimises the sum of squared vertical distances from each observed point to the plane. These vertical distances are the residuals — the gaps between what the model predicts and what actually happened.

Imagine you have house prices and square footage plotted as a scatter. Dozens of lines could pass through that cloud. Linear regression finds the unique line where the sum of all the squared vertical gaps is smaller than for any other line. Every other possible line you could draw would have a higher total squared error. That is the entire optimisation problem.

Info

Why square the errors instead of taking absolute values? Squared errors are mathematically convenient — the squared loss is smooth and differentiable everywhere, which means gradient-based solvers work cleanly. Absolute errors create a kink at zero that complicates optimisation. The practical consequence is that MSE-optimised models are more sensitive to outliers than MAE-optimised ones. That sensitivity is a bug in some situations and a feature in others.

The OLS Solution

For linear regression specifically, the minimum is reachable in closed form. The optimal coefficients are:

β = (XᵀX)⁻¹ Xᵀy

You do not need to memorise this. What matters is what it implies:

1. If XᵀX is not invertible (perfect multicollinearity), the system has no unique solution — Ridge regression fixes this by adding a small constant to the diagonal before inverting
2. The solution is computed exactly, in one step — no iterative gradient descent needed
3. Adding more features always reduces training error (because the model has more degrees of freedom) — but this is why we evaluate on held-out test data

---

Fitting with scikit-learn

import pandas as pd
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error, r2_score

# Load a realistic regression dataset
housing = fetch_california_housing(as_frame=True)
X = housing.data
y = housing.target  # median house value in units of $100,000

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Pipeline: scale then fit (good habit, required for regularised models)
pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("model", LinearRegression())
])

pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)

mae  = mean_absolute_error(y_test, y_pred)
r2   = r2_score(y_test, y_pred)

print(f"MAE:  {mae:.3f}  (units of $100k → ~${mae * 100_000:,.0f} average error)")
print(f"R²:   {r2:.3f}")
# Output:
# MAE:  0.533  (units of $100k → ~$53,300 average error)
# R²:   0.576

An R² of 0.576 means the model explains about 58% of the variance in house prices. That is a reasonable starting point for a raw linear model on this dataset — not good enough for a production pricing tool, but a meaningful baseline.

---

Interpreting Coefficients

# Extract coefficients from the pipeline
lr_model = pipe.named_steps["model"]
scaler   = pipe.named_steps["scaler"]

coef_series = pd.Series(lr_model.coef_, index=X.columns).sort_values()
print(coef_series)
# Output (approximate, scaled units):
# AveOccup    -0.196
# Latitude    -0.899
# ...
# MedInc       0.855

print(f"Intercept: {lr_model.intercept_:.3f}")

Warning

When you use StandardScaler, coefficients are in standardised units, not original feature units. The coefficient for MedInc (median income) of 0.855 means: a one-standard-deviation increase in median income predicts a 0.855-unit increase in the median house value (in $100k units). To compare feature importance, standardised coefficients work. To communicate "a $10,000 raise in income predicts..." you need to back-transform using the scaler's mean_ and scale_ attributes.

Coefficients with raw (unscaled) features:

# Refit without scaling just to show raw coefficient interpretation
raw_model = LinearRegression().fit(X_train, y_train)
raw_coefs = pd.Series(raw_model.coef_, index=X.columns).sort_values()
print(raw_coefs)
# Output (approximate):
# AveOccup   -0.004
# Latitude   -0.419
# ...
# MedInc      0.436

MedInc coefficient of 0.436: for every one-unit increase in median income (recall the dataset's income unit is tens of thousands of dollars), predicted house value increases by $43,600, holding other features constant.

---

The Four Assumptions — and How to Check Them

Linear regression has four key assumptions. Violating them does not necessarily mean the model is wrong — it means the coefficient estimates and confidence intervals may be unreliable, and interpretation requires care.

1. Linearity

The relationship between features and target is (approximately) linear.

How to check: Plot each feature against the residuals. A random scatter means linearity holds. A curved pattern (like a parabola) means the relationship is nonlinear and linear regression is missing it.

import matplotlib.pyplot as plt

residuals = y_test - y_pred

# Residuals vs. a specific feature
plt.figure(figsize=(8, 4))
plt.scatter(X_test["MedInc"], residuals, alpha=0.3, s=10)
plt.axhline(0, color="red", linewidth=1)
plt.xlabel("Median Income")
plt.ylabel("Residual")
plt.title("Residuals vs. MedInc — check for curvature")
plt.tight_layout()
plt.show()

2. Independence

Observations are independent of each other. This assumption is mostly about data collection — if you have time series data or repeated measurements on the same subjects, residuals are likely correlated.

How to check: Plot residuals in time order. Patterns (trending up or down, oscillating) suggest autocorrelation. For time series, use the Durbin-Watson test.

3. Homoscedasticity (Constant Variance)

The spread of residuals is the same across all predicted values. When residuals fan out as predictions grow, you have heteroscedasticity — a sign that the model is less reliable for high-value predictions.

How to check: Plot residuals vs. predicted values (the most important diagnostic plot).

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Residuals vs. Fitted
axes[0].scatter(y_pred, residuals, alpha=0.3, s=10)
axes[0].axhline(0, color="red", linewidth=1)
axes[0].set_xlabel("Predicted Values")
axes[0].set_ylabel("Residuals")
axes[0].set_title("Residuals vs. Fitted")

# Distribution of residuals
axes[1].hist(residuals, bins=50, edgecolor="black")
axes[1].set_xlabel("Residual")
axes[1].set_title("Residual Distribution")

plt.tight_layout()
plt.show()

A "good" residuals-vs-fitted plot looks like a random cloud centred on zero with consistent vertical spread across the horizontal axis. If the cloud fans outward (heteroscedastic), consider log-transforming the target.

4. Normality of Residuals

Residuals are approximately normally distributed. This matters most for statistical inference (p-values, confidence intervals) rather than for prediction accuracy.

How to check: Q-Q plot of residuals.

import scipy.stats as stats

stats.probplot(residuals, dist="norm", plot=plt)
plt.title("Q-Q Plot of Residuals")
plt.show()

Points following the diagonal line closely indicate normality. Heavy tails or S-curves suggest skewed residuals.

Tip

In practice, the most important diagnostic is the residuals vs. fitted plot. If you only look at one plot, make it that one. It simultaneously reveals non-linearity, heteroscedasticity, and outliers. A pattern in this plot tells you the model is systematically wrong in some region of the prediction space.

---

The Full Diagnostic Workflow

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error, r2_score
import scipy.stats as stats

housing = fetch_california_housing(as_frame=True)
X, y = housing.data, housing.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

pipe = Pipeline([("scaler", StandardScaler()), ("model", LinearRegression())])
pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)
residuals = y_test.values - y_pred

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 1. Residuals vs. Fitted
axes[0].scatter(y_pred, residuals, alpha=0.2, s=8)
axes[0].axhline(0, color="red", linewidth=1)
axes[0].set_xlabel("Predicted")
axes[0].set_ylabel("Residual")
axes[0].set_title("Residuals vs. Fitted")

# 2. Residual histogram
axes[1].hist(residuals, bins=60, edgecolor="black")
axes[1].set_xlabel("Residual")
axes[1].set_title("Residual Distribution")

# 3. Q-Q plot
stats.probplot(residuals, dist="norm", plot=axes[2])
axes[2].set_title("Q-Q Plot")

plt.tight_layout()
plt.show()

# Summary metrics
print(f"MAE:  {mean_absolute_error(y_test, y_pred):.3f}")
print(f"R²:   {r2_score(y_test, y_pred):.3f}")
print(f"Residual std: {np.std(residuals):.3f}")
print(f"Residual mean: {np.mean(residuals):.4f}  (should be near 0)")
# Output:
# MAE:  0.533
# R²:   0.576
# Residual std: 0.745
# Residual mean: 0.0000  (should be near 0)

---

When Linear Regression Is the Right Choice

Use linear regression when:

• You need an interpretable model — stakeholders will ask "why did you predict $450,000?"
• The relationship is plausibly linear (or can be made linear with log transforms)
• You have limited data and a complex model would overfit
• You are building a baseline — every project should have one

Warning

The R² trap. An R² of 0.95 does not mean your model is correct. It means your model explains 95% of the variance in the training data. If the residuals vs. fitted plot shows a clear U-shape, your model is systematically wrong in a way R² completely hides. Always check residuals. Always.

Use a different model when:

• Relationships are highly nonlinear (tree-based models)
• You have many features with little data (Ridge or Lasso)
• You care about ranking rather than precise values (consider rank-based metrics)
• The target is a count, probability, or duration (consider GLMs)

---

Common Coefficient Interpretation Mistakes

Warning

Multicollinearity makes individual coefficients unreliable. If square_footage and total_rooms are highly correlated, the model may assign a large positive coefficient to one and a large negative coefficient to the other — not because either truly drives price in that direction, but because the model is splitting the shared variance between them arbitrarily. Check the Variance Inflation Factor (VIF) for any features with suspicious coefficients.

from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_data = pd.DataFrame({
    "Feature": X_train.columns,
    "VIF": [
        variance_inflation_factor(X_train.values, i)
        for i in range(X_train.shape[1])
    ]
}).sort_values("VIF", ascending=False)

print(vif_data)
# VIF > 10 signals problematic multicollinearity
# VIF > 5 is worth investigating

Success

Linear regression's greatest strength is its transparency. When a model says a feature has a coefficient of -0.42, you can ask "does it make sense for this variable to decrease the prediction?" and answer that question without opening a black box. This interpretability makes linear regression the default choice in regulated industries (credit, healthcare) where predictions must be explainable.

---

---

What's Next

You've covered OLS linear regression, coefficient interpretation, the four regression assumptions (linearity, independence, homoscedasticity, normality of residuals), multicollinearity detection with VIF, the residuals-vs-fitted diagnostic workflow, and when to use a different model. Next up: 03-ridge-lasso-elasticnet — where you'll learn how L1 and L2 regularisation penalties shrink and zero-out coefficients, how to use RidgeCV and LassoCV to automatically find the optimal alpha, and when ElasticNet outperforms either alone.

Optional Deep Dive

Read "ISLR" Chapter 6 (Shrinkage Methods) — it provides the geometric and algebraic intuition for why Ridge shrinks all coefficients toward zero while Lasso zeros them out exactly, using the constraint-based view of regularisation that makes the difference immediately clear.

00-agenda | 03-ridge-lasso-elasticnet