Moving Beyond Linearity

Polynomials

\[y_i = \beta_0 + \beta_1x_i + \beta_2x_i^2+\beta_3x_i^3 \dots + \beta_dx_i^d+\epsilon_i\]

ggplot(co, aes(age, pop)) + 
  geom_point() + 
  geom_smooth(formula = y ~ poly(x, 4), method = "lm")

How is it done?

• New variables are created ( \(X_1 = X\), \(X_2 = X^2\), \(X_3 = X^3\), etc) and treated as multiple linear regression
• We are not interested in the individual coefficients, we are interested in how a specific \(x\) value behaves
• \(\hat{f}(x_0) = \hat\beta_0 + \hat\beta_1x_0 + \hat\beta_2x_0^2 + \hat\beta_3x_0^3 + \hat\beta_4x_0^4\)
• or more often a change between two values, \(a\) and \(b\)
• \(\hat{f}(b) -\hat{f}(a) = \hat\beta_1b + \hat\beta_2b^2 + \hat\beta_3b^3 + \hat\beta_4b^4 - \hat\beta_1a - \hat\beta_2a^2 - \hat\beta_3a^3 -\hat\beta_4a^4\)
• \(\hat{f}(b) -\hat{f}(a) =\hat\beta_1(b-a) + \hat\beta_2(b^2-a^2)+\hat\beta_3(b^3-a^3)+\hat\beta_4(b^4-a^4)\)

Polynomial Regression

\[\hat{f}(b) -\hat{f}(a) =\hat\beta_1(b-a) + \hat\beta_2(b^2-a^2)+\hat\beta_3(b^3-a^3)+\hat\beta_4(b^4-a^4)\]

• If given no other information, a sensible choice may be the 25th and 75th percentiles of \(x\)

Choosing \(d\)

\[y_i = \beta_0 + \beta_1x_i + \beta_2x_i^2+\beta_3x_i^3 \dots + \beta_dx_i^d+\epsilon_i\]

Either:

Step functions

• Create dummy variables for each group
• Include each of these variables in multiple regression
• The choice of cutpoints or knots can be problematic (and make a big difference!)

Step functions

\[C_1(X) = I(X < 35), C_2(X) = I(35\leq X<65), C_3(X) = I(X \geq 65)\]

mod <- lm(pop ~ I(age < 35) + I(age >=35 & age < 65) + I(age >= 65), data = co)
p <- predict(mod)
ggplot(co, aes(x = age, y = pop)) +
  geom_point() +
  geom_line(aes(x = age, y = p), color = "blue")

Step functions

\[C_1(X) = I(X < 15), C_2(X) = I(15\leq X<65), C_3(X) = I(X \geq 65)\]

mod <- lm(pop ~ I(age < 15) + I(age >=15 & age < 65) + I(age >= 65), data = co)
p <- predict(mod)
ggplot(co, aes(x = age, y = pop)) +
  geom_point() +
  geom_line(aes(x = age, y = p), color = "blue")

Piecewise polynomials

Instead of a single polynomial in \(X\) over it’s whole domain, we can use different polynomials in regions defined by knots

\[y_i = \begin{cases}\beta_{01}+\beta_{11}x_i + \beta_{21}x^2_i+\beta_{31}x^3_i+\epsilon_i& \textrm{if } x_i < c\\ \beta_{02}+\beta_{12}x_i + \beta_{22}x_i^2 + \beta_{32}x_{i}^3+\epsilon_i&\textrm{if }x_i\geq c\end{cases}\]

• It would be nice to have constraints (like continuity!)
• Insert splines!

Linear Splines

A linear spline with knots at \(\xi_k\), \(k = 1,\dots, K\) is a piecewise linear polynomial continuous at each knot

\[y_i = \beta_0 + \beta_1b_1(x_i)+\beta_2b_2(x_i)+\dots+\beta_{K+1}b_{K+1}(x_i)+\epsilon_i\]

• \(b_k\) are basis functions
• \(\begin{align}b_1(x_i)&=x_i\\ b_{k+1}(x_i)&=(x_i-\xi_k)_+,k=1,\dots,K\end{align}\)
• Here \(()_+\) means the positive part
• \((x_i-\xi_k)_+=\begin{cases}x_i-\xi_k & \textrm{if } x_i>\xi_k\\0&\textrm{otherwise}\end{cases}\)

Application Exercise

Let’s create data set to fit a linear spline with 2 knots: 35 and 65.

Application Exercise

Below is a linear regression model fit to include the 3 bases you just created with 2 knots: 35 and 65. Use the information here to draw the relationship between \(x\) and \(y\).

set.seed(1)
d <- tibble(
  b1 = c(4, 15, 25, 37, 49, 66, 70, 80),
  b2 = ifelse(b1 < 35, 0, b1 - 35),
  b3 = ifelse(b1 < 65, 0, b1 - 65),
  y = 2 * b1 + -2 * b2 -3 * b3 + rnorm(8, sd = 0.25)
)
lm(y ~ b1 + b2 + b3, data = d) |>
  tidy() |>
  knitr::kable(digits = 1)

Linear Splines

Linear Splines

p_ <- predict(lm(y ~ b1 + b2 + b3, data = d))

ggplot(d, aes(x = b1, y = p_)) +
  geom_point() +
  labs(x = "X",
       y = expression(hat(y)))

Linear Splines

newdat <- tibble(
  b1 = 4:80,
  b2 = ifelse(b1 > 35, b1 - 35, 0),
  b3 = ifelse(b1 > 65, b1 - 65, 0)
)
p <- predict(lm(y ~ b1 + b2 + b3, data = d),
             newdata = newdat)

ggplot(newdat, aes(x = b1, y = p)) +
  geom_point() + 
  labs(x = "X",
       y = expression(hat(y)))

Cubic Splines

A cubic splines with knots at \(\xi_i, k = 1, \dots, K\) is a piecewise cubic polynomial with continuous derivatives up to order 2 at each knot.

Again we can represent this model with truncated power functions

\[y_i = \beta_0 + \beta_1b_1(x_i)+\beta_2b_2(x_i)+\dots+\beta_{K+3}b_{K+3}(x_i) + \epsilon_i\]

\[\begin{align}b_1(x_i)&=x_i\\b_2(x_i)&=x_i^2\\b_3(x_i)&=x_i^3\\b_{k+3}(x_i)&=(x_i-\xi_k)^3_+, k = 1,\dots,K\end{align}\]

where

\[(x_i-\xi_k)^{3}_+=\begin{cases}(x_i-\xi_k)^3&\textrm{if }x_i>\xi_k\\0&\textrm{otherwise}\end{cases}\]

Application Exercise

Let’s create data set to fit a cubic spline with 2 knots: 35 and 65.

Cubic Spline

lm(y ~ b1 + b2 + b3 + b4 + b5, data = d) |>
  tidy() |>
  knitr::kable(digits = 3)

Cubic Splines

p_ <- predict(lm(y ~ b1 + b2 + b3 + b4 + b5, data = d))

ggplot(d, aes(x = b1, y = p_)) +
  geom_point() +
  labs(x = "X",
       y = expression(hat(y)))

Cubic Splines

newdat <- tibble(
  b1 = 4:80,
  b2 = b1^2,
  b3 = b1^3,
  b4 = ifelse(b1 > 35, (b1 - 35)^3, 0),
  b5 = ifelse(b1 > 65, (b1 - 65)^3, 0)
)
p <- predict(lm(y ~ b1 + b2 + b3 + b4 + b5, data = d),
             newdata = newdat)

ggplot(newdat, aes(x = b1, y = p)) +
  geom_point() + 
  labs(x = "X",
       y = expression(hat(y)))

Cubic Splines

newdat <- tibble(
  b1 = -100:100,
  b2 = b1^2,
  b3 = b1^3,
  b4 = ifelse(b1 > 35, (b1 - 35)^3, 0),
  b5 = ifelse(b1 > 65, (b1 - 65)^3, 0)
)
p <- predict(lm(y ~ b1 + b2 + b3 + b4 + b5, data = d),
             newdata = newdat)

ggplot(newdat, aes(x = b1, y = p)) +
  geom_point() + 
  geom_vline(xintercept = c(4, 80), lty = 2) +
  labs(x = "X",
       y = expression(hat(y)))

Natural cubic splines

A natural cubic spline extrapolates linearly beyond the boundary knots

This adds 4 extra constraints and allows us to put more internal knots for the same degrees of freedom as a regular cubic spline

Natural Cubic Splines

p <- predict(lm(y ~ splines::ns(b1, knots = c(35, 65)), data = d),
             newdata = newdat)

ggplot(newdat, aes(x = b1, y = p)) +
  geom_point() + 
  geom_vline(xintercept = c(4, 80), lty = 2) +
  labs(x = "X",
       y = expression(hat(y)))

Natural Cubic Splines

da <- tibble(
  x = newdat$b1,
  ns = predict(lm(y ~ splines::ns(b1, knots = c(35, 65)), data = d),
             newdata = newdat),
  cubic = predict(lm(y ~ b1 + b2 + b3 + b4 + b5, data = d), 
                  newdata = newdat),
  linear = predict(lm(y ~ b1 + ifelse(b1>35, b1 - 35, 0) + ifelse(b1>65, b1 - 65, 0), data = d),
                   newdata = newdat)
) |>
  pivot_longer(ns:linear)

da |>
  filter(name != "linear") |>
ggplot(aes(x = x, y = value, color = name)) +
  geom_point(alpha = 0.5) + 
  geom_vline(xintercept = c(4, 80), lty = 2) +
  labs(x = "X",
       y = expression(hat(y)),
       color = "Spline")

Natural Splines

ggplot(da, aes(x = x, y = value, color = name)) +
  geom_point(alpha = 0.5) + 
  geom_vline(xintercept = c(4, 80), lty = 2) +
  labs(x = "X",
       y = expression(hat(y)),
       color = "Spline")

Knot placement

• One strategy is to decide \(K\) (the number of knots) in advance and then place them at appropriate quantiles of the observed \(X\)
• A cubic spline with \(K\) knots has \(K+3\) parameters (or degrees of freedom!)
• A natural spline with \(K\) knots has \(K-1\) degrees of freedom

Knot placement

Here is a comparison of a degree-14 polynomial and natural cubic spline (both have 15 degrees of freedom)