M44: OLS geometry & assumptions

CSCI 1109 — Practical Data Science

Frank Rudzicz

Learning outcomes

By the end of this module, you should be able to:

  • Explain linear regression as projecting a target vector onto a feature space.
  • Write the matrix form \(y = X\beta + \varepsilon\) and interpret \(\beta, \varepsilon\).
  • Derive (lightly) the normal equations and see why residuals are orthogonal.
  • Diagnose basic violations of OLS assumptions using residual plots and context.
  • Connect these assumptions to honest evaluation and responsible decisions in higher-impact settings.

Warm-up: regression in words

We want a model that predicts a numeric outcome \(y\):

  • Examples:
    • Disease progression score from clinical measurements
    • Daily air pollution level from traffic + weather
    • Vehicle fuel efficiency from weight + horsepower

Linear regression says:

“Predict \(y\) as a linear combination of features \(x_1, x_2, \dots, x_p\).”

\[ \hat y = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p \]

We choose coefficients \(\beta\) to make predictions “close” to the observed \(y\) on the training data.

What does “close” mean? → Least squares: minimize the sum of squared errors.

Matrix form: \(y = X\beta + \varepsilon\)

For \(n\) data points, with \(p\) features:

  • \(y \in \mathbb{R}^n\) – vector of outcomes
  • \(X \in \mathbb{R}^{n \times (p+1)}\) – design matrix (input features)

\[ X = \begin{bmatrix} 1 & x_{11} & x_{12} & \dots & x_{1p} \\ 1 & x_{21} & x_{22} & \dots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \dots & x_{np} \end{bmatrix} \]

  • First column of 1’s = intercept
  • \(\beta \in \mathbb{R}^{p+1}\) – coefficients

Model of reality:

\[ y = X\beta + \varepsilon \]

where:

  • \(X\beta\) – model’s predictions (a vector in \(\mathbb{R}^n\))
  • \(\varepsilon\) – residuals/noise (also in \(\mathbb{R}^n\))

OLS objective

Ordinary least squares chooses \(\hat\beta\) to minimize squared residuals:

\[ \hat\beta = \arg\min_{\beta} \, \|y - X\beta\|^2 \]

  • \(\| \cdot \|\) is the usual Euclidean norm
  • \(\|y - X\beta\|^2 = (y - X\beta)^\top(y - X\beta)\)

OLS gives one best-fitting vector of predictions \(X\hat\beta\) among all linear combinations of the columns of \(X\).

Key idea: the columns of \(X\) span a subspace of \(\mathbb{R}^n\).
OLS finds the closest point in that subspace to the actual data vector \(y\).

🔬 Geometry and projection

Think of \(\mathbb{R}^n\) as one big space of length-\(n\) vectors (for \(n\) data points):

  • Our target values (one number per data point) form a single vector \(y\) in this space.
  • Our features give us columns of \(X\) — each column is also a length-\(n\) vector.

Now look at all the vectors you can get by taking linear combinations of the columns of \(X\) (things of the form \(X\beta\)).

Those all live together in a flat “sheet” inside \(\mathbb{R}^n\) called the column space of \(X\):

  • When \(p = 1\), that sheet is just a line.
  • When \(p = 2\), it’s a plane.
  • In general (if features are independent), it’s a \((p+1)\)-dimensional flat surface.

OLS says: “Among all those \(X\beta\) on that sheet, pick the one closest to \(y\) in straight-line distance.”

  • \(X\hat\beta\) is that closest point — the projection of \(y\) onto the column space of \(X\).
  • The residual vector \(r = y - X\hat\beta\) is the “gap” from the sheet up to \(y\), and it ends up 📐 perpendicular to that sheet.

Regression as projection

Note

We have lots of possible straight-line prediction vectors (the plane), and one true outcome vector \(y\). OLS just picks the prediction vector on that plane that’s closest to \(y\), like dropping \(y\) straight down onto the plane and using its shadow.

OLS as orthogonal projection: projecting y onto the column space of X.

🔬 Light derivation: normal equations

We want to choose \(\beta\) to minimize:

\[ L(\beta) = \|y - X\beta\|^2 = (y - X\beta)^\top (y - X\beta) \]

Code 
import numpy as np
import matplotlib.pyplot as plt

beta_vals = np.linspace(-3, 3, 200)
# toy 1D quadratic loss: L(β) = (β - 1)^2 + 2
loss = (beta_vals - 1)**2 + 2

plt.plot(beta_vals, loss)
plt.axvline(1, linestyle="--")          # true minimizer β* = 1
plt.scatter([1], [2], zorder=3)         # point at the bottom
plt.xlabel("β")
plt.ylabel("L(β)")
plt.title("Quadratic 'bowl': minimum where gradient = 0")
plt.tight_layout()
plt.show()

A simple quadratic loss: minimum where the slope (gradient) is zero.

Because \(L(\beta)\) is a smooth quadratic bowl (a convex function of \(\beta\)), its unique minimum is exactly where the gradient is zero, so solving \(\nabla_\beta L(\beta) = 0\) gives us the best \(\hat\beta\).

  1. Expand \(L(\beta) = \|y - X\beta\|^2\) and differentiate:

\[ \nabla_\beta L = -2X^\top(y - X\beta) \]

  1. Set \(\nabla_\beta L = 0\):

\[ X^\top y - X^\top X \beta = 0 \]

  1. Rearrange:

\[ X^\top X \hat\beta = X^\top y \quad \text{(normal equations)} \]

If \(X^\top X\) is invertible:

\[ \boxed{\hat\beta = (X^\top X)^{-1} X^\top y} \]

Orthogonality condition

From the normal equations:

\[ X^\top (y - X\hat\beta) = 0 \]

Let \(r = y - X\hat\beta\) (residual vector). Then:

\[ X^\top r = 0 \]

Meaning:

  • Each column of \(X\) is orthogonal to the residual vector.

  • In other words:

    Once you’ve chosen the OLS line/plane, there’s no remaining linear pattern between residuals and any feature in \(X\).

Worked example 1

Diabetes progression (1D)

Dataset: classic clinical dataset on diabetes progression.

  • Each row = one patient
  • Features = baseline measurements (BMI, blood pressure, etc.)
  • Target = disease progression score one year later

We’ll start with a simple 1D regression:

Predict disease progression from BMI alone.

Code 
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression

diab = load_diabetes(as_frame=True)
df_d = diab.frame.copy()
df_d = df_d.rename(columns={"target": "progression"})

df_d[["bmi", "progression"]].head()
bmi progression
0 0.061696 151.0
1 -0.051474 75.0
2 0.044451 141.0
3 -0.011595 206.0
4 -0.036385 135.0

Visualizing the relationship (scatter)

Code 
plt.scatter(df_d["bmi"], df_d["progression"])
plt.xlabel("BMI (standardized)")
plt.ylabel("Disease progression (1-year)")
plt.title("Diabetes progression vs BMI")
plt.tight_layout()
plt.show()

✏️ Questions

  • Does a straight line look at least plausible?
  • Any obvious non-linear pattern?

Fit the OLS line (progression ~ bmi)

Code 
X = df_d[["bmi"]]
y = df_d["progression"]

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

df_d["pred"] = model.predict(X)

print("Intercept (beta_0):", model.intercept_)
print("Slope (beta_1):    ", model.coef_[0])
df_d.head()
Intercept (beta_0): 152.13348416289617
Slope (beta_1):     949.4352603840372
age sex bmi bp s1 s2 s3 s4 s5 s6 progression pred
0 0.038076 0.050680 0.061696 0.021872 -0.044223 -0.034821 -0.043401 -0.002592 0.019907 -0.017646 151.0 210.710038
1 -0.001882 -0.044642 -0.051474 -0.026328 -0.008449 -0.019163 0.074412 -0.039493 -0.068332 -0.092204 75.0 103.262195
2 0.085299 0.050680 0.044451 -0.005670 -0.045599 -0.034194 -0.032356 -0.002592 0.002861 -0.025930 141.0 194.337033
3 -0.089063 -0.044642 -0.011595 -0.036656 0.012191 0.024991 -0.036038 0.034309 0.022688 -0.009362 206.0 141.124769
4 0.005383 -0.044642 -0.036385 0.021872 0.003935 0.015596 0.008142 -0.002592 -0.031988 -0.046641 135.0 117.588574

Visualizing the fit (line + data)

Code 
plt.scatter(df_d["bmi"], df_d["progression"], label="Observed")
# Sort for a clean line
order = np.argsort(df_d["bmi"].values)
plt.plot(df_d["bmi"].values[order], df_d["pred"].values[order], label="OLS fit")
plt.xlabel("BMI (standardized)")
plt.ylabel("Disease progression (1-year)")
plt.title("Diabetes progression vs BMI (OLS fit)")
plt.legend()
plt.tight_layout()
plt.show()

Interpretation

  • Slope: how much progression score changes per unit increase in standardized BMI
  • Intercept: predicted progression when BMI = 0 (average baseline)

Residuals vs BMI (check orthogonality)

Code 
df_d["resid"] = df_d["progression"] - df_d["pred"]

plt.scatter(df_d["bmi"], df_d["resid"])
plt.axhline(0, linestyle="--")
plt.xlabel("BMI (standardized)")
plt.ylabel("Residual (progression - prediction)")
plt.title("Residuals vs BMI")
plt.tight_layout()
plt.show()

Residuals vs BMI — look for patterns.

👀 Visual check:

  • Residuals should be:
    • roughly centered around 0 across BMI values
    • no obvious U-shape / non-linear curve
  • This is the visual version of “no remaining linear pattern with BMI”.

Residual distribution (Normal?)

Code 
plt.hist(df_d["resid"], bins=20)
plt.xlabel("Residual")
plt.ylabel("Count")
plt.title("Residual distribution: progression ~ BMI")
plt.tight_layout()
plt.show()

Histogram of residuals for diabetes progression model.

We’re looking for:

  • Rough symmetry around 0
  • No giant spikes or extreme tails

This matters more when we want confidence intervals / p-values, less for pure prediction.

OLS assumptions (informal)

In practice, OLS is used when these are at least roughly true:

  1. Linearity & additivity
    • The true relationship is approximately linear in the features.
  2. Independent errors
    • Residuals aren’t strongly dependent on each other (no serious autocorrelation).
  3. Homoskedasticity (constant variance)
    • Residual variance is roughly constant across fitted values / key features.
  4. Approximate normality of residuals
    • For inference (confidence intervals, p-values).
  5. No extreme multicollinearity
    • Features aren’t nearly perfect linear combos of each other.

Plus earlier good practice:

  • No major leakage.
  • Features and outcome measured sensibly (units, scales, etc.).

Assumption checks: what we can see

Simple visual diagnostics:

  • Linearity:
    • Residuals vs fitted values: should look like random scatter around 0.
  • Homoskedasticity:
    • Same plot: roughly equal vertical spread across the x-axis.
  • Outliers / leverage:
    • Look for points far from the cloud or with big residuals.
  • Normality (if we care about inference):
    • Histogram or Q–Q plot of residuals.

We saw a mostly okay example (progression ~ BMI).
Next: a higher-impact anti-example from energy & climate.

Worked example 2

Fuel efficiency & climate impact

Dataset: car fuel efficiency dataset (mpg).

  • Each row = one car model
  • Features = horsepower, weight, cylinders, etc.
  • Target = fuel efficiency (miles per gallon, mpg)

Better mpg → lower fuel consumption → lower emissions.

Code 
import seaborn as sns

mpg = sns.load_dataset("mpg").dropna()
mpg[["mpg", "weight", "horsepower"]].head()
mpg weight horsepower
0 18.0 3504 130.0
1 15.0 3693 165.0
2 18.0 3436 150.0
3 16.0 3433 150.0
4 17.0 3449 140.0

Visualize mpg vs weight

Code 
plt.scatter(mpg["weight"], mpg["mpg"], alpha=0.6)
plt.xlabel("Vehicle weight")
plt.ylabel("Fuel efficiency (mpg)")
plt.title("mpg vs weight")
plt.tight_layout()
plt.show()

✏️ Questions

  • Is the relationship roughly linear?
  • Does the spread change as weight increases?

Fit OLS: mpg ~ weight

Code 
from sklearn.linear_model import LinearRegression

X_m = mpg[["weight"]]
y_m = mpg["mpg"]

lin_m = LinearRegression()
lin_m.fit(X_m, y_m)

mpg["pred_mpg"] = lin_m.predict(X_m)

print("Intercept:", lin_m.intercept_)
print("Slope   :", lin_m.coef_[0])
Intercept: 46.216524549017564
Slope   : -0.007647342535779576

Plot the fit: line + data

Code 
plt.scatter(mpg["weight"], mpg["mpg"], alpha=0.6, label="Observed")
order = np.argsort(mpg["weight"].values)
plt.plot(mpg["weight"].values[order], mpg["pred_mpg"].values[order], label="OLS fit")
plt.xlabel("Vehicle weight")
plt.ylabel("Fuel efficiency (mpg)")
plt.title("mpg vs weight (OLS fit)")
plt.legend()
plt.tight_layout()
plt.show()

This looks somewhat reasonable… but is it good enough?

Residuals vs fitted mpg: looking for trouble

Code 
mpg["resid_mpg"] = mpg["mpg"] - mpg["pred_mpg"]

plt.scatter(mpg["pred_mpg"], mpg["resid_mpg"], alpha=0.6)
plt.axhline(0, linestyle="--")
plt.xlabel("Fitted mpg")
plt.ylabel("Residual")
plt.title("Residuals vs fitted mpg")
plt.tight_layout()
plt.show()

Residuals vs fitted mpg — check linearity & variance.

Watch for:

  • Curves / systematic bends → violated linearity
  • “Fanning out” (bigger spread at some fitted values) → violated homoskedasticity
  • Individual huge residuals → potential outliers

Residuals vs weight: alternative view

Code 
plt.scatter(mpg["weight"], mpg["resid_mpg"], alpha=0.6)
plt.axhline(0, linestyle="--")
plt.xlabel("Vehicle weight")
plt.ylabel("Residual")
plt.title("Residuals vs weight")
plt.tight_layout()
plt.show()

Warning

If we see:

  • Higher variance for heavier cars, or
  • Curvature for very light/heavy cars,

then a simple straight line may be too naïve for some ranges.

Residual distribution: mpg model

Code 
plt.hist(mpg["resid_mpg"], bins=20)
plt.xlabel("Residual")
plt.ylabel("Count")
plt.title("Residual distribution: mpg ~ weight")
plt.tight_layout()
plt.show()

Histogram of residuals for mpg model.

Warning

Again: we’re looking for rough symmetry and no wild tails.
But remember:

Even a nice residual histogram doesn’t fix a systematically wrong relationship.

Summary

❌ Anti-example: misusing linear models

Two common misuses in high-impact settings:

  1. Non-linear relationship forced into a line

    • For example, fuel efficiency might flatten out beyond a certain vehicle weight.
    • Residual plot shows a clear U-shape or S-shape.

  2. Wrong data type / units in serious decisions

    • Using a numeric code (e.g., pollution monitoring station ID) as if it were a meaningful continuous feature.
    • Year encoded as 1970, 1980, … but treated as “distance” in a way that doesn’t match the scientific question.

When this happens:

  • OLS will still give you coefficients and R²…
  • …but the geometry story is misaligned with reality: the real signal doesn’t live in the subspace spanned by your chosen features.

Regression metrics & decision context

From earlier ML modules (M38), for regression we often use:

  • MAE (mean absolute error)
    • “Typical absolute error”
  • RMSE (root mean squared error)
    • Penalizes large errors more; in same units as outcome
  • R² (coefficient of determination)
    • Fraction of variance in \(y\) explained by the model

Example (mpg):

Code 
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

mae = mean_absolute_error(y_m, mpg["pred_mpg"])
rmse = np.sqrt(mean_squared_error(y_m, mpg["pred_mpg"]))
r2 = r2_score(y_m, mpg["pred_mpg"])

print("MAE :", mae)
print("RMSE:", rmse)
print("R^2 :", r2)
MAE : 3.2787024115118464
RMSE: 4.321645126270701
R^2 : 0.6926304331206254

Warning

But remember:

Good metrics do not fix broken assumptions or bad causal reasoning.

Visualizing before you trust a model

Code 
fig, axes = plt.subplots(2, 2, figsize=(10, 8))

# 1. Data + fit
axes[0, 0].scatter(mpg["weight"], mpg["mpg"], alpha=0.6)
order = np.argsort(mpg["weight"].values)
axes[0, 0].plot(mpg["weight"].values[order], mpg["pred_mpg"].values[order])
axes[0, 0].set_xlabel("Weight")
axes[0, 0].set_ylabel("mpg")
axes[0, 0].set_title("mpg vs weight (fit)")

# 2. Residuals vs fitted
axes[0, 1].scatter(mpg["pred_mpg"], mpg["resid_mpg"], alpha=0.6)
axes[0, 1].axhline(0, linestyle="--")
axes[0, 1].set_xlabel("Fitted mpg")
axes[0, 1].set_ylabel("Residual")
axes[0, 1].set_title("Residuals vs fitted")

# 3. Residuals vs weight
axes[1, 0].scatter(mpg["weight"], mpg["resid_mpg"], alpha=0.6)
axes[1, 0].axhline(0, linestyle="--")
axes[1, 0].set_xlabel("Weight")
axes[1, 0].set_ylabel("Residual")
axes[1, 0].set_title("Residuals vs weight")

# 4. Histogram
axes[1, 1].hist(mpg["resid_mpg"], bins=15)
axes[1, 1].set_xlabel("Residual")
axes[1, 1].set_ylabel("Count")
axes[1, 1].setTitle = ("Residual histogram")

plt.tight_layout()
plt.show()

Four basic diagnostic views for mpg model.

These four views give a compact visual check of OLS assumptions.

Diagnostics checklist (high-impact decisions)

Before you brag about an \(R^2\) value (in any high-impact domain), ask:

  1. Linearity:
    • Do residuals vs fitted show a cloud or a pattern?
  2. Homoskedasticity:
    • Is the vertical spread roughly constant?
  3. Outliers & leverage:
    • Any points dominating the fit?
  4. Data collection & leakage:
    • Are any features “peeking at the future” or encoding sensitive info?
  5. Causal questions:
    • Are you predicting or trying to claim cause and effect?

These are lightweight but powerful habits.

Ethics & risks: medical & climate contexts

OLS feels simple, but in high-impact domains:

  • Health care:
    • A mis-specified regression for disease progression can under-estimate risk for certain groups.
    • Including variables that act as proxies for protected attributes (e.g., postal code, income) can encode structural inequities.
  • Climate & energy:
    • Poorly modeled fuel efficiency or emissions can underestimate environmental impact.
    • Overconfident predictions can be used to justify weak policy decisions.

Mitigations at this level:

  • Think carefully about what features you include and what they proxy.
  • Separate prediction from intervention claims (“if we do X, Y will change”).
  • Report limitations and uncertainty, not just a single R² number.

Quick check

Q1: Which statement best captures the geometric view of OLS?

A. OLS chooses coefficients so that the predictions equal the observed values exactly.
B. OLS finds the projection of \(y\) onto the column space of \(X\).
C. OLS chooses coefficients that minimize the absolute residuals.
D. OLS finds the largest possible R² on the test set.

Q: You fit a linear model for fuel efficiency and plot residuals vs fitted values. You see a clear U-shaped pattern.

1–2 sentences: what does this suggest, and what might you do?

Q3: Fill in the blanks: The OLS estimate \(\hat\beta\) satisfies the normal equations \[ \_\_\_\_\_\_\_\_ X\hat\beta = X^\top y. \]

Q4: You have two regression models for disease progression:

  • Model A: lower RMSE but slightly worse R²
  • Model B: higher R² but RMSE is much larger, and clinicians care most about typical prediction error in the original units.

Which model is more appropriate, and why?