Statistical Rethinking

4.4 Adding a predictor

Adding a predictor to the mean turns a simple Gaussian model into a linear regression. The expected value of the outcome becomes a linear function of the predictor, while the observations scatter around this expectation with residual variation.

4.4.1 The linear model

The linear regression model for height with weight as a predictor is:

$h_i \sim \text{Normal}(\mu_i, \sigma)$ $\mu_i = \alpha + \beta x_i$ $\alpha \sim \text{Normal}(178, 20)$ $\beta \sim \text{Log-Normal}(0, 1)$ $\sigma \sim \text{Uniform}(0, 50)$

The prior for $\beta$ is log-normal, which ensures that $\beta > 0$ (taller people weigh more). The log-normal prior has median 1, meaning that we expect about 1 cm of height per kg of weight.

# Simulate the Howell data with height and weight
set.seed(42)
n <- 352
weight <- rnorm(n, mean = 35, sd = 5)
height <- 100 + 2.5 * weight + rnorm(n, sd = 4)

# Plot the data

plot(weight, height, pch = 16, col = "steelblue",
main = "Howell data: height vs weight",
xlab = "Weight (kg)", ylab = "Height (cm)")

# Fit linear regression

fit <- lm(height ~ weight)
abline(fit, lwd = 2, col = "firebrick")
cat("Intercept (alpha):", round(coef(fit)[1], 2), "\n")
cat("Slope (beta):", round(coef(fit)[2], 2), "\n")
cat("R-squared:", round(summary(fit)$r.squared, 4), "\n")

Click Run to execute this code in your browser with WebR

4.4.2 Prior predictive simulation

Before fitting, we check what the priors imply about the relationship.

# Prior predictive simulation for the linear regression
set.seed(42)

# Priors

n_draws <- 50
alpha_prior <- rnorm(n_draws, 178, 20)
beta_prior <- rlnorm(n_draws, 0, 1) # log-normal: always positive

# Plot prior predictive lines

plot(NULL, xlim = c(0, 100), ylim = c(0, 300),
xlab = "Weight (kg)", ylab = "Height (cm)",
main = "Prior predictive: 50 lines from the prior")
for (i in 1:n_draws) {
abline(a = alpha_prior[i], b = beta_prior[i],
col = rgb(0.27, 0.51, 0.71, 0.2))
}
abline(h = 0, lty = 2, col = "gray50")
abline(h = 272, lty = 2, col = "gray50") # tallest person ever

cat("Prior for alpha: Normal(178, 20)\n")
cat(" Range:", round(range(alpha_prior), 1), "\n")
cat("Prior for beta: Log-Normal(0, 1)\n")
cat(" Range:", round(range(beta_prior), 1), "\n")
cat("Median of beta prior:", round(median(beta_prior), 2), "\n")
cat("\nSome prior lines are unreasonable (e.g., 300 cm at 50 kg)\n")
cat("But the data will constrain the posterior\n")

Click Run to execute this code in your browser with WebR

4.4.3 Fitting the model

We fit the model using quadratic approximation. The posterior gives us the joint distribution of $\alpha$ , $\beta$ , and $\sigma$ .

# Fit the linear regression using lm (quadratic approximation equivalent)
set.seed(42)
n <- 352
weight <- rnorm(n, mean = 35, sd = 5)
height <- 100 + 2.5 * weight + rnorm(n, sd = 4)

fit <- lm(height ~ weight)

# Posterior means and standard deviations

cat("Posterior means (quadratic approximation):\n")
cat("alpha:", round(coef(fit)[1], 3), "\n")
cat("beta:", round(coef(fit)[2], 3), "\n")
cat("sigma:", round(summary(fit)$sigma, 3), "\n")

# 95% credible intervals (approximate)

conf <- confint(fit)
cat("\n95% intervals:\n")
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")

Click Run to execute this code in your browser with WebR

4.4.4 Interpreting the posterior

The posterior for $\beta$ tells us how much height changes per unit of weight. The posterior for $\alpha$ tells us the expected height when weight is zero (which is not meaningful — it is an extrapolation).

# Interpret the posterior
set.seed(42)
n <- 352
weight <- rnorm(n, mean = 35, sd = 5)
height <- 100 + 2.5 * weight + rnorm(n, sd = 4)

fit <- lm(height ~ weight)

# Posterior samples (using the multivariate normal approximation)

library(MASS)
post_samples <- mvrnorm(1e4, mu = coef(fit), Sigma = vcov(fit))
alpha_samples <- post_samples[, 1]
beta_samples <- post_samples[, 2]
sigma_samples <- summary(fit)$sigma * sqrt((n - 2) / rchisq(1e4, n - 2))

# Posterior summaries

cat("Posterior for beta (cm per kg):\n")
cat(" Mean:", round(mean(beta_samples), 3), "\n")
cat(" 95% interval: [", round(quantile(beta_samples, 0.025), 3), ",",
round(quantile(beta_samples, 0.975), 3), "]\n")

# Plot posterior for beta

hist(beta_samples, breaks = 50, col = "steelblue", border = "white",
freq = FALSE, xlab = "beta (cm per kg)",
main = "Posterior distribution of beta")
abline(v = quantile(beta_samples, c(0.025, 0.975)), lty = 2, col = "firebrick", lwd = 2)

Click Run to execute this code in your browser with WebR

4.4.5 Posterior predictive lines

We can visualize the uncertainty in the regression line by plotting many lines sampled from the posterior.

# Posterior predictive lines
set.seed(42)
n <- 352
weight <- rnorm(n, mean = 35, sd = 5)
height <- 100 + 2.5 * weight + rnorm(n, sd = 4)

fit <- lm(height ~ weight)
library(MASS)
post_samples <- mvrnorm(100, mu = coef(fit), Sigma = vcov(fit))

# Plot data and posterior lines

plot(weight, height, pch = 16, col = rgb(0.27, 0.51, 0.71, 0.5),
main = "Posterior predictive lines",
xlab = "Weight (kg)", ylab = "Height (cm)")
for (i in 1:100) {
abline(a = post_samples[i, 1], b = post_samples[i, 2],
col = rgb(0.8, 0.2, 0.2, 0.1))
}
abline(fit, lwd = 2, col = "darkred")

Click Run to execute this code in your browser with WebR

Rethinking

The intercept is usually meaningless. The intercept $\alpha$ is the expected value of the outcome when the predictor is zero. But zero weight is not a meaningful value — it is an extrapolation far outside the data. To make the intercept meaningful, center the predictor: subtract the mean weight from each observation. Then $\alpha$ becomes the expected height at the average weight, which is a meaningful quantity. Centering also makes the posterior of $\alpha$ and $\beta$ less correlated, which improves computation.

Overthinking

Centering and standardization. Centering the predictor ( $x_i - \bar{x}$ ) makes the intercept interpretable. Standardization ( $\frac{x_i - \bar{x}}{s_x}$ ) goes further: it makes the slope interpretable as the expected change in $y$ per standard deviation of $x$ . This is useful for comparing the effects of different predictors measured on different scales. In the rethinking package, the function standardize() does this automatically.

Summary of Section 4.4

Adding a predictor turns the Gaussian model into a linear regression: $\mu_i = \alpha + \beta x_i$ .
The log-normal prior for $\beta$ ensures a positive relationship (taller people weigh more).
Prior predictive simulation checks whether the priors imply reasonable regression lines.
The posterior for $\beta$ tells us how much the outcome changes per unit of the predictor.
Centering the predictor makes the intercept interpretable and improves computation.
Posterior predictive lines visualize the uncertainty in the regression relationship.

A Gaussian model of height

Polynomial regression