5.5 Ordinary least squares and lm

Ordinary least squares (OLS) and the lm function provide a fast and widely used frequentist approach to linear regression. This section shows how the Bayesian analogue relates to OLS and when the two approaches agree.

5.5.1 The OLS estimator

The OLS estimator minimizes the sum of squared residuals:

$\hat{\beta}_{\text{OLS}} = (X^T X)^{-1} X^T y$

This is the maximum likelihood estimate under the Gaussian model. It is also the posterior mode (MAP) when the priors are flat (uniform).

# Demonstrate OLS and its Bayesian equivalent
set.seed(42)
n <- 100
x <- runif(n, 0, 10)
y <- 2 + 1.5 * x + rnorm(n, sd = 2)

# OLS fit

fit <- lm(y ~ x)
cat("OLS estimates:\n")
cat(" Intercept:", round(coef(fit)[1], 3), "\n")
cat(" Slope:", round(coef(fit)[2], 3), "\n")
cat(" Residual SD:", round(summary(fit)$sigma, 3), "\n")

# OLS as matrix computation

X <- cbind(1, x)
beta*ols <- solve(t(X) %*% X) %_% t(X) %*% y
cat("\nMatrix computation:\n")
cat(" Intercept:", round(beta_ols[1], 3), "\n")
cat(" Slope:", round(beta_ols[2], 3), "\n")

# They are identical

cat("\nBoth methods give the same answer.\n")

Click Run to execute this code in your browser with WebR

5.5.2 Bayesian vs frequentist

The Bayesian model with flat priors gives the same point estimates as OLS. The difference is that the Bayesian model also provides the full posterior distribution, not just point estimates and standard errors.

# Compare Bayesian (flat prior) with OLS
set.seed(42)
n <- 100
x <- runif(n, 0, 10)
y <- 2 + 1.5 * x + rnorm(n, sd = 2)

# OLS

fit <- lm(y ~ x)

# Bayesian posterior (flat prior = uniform)

# The posterior for (alpha, beta) is multivariate normal:

# (alpha, beta) ~ Normal(beta_OLS, sigma^2 * (X'X)^{-1})

library(MASS)
post <- mvrnorm(1e4, mu = coef(fit), Sigma = vcov(fit))

cat("Bayesian posterior (flat prior):\n")
cat(" alpha: mean =", round(mean(post[,1]), 3),
", 95% CI = [", round(quantile(post[,1], 0.025), 3), ",",
round(quantile(post[,1], 0.975), 3), "]\n")
cat(" beta: mean =", round(mean(post[,2]), 3),
", 95% CI = [", round(quantile(post[,2], 0.025), 3), ",",
round(quantile(post[,2], 0.975), 3), "]\n")

cat("\nOLS confidence intervals:\n")
conf <- confint(fit)
cat(" alpha: [", round(conf[1,1], 3), ",", round(conf[1,2], 3), "]\n")
cat(" beta: [", round(conf[2,1], 3), ",", round(conf[2,2], 3), "]\n")

cat("\nWith flat priors, Bayesian and frequentist intervals coincide.\n")

Click Run to execute this code in your browser with WebR

5.5.3 When priors matter

With informative priors, the Bayesian model differs from OLS. The priors pull the estimates toward the prior mean, which is especially useful with small samples or collinear predictors.

Rethinking

OLS is a special case, not the default. OLS is the Bayesian model with flat priors. Flat priors are not "objective" — they assign equal weight to absurd and reasonable values. In most cases, weakly informative priors are better: they regularize the estimates and prevent overfitting. The lm function is fast and convenient, but it is a special case of the Bayesian model, not the other way around.

Summary of Section 5.5

OLS minimizes the sum of squared residuals: $\hat{\beta} = (X^T X)^{-1} X^T y$ .
OLS is the MAP estimate of the Bayesian model with flat priors.
With flat priors, Bayesian and frequentist intervals coincide.
Informative priors make the Bayesian model differ from OLS, usually for the better.
OLS is a special case of the Bayesian model, not the default.

Categorical variables

Summary