4.5 Polynomial regression

Linear relationships are not always enough. Polynomial regression adds powers of a predictor to the model, allowing the expected outcome to curve. This section discusses how to choose the degree of the polynomial and why too much flexibility leads to overfitting.

4.5.1 The polynomial model

The polynomial regression model adds squared, cubed, or higher-order terms:

$\mu_i = \alpha + \beta_1 x_i + \beta_2 x_i^2 + \beta_3 x_i^3$

The quadratic term ( $\beta_2 x_i^2$ ) allows the relationship to curve. The cubic term ( $\beta_3 x_i^3$ ) allows for an inflection point. Higher-order terms add more flexibility but also more risk of overfitting.

# Polynomial regression: height vs weight
# The relationship is slightly nonlinear for the full age range

set.seed(42)
n <- 200
weight <- runif(n, 10, 65)

# True relationship: nonlinear (growth curve)

height <- 50 + 3 _ weight - 0.02 _ weight^2 + rnorm(n, sd = 3)

# Fit linear, quadratic, and cubic models

m1 <- lm(height ~ weight)
m2 <- lm(height ~ weight + I(weight^2))
m3 <- lm(height ~ weight + I(weight^2) + I(weight^3))

# Compare R-squared

cat("Linear: R-squared =", round(summary(m1)$r.squared, 4), "\n")
cat("Quadratic: R-squared =", round(summary(m2)$r.squared, 4), "\n")
cat("Cubic: R-squared =", round(summary(m3)$r.squared, 4), "\n")

# Plot

plot(weight, height, pch = 16, col = "steelblue",
main = "Polynomial regression: linear vs quadratic vs cubic",
xlab = "Weight (kg)", ylab = "Height (cm)")
x_seq <- seq(min(weight), max(weight), length.out = 100)
lines(x_seq, predict(m1, data.frame(weight = x_seq)), lwd = 2, col = "firebrick")
lines(x_seq, predict(m2, data.frame(weight = x_seq)), lwd = 2, col = "darkgreen")
lines(x_seq, predict(m3, data.frame(weight = x_seq)), lwd = 2, col = "purple")
legend("bottomright",
legend = c("Linear", "Quadratic", "Cubic"),
col = c("firebrick", "darkgreen", "purple"), lwd = 2)

Click Run to execute this code in your browser with WebR

4.5.2 Standardizing the predictor

When using polynomial terms, the predictor values can become very large (e.g., $65^3 = 274625$ ). This causes numerical problems. The solution is to standardize the predictor before computing powers:

$z_i = \frac{x_i - \bar{x}}{s_x}$

Then the model becomes:

$\mu_i = \alpha + \beta_1 z_i + \beta_2 z_i^2 + \beta_3 z_i^3$

# Standardize the predictor before polynomial regression
set.seed(42)
n <- 200
weight <- runif(n, 10, 65)
height <- 50 + 3 * weight - 0.02 * weight^2 + rnorm(n, sd = 3)

# Standardize weight

weight_z <- (weight - mean(weight)) / sd(weight)

# Fit models with standardized predictor

m1z <- lm(height ~ weight_z)
m2z <- lm(height ~ weight_z + I(weight_z^2))
m3z <- lm(height ~ weight_z + I(weight_z^2) + I(weight_z^3))

cat("With standardized predictor:\n")
cat("Linear: R-squared =", round(summary(m1z)$r.squared, 4), "\n")
cat("Quadratic: R-squared =", round(summary(m2z)$r.squared, 4), "\n")
cat("Cubic: R-squared =", round(summary(m3z)$r.squared, 4), "\n")

# Coefficients are now on a comparable scale

cat("\nQuadratic coefficients:\n")
cat("alpha:", round(coef(m2z)[1], 2), "\n")
cat("beta1:", round(coef(m2z)[2], 2), "\n")
cat("beta2:", round(coef(m2z)[3], 2), "\n")

Click Run to execute this code in your browser with WebR

4.5.3 The danger of overfitting

Higher-order polynomials can wiggle wildly between data points, fitting noise rather than signal. This is overfitting — the model fits the sample well but generalizes poorly.

# Demonstrate overfitting with high-degree polynomials
set.seed(42)
n <- 30
x <- runif(n, 0, 10)
y <- 2 + 1.5 * x + rnorm(n, sd = 2)

# Fit polynomials of increasing degree

par(mfrow = c(2, 2))
for (deg in c(1, 3, 5, 10)) {
fit <- lm(y ~ poly(x, degree = deg, raw = TRUE))

plot(x, y, pch = 16, col = "steelblue",
main = sprintf("Degree %d (R² = %.3f)", deg, summary(fit)$r.squared),
xlab = "x", ylab = "y", ylim = c(-5, 25))
x_seq <- seq(0, 10, length.out = 200)
lines(x_seq, predict(fit, data.frame(x = x_seq)), lwd = 2, col = "firebrick")
}
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

The degree-10 polynomial fits the data almost perfectly, but it wiggles wildly between data points. This model would generalize terribly to new data. The degree-1 (linear) model captures the true relationship without overfitting.

Rethinking

Simplicity is a virtue. The golem does not know when to stop adding parameters. It will happily fit a degree-20 polynomial if you tell it to. It is the modeler's job to choose the simplest model that captures the important features of the data. This is the principle of parsimony: use the least complexity necessary. Information criteria (Chapter 6) provide a principled way to balance fit and complexity.

Overthinking

Splines and natural polynomials. Polynomial regression is the simplest form of nonlinear regression, but it has limitations: high-degree polynomials are unstable at the edges of the data. Splines are piecewise polynomials that are smoother and more stable. A spline divides the range of the predictor into segments and fits a low-degree polynomial to each segment, joining them smoothly at knots. The rethinking package provides spline() and related functions. Splines are more flexible than polynomials and less prone to edge effects, but they require choosing the number and placement of knots.

Summary of Section 4.5

Polynomial regression adds powers of the predictor: $\mu_i = \alpha + \beta_1 x_i + \beta_2 x_i^2 + \ldots$
Standardize the predictor before computing powers to avoid numerical problems.
High-degree polynomials overfit: they wiggle wildly between data points.
Use the simplest model that captures the important features — parsimony is a virtue.
Splines are a more stable alternative to high-degree polynomials.

Adding a predictor

Summary