Penalised regression splines from scratch

GAMs

smoothing

regression

A cubic B-spline basis with a difference penalty, coded by hand and checked against mgcv, shows why the penalty, not the polynomial degree, decides how wiggly a smooth becomes.

Author

Tidy Ecology

Published

2026-08-07

Fitting a smooth curve to an ecological gradient looks like a choice of flexibility: a straight line is too stiff, a high-order polynomial bends to every point. Neither extreme is satisfying. A straight line misses real structure; a flexible polynomial chases noise and behaves wildly near the edges of the data. The usual reflex is to tune the polynomial degree, but degree is a coarse dial, and every step adds global wiggliness that shows up worst where you have least data.

Penalised regression splines take a different route. You lay down a generous set of local basis functions, more than you think you need, and then discourage the fitted coefficients from varying too quickly between neighbours. A single continuous penalty replaces the discrete degree, and the data decide how much of the available flexibility to spend. This post builds the machinery from a B-spline basis and a difference penalty, selects the penalty by generalised cross-validation, and confirms the hand-coded fit against mgcv::gam.

A basis plus a penalty

We simulate a response along a gradient x on the unit interval. The truth is a growing-amplitude cycle, smooth but not a polynomial, so a fixed-degree polynomial is always slightly wrong.

set.seed(4206)
n <- 200
x <- sort(runif(n))
f <- function(x) sin(2 * pi * x) * (0.8 + 0.9 * x)
sigma <- 0.45
y <- f(x) + rnorm(n, 0, sigma)

## cubic B-spline basis with K functions, equally spaced knots
K <- 20; deg <- 3
rng <- range(x)
nik <- K - deg - 1
dx  <- diff(rng) / (nik + 1)
knots <- c(rng[1] - dx * (deg:1),
           seq(rng[1], rng[2], length = nik + 2),
           rng[2] + dx * (1:deg))
B <- splineDesign(knots, x, ord = deg + 1, outer.ok = TRUE)

## second-order difference penalty on adjacent coefficients
D <- diff(diag(K), differences = 2)
S <- crossprod(D)

The basis B has 20 columns, one per local bump. The penalty matrix S measures how fast the coefficients change from one bump to the next. Fitting minimises the usual residual sum of squares plus a multiple of that roughness, so the smoothing parameter lambda trades fit against smoothness continuously.

BtB <- crossprod(B); Bty <- crossprod(B, y)
fit_ps <- function(lam) {
  A  <- BtB + lam * S
  Ai <- solve(A)
  beta <- Ai %*% Bty
  fit  <- as.vector(B %*% beta)
  H    <- B %*% Ai %*% t(B)          # influence matrix
  edf  <- sum(diag(H))               # effective degrees of freedom
  rss  <- sum((y - fit)^2)
  gcv  <- n * rss / (n - edf)^2      # generalised cross-validation
  list(beta = beta, fit = fit, edf = edf, gcv = gcv)
}

## select lambda by GCV over a log-spaced grid
lgrid <- 10^seq(-4, 6, length = 200)
gcvv  <- sapply(lgrid, function(l) fit_ps(l)$gcv)
p6_lam    <- lgrid[which.min(gcvv)]
ps_best   <- fit_ps(p6_lam)
ps_unpen  <- fit_ps(1e-8)
p6_edf    <- ps_best$edf
p6_edfun  <- ps_unpen$edf

Generalised cross-validation picks a penalty of 7.49. At that penalty the fit spends 7.1 effective degrees of freedom, well short of the 20 it would use with the penalty switched off. Effective degrees of freedom count how many free parameters the fit behaves as if it has: the penalty has folded twenty basis functions down to about seven worth of flexibility.

Scatter plot of the simulated data with three curves: a dotted line for the true function, an orange polynomial curve that swings sharply upward at the right edge, and a green penalised-spline curve that stays close to the dotted line throughout. — Figure 1: A degree-12 polynomial (orange) overshoots near both edges, while the penalised spline (green) tracks the underlying curve (dotted). Both fit the interior points, but the polynomial pays for its global flexibility at the boundaries.

Does the hand-coded fit match mgcv?

The point of coding the penalty by hand is to see that nothing is hidden. To check it, we fit the same P-spline through mgcv, which places the knots slightly differently and applies its own identifiability constraint, then compare fitted values.

d <- data.frame(x = x, y = y)
m_gcv  <- gam(y ~ s(x, bs = "ps", k = 20, m = c(2, 2)), data = d, method = "GCV.Cp")
m_reml <- gam(y ~ s(x, bs = "ps", k = 20, m = c(2, 2)), data = d, method = "REML")

p6_edfmg   <- sum(m_gcv$edf)                 # total, including intercept
p6_edfreml <- summary(m_reml)$edf            # smooth term, REML
p6_maxdiff <- max(abs(ps_best$fit - as.vector(predict(m_gcv))))

The two fits agree to 0.0035 at every point: the hand-coded total of 7.1 effective degrees of freedom sits alongside mgcv’s 7.2. Swapping generalised cross-validation for restricted maximum likelihood, the penalty selection recommended by Wood (2011) because it is less prone to under-smoothing, gives a slightly smoother 7.4 degrees of freedom for the smooth term. The machinery is the same either way; only the criterion for lambda differs.

Too stiff, too flexible, or penalised

The three fitting strategies fail and succeed in instructive ways. We compare each against the known truth by root mean squared error, and we ask what each does just outside the data, at x equal to 1.15.

poly4 <- lm(y ~ poly(x, 4)); poly12 <- lm(y ~ poly(x, 12))
truef <- f(x)
rmse  <- function(fit) sqrt(mean((fit - truef)^2))

p6_r_p4  <- rmse(fitted(poly4))
p6_r_p12 <- rmse(fitted(poly12))
p6_r_un  <- rmse(ps_unpen$fit)
p6_r_pen <- rmse(ps_best$fit)

xe <- 1.15
Be <- splineDesign(knots, xe, ord = deg + 1, outer.ok = TRUE)
p6_ext_pen <- as.vector(Be %*% ps_best$beta)
p6_ext_p12 <- as.vector(predict(poly12, data.frame(x = xe)))

Against the truth, the penalised spline reaches an error of 0.082, beating both the low-order polynomial (0.118) and the unpenalised spline (0.152). The ranking is the whole story in miniature. The degree-4 polynomial is too stiff and misses the shape; the unpenalised spline is too flexible and fits the noise; the penalty lands between them without anyone choosing where. The high-order polynomial is worst of all outside the data: at x of 1.15 it predicts -147, a runaway swing, while the spline stays bounded near 0.00.

A U-shaped curve of the GCV score against the log penalty, with a vertical dashed line marking the minimum near a penalty of about seven, annotated with the effective degrees of freedom at that point. — Figure 2: The generalised cross-validation score as a function of the penalty. The minimum (marked) selects the smoothing parameter; nobody sets the wiggliness by hand.

What the penalty does and does not buy

The penalty turns a discrete modelling decision into an estimated quantity. You no longer commit to a degree; you provide enough basis functions and let generalised cross-validation or restricted maximum likelihood set the effective flexibility. That removes the worst failure of global polynomials, the wild boundary behaviour, and it adapts the fit to the data at hand.

It does not make extrapolation safe. The penalised spline stayed bounded outside the data, but bounded is not correct: near the boundary it drifts towards a straight line, not towards the true continuation. The penalty controls behaviour inside the data, where there is information; outside, the fit follows the basis, not the process. Nothing here tells you how the curve behaves where you have no points, and the smoothness itself is only weakly identified from a single sample. Those limits are the subject of the rest of this cluster: how large the basis should be, what goes wrong when smooth terms compete for the same signal, and how to check a fitted smooth before trusting it.

References

Eilers PHC, Marx BD 1996. Statistical Science 11(2):89-121 (10.1214/ss/1038425655).

Hastie T, Tibshirani R 1986. Statistical Science 1(3):297-310 (10.1214/ss/1177013604).

Craven P, Wahba G 1979. Numerische Mathematik 31(4):377-403 (10.1007/BF01404567).

Wood SN 2011. Journal of the Royal Statistical Society Series B 73(1):3-36 (10.1111/j.1467-9868.2010.00749.x).

Wood SN 2017. 2nd edn. CRC Press. ISBN 978-1-4987-2833-1.

--- title: "Penalised regression splines from scratch" description: "A cubic B-spline basis with a difference penalty, coded by hand and checked against mgcv, shows why the penalty, not the polynomial degree, decides how wiggly a smooth becomes." date: "2026-08-07 09:00" categories: [GAMs, smoothing, regression] image: thumbnail.png image-alt: "Scatter of points along a gradient with a wiggly polynomial overshooting at the edges and a penalised spline tracking the underlying curve." --- Fitting a smooth curve to an ecological gradient looks like a choice of flexibility: a straight line is too stiff, a high-order polynomial bends to every point. Neither extreme is satisfying. A straight line misses real structure; a flexible polynomial chases noise and behaves wildly near the edges of the data. The usual reflex is to tune the polynomial degree, but degree is a coarse dial, and every step adds global wiggliness that shows up worst where you have least data. Penalised regression splines take a different route. You lay down a generous set of local basis functions, more than you think you need, and then discourage the fitted coefficients from varying too quickly between neighbours. A single continuous penalty replaces the discrete degree, and the data decide how much of the available flexibility to spend. This post builds the machinery from a B-spline basis and a difference penalty, selects the penalty by generalised cross-validation, and confirms the hand-coded fit against `mgcv::gam`. ## A basis plus a penalty We simulate a response along a gradient `x` on the unit interval. The truth is a growing-amplitude cycle, smooth but not a polynomial, so a fixed-degree polynomial is always slightly wrong. ```{r} #| label: setup #| echo: false #| message: false #| warning: false library(splines) library(mgcv) library(ggplot2) theme_te <- function() { theme_minimal(base_size = 12) + theme( panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e5e4dc", linewidth = 0.3), plot.background = element_rect(fill = "#f5f4ee", colour = NA), panel.background = element_rect(fill = "#f5f4ee", colour = NA), text = element_text(colour = "#16241d"), axis.text = element_text(colour = "#16241d"), plot.title = element_text(face = "bold") ) } col_truth <- "#16241d"; col_naive <- "#c0622d" col_pen <- "#275139"; col_sage <- "#93a87f" ``` ```{r} #| label: simulate-and-basis #| echo: true set.seed(4206) n <- 200 x <- sort(runif(n)) f <- function(x) sin(2 * pi * x) * (0.8 + 0.9 * x) sigma <- 0.45 y <- f(x) + rnorm(n, 0, sigma) ## cubic B-spline basis with K functions, equally spaced knots K <- 20; deg <- 3 rng <- range(x) nik <- K - deg - 1 dx <- diff(rng) / (nik + 1) knots <- c(rng[1] - dx * (deg:1), seq(rng[1], rng[2], length = nik + 2), rng[2] + dx * (1:deg)) B <- splineDesign(knots, x, ord = deg + 1, outer.ok = TRUE) ## second-order difference penalty on adjacent coefficients D <- diff(diag(K), differences = 2) S <- crossprod(D) ``` The basis `B` has `r ncol(B)` columns, one per local bump. The penalty matrix `S` measures how fast the coefficients change from one bump to the next. Fitting minimises the usual residual sum of squares plus a multiple of that roughness, so the smoothing parameter `lambda` trades fit against smoothness continuously. ```{r} #| label: penalised-fit #| echo: true BtB <- crossprod(B); Bty <- crossprod(B, y) fit_ps <- function(lam) { A <- BtB + lam * S Ai <- solve(A) beta <- Ai %*% Bty fit <- as.vector(B %*% beta) H <- B %*% Ai %*% t(B) # influence matrix edf <- sum(diag(H)) # effective degrees of freedom rss <- sum((y - fit)^2) gcv <- n * rss / (n - edf)^2 # generalised cross-validation list(beta = beta, fit = fit, edf = edf, gcv = gcv) } ## select lambda by GCV over a log-spaced grid lgrid <- 10^seq(-4, 6, length = 200) gcvv <- sapply(lgrid, function(l) fit_ps(l)$gcv) p6_lam <- lgrid[which.min(gcvv)] ps_best <- fit_ps(p6_lam) ps_unpen <- fit_ps(1e-8) p6_edf <- ps_best$edf p6_edfun <- ps_unpen$edf ``` Generalised cross-validation picks a penalty of `r sprintf("%.2f", p6_lam)`. At that penalty the fit spends `r sprintf("%.1f", p6_edf)` effective degrees of freedom, well short of the `r round(p6_edfun)` it would use with the penalty switched off. Effective degrees of freedom count how many free parameters the fit behaves as if it has: the penalty has folded twenty basis functions down to about seven worth of flexibility. ```{r} #| label: fig-fit #| echo: false #| fig-width: 7 #| fig-height: 4.5 #| fig-cap: "A degree-12 polynomial (orange) overshoots near both edges, while the penalised spline (green) tracks the underlying curve (dotted). Both fit the interior points, but the polynomial pays for its global flexibility at the boundaries." #| fig-alt: "Scatter plot of the simulated data with three curves: a dotted line for the true function, an orange polynomial curve that swings sharply upward at the right edge, and a green penalised-spline curve that stays close to the dotted line throughout." poly12 <- lm(y ~ poly(x, 12)) grid <- data.frame(x = seq(min(x), max(x), length = 400)) Bg <- splineDesign(knots, grid$x, ord = deg + 1, outer.ok = TRUE) grid$pen <- as.vector(Bg %*% ps_best$beta) grid$pl12 <- predict(poly12, grid) grid$tru <- f(grid$x) ggplot() + geom_point(data = data.frame(x, y), aes(x, y), colour = col_sage, alpha = 0.6, size = 1.4) + geom_line(data = grid, aes(x, tru), colour = col_truth, linetype = "dotted", linewidth = 0.8) + geom_line(data = grid, aes(x, pl12), colour = col_naive, linewidth = 0.9) + geom_line(data = grid, aes(x, pen), colour = col_pen, linewidth = 1) + coord_cartesian(ylim = c(-2.4, 2.4)) + labs(x = "gradient x", y = "response y", title = "Global polynomial versus penalised spline") + theme_te() ``` ## Does the hand-coded fit match mgcv? The point of coding the penalty by hand is to see that nothing is hidden. To check it, we fit the same P-spline through `mgcv`, which places the knots slightly differently and applies its own identifiability constraint, then compare fitted values. ```{r} #| label: mgcv-check #| echo: true d <- data.frame(x = x, y = y) m_gcv <- gam(y ~ s(x, bs = "ps", k = 20, m = c(2, 2)), data = d, method = "GCV.Cp") m_reml <- gam(y ~ s(x, bs = "ps", k = 20, m = c(2, 2)), data = d, method = "REML") p6_edfmg <- sum(m_gcv$edf) # total, including intercept p6_edfreml <- summary(m_reml)$edf # smooth term, REML p6_maxdiff <- max(abs(ps_best$fit - as.vector(predict(m_gcv)))) ``` The two fits agree to `r sprintf("%.4f", p6_maxdiff)` at every point: the hand-coded total of `r sprintf("%.1f", p6_edf)` effective degrees of freedom sits alongside mgcv's `r sprintf("%.1f", p6_edfmg)`. Swapping generalised cross-validation for restricted maximum likelihood, the penalty selection recommended by Wood (2011) because it is less prone to under-smoothing, gives a slightly smoother `r sprintf("%.1f", p6_edfreml)` degrees of freedom for the smooth term. The machinery is the same either way; only the criterion for `lambda` differs. ## Too stiff, too flexible, or penalised The three fitting strategies fail and succeed in instructive ways. We compare each against the known truth by root mean squared error, and we ask what each does just outside the data, at `x` equal to 1.15. ```{r} #| label: compare #| echo: true poly4 <- lm(y ~ poly(x, 4)); poly12 <- lm(y ~ poly(x, 12)) truef <- f(x) rmse <- function(fit) sqrt(mean((fit - truef)^2)) p6_r_p4 <- rmse(fitted(poly4)) p6_r_p12 <- rmse(fitted(poly12)) p6_r_un <- rmse(ps_unpen$fit) p6_r_pen <- rmse(ps_best$fit) xe <- 1.15 Be <- splineDesign(knots, xe, ord = deg + 1, outer.ok = TRUE) p6_ext_pen <- as.vector(Be %*% ps_best$beta) p6_ext_p12 <- as.vector(predict(poly12, data.frame(x = xe))) ``` Against the truth, the penalised spline reaches an error of `r sprintf("%.3f", p6_r_pen)`, beating both the low-order polynomial (`r sprintf("%.3f", p6_r_p4)`) and the unpenalised spline (`r sprintf("%.3f", p6_r_un)`). The ranking is the whole story in miniature. The degree-4 polynomial is too stiff and misses the shape; the unpenalised spline is too flexible and fits the noise; the penalty lands between them without anyone choosing where. The high-order polynomial is worst of all outside the data: at `x` of 1.15 it predicts `r sprintf("%.0f", p6_ext_p12)`, a runaway swing, while the spline stays bounded near `r sprintf("%.2f", p6_ext_pen)`. ```{r} #| label: fig-gcv #| echo: false #| fig-width: 7 #| fig-height: 4.5 #| fig-cap: "The generalised cross-validation score as a function of the penalty. The minimum (marked) selects the smoothing parameter; nobody sets the wiggliness by hand." #| fig-alt: "A U-shaped curve of the GCV score against the log penalty, with a vertical dashed line marking the minimum near a penalty of about seven, annotated with the effective degrees of freedom at that point." gdf <- data.frame(loglam = log10(lgrid), gcv = gcvv) ggplot(gdf, aes(loglam, gcv)) + geom_line(colour = col_pen, linewidth = 1) + geom_vline(xintercept = log10(p6_lam), linetype = "dashed", colour = col_naive, linewidth = 0.7) + annotate("text", x = log10(p6_lam), y = max(gcvv) * 0.96, label = sprintf("edf = %.1f", p6_edf), hjust = -0.1, colour = col_truth, size = 3.6) + labs(x = "log10 penalty", y = "GCV score", title = "The penalty is selected, not chosen") + theme_te() ``` ## What the penalty does and does not buy The penalty turns a discrete modelling decision into an estimated quantity. You no longer commit to a degree; you provide enough basis functions and let generalised cross-validation or restricted maximum likelihood set the effective flexibility. That removes the worst failure of global polynomials, the wild boundary behaviour, and it adapts the fit to the data at hand. It does not make extrapolation safe. The penalised spline stayed bounded outside the data, but bounded is not correct: near the boundary it drifts towards a straight line, not towards the true continuation. The penalty controls behaviour inside the data, where there is information; outside, the fit follows the basis, not the process. Nothing here tells you how the curve behaves where you have no points, and the smoothness itself is only weakly identified from a single sample. Those limits are the subject of the rest of this cluster: how large the basis should be, what goes wrong when smooth terms compete for the same signal, and how to check a fitted smooth before trusting it. ## References Eilers PHC, Marx BD 1996. Statistical Science 11(2):89-121 (10.1214/ss/1038425655). Hastie T, Tibshirani R 1986. Statistical Science 1(3):297-310 (10.1214/ss/1177013604). Craven P, Wahba G 1979. Numerische Mathematik 31(4):377-403 (10.1007/BF01404567). Wood SN 2011. Journal of the Royal Statistical Society Series B 73(1):3-36 (10.1111/j.1467-9868.2010.00749.x). Wood SN 2017. 2nd edn. CRC Press. ISBN 978-1-4987-2833-1. ## Related tutorials - [Choosing the basis dimension k in mgcv](../choosing-basis-dimension-k/) - [Concurvity in additive models](../concurvity-in-additive-models/) - [Checking a generalised additive model](../checking-a-gam/) - [Species response curves with GAMs](../gam-species-response-curves/)