Penalised regression splines from scratch

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

Updated 26 September 2026: a new section, GCV or REML over many samples, refits both smoothing criteria to fresh samples of this curve and of a straight and an empty term, and measures how often each adds curvature the truth does not have and how much worse than the other each can fit.

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.

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"

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.

Figure 1: A degree-12 polynomial (orange) and the penalised spline (green) against the truth (dotted). Inside the data the two are close, the polynomial drifting furthest at the left edge, by about 0.25. Its real failure sits just outside this window, where it runs to -147 at an x of 1.15.

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_edfsm   <- summary(m_gcv)$edf             # smooth term only, GCV
p6_edfreml <- summary(m_reml)$edf            # smooth term only, 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, both counts including the intercept. Compare the smooth term on its own and the criteria separate: 6.2 degrees of freedom under generalised cross-validation against 7.4 under restricted maximum likelihood, the selection Wood (2011) recommends because its criterion has a better defined optimum. On this sample REML smooths a little less, not more. The machinery is the same either way; only the criterion for lambda differs.

GCV or REML over many samples

The comparison above rests on one sample, and gam() does not ask which criterion you want: without a method argument it uses generalised cross-validation, method = "GCV.Cp" in mgcv 1.9-4, the version this post was knitted with (the chunk below stops if a later version changes the default). To see whether “REML smooths a little less” belongs to the criterion or to this draw, the chunk refits both criteria to 400 fresh samples of the same curve, with the same n, noise standard deviation and P-spline basis. A second arm fits y ~ s(x) + s(z) to data in which x has a straight effect of slope 1.5 and z has none, the same basis for both terms. The second-order difference penalty leaves a straight line unpenalised, so a term fitted as a line has about one effective degree of freedom, and more than two means curvature the truth does not have.

stopifnot(identical(formals(gam)$method, "GCV.Cp"),   # the default measured here,
          gam(round(3 * exp(f(x))) ~ s(x), family = poisson)$method == "UBRE")  # UBRE for counts
crit_reps <- 400; crit_slope <- 1.5
crit_ps <- function(v) sprintf("s(%s, bs = 'ps', k = 20, m = c(2, 2))", v)
crit_pair <- function(fm, dat) list(g = gam(fm, data = dat, method = "GCV.Cp"),
                                    r = gam(fm, data = dat, method = "REML"))
crit_wiggly <- function() {                 # this post's curve, n and sigma
  dw <- data.frame(xw = sort(runif(n))); dw$yw <- f(dw$xw) + rnorm(n, 0, sigma)
  fp <- crit_pair(as.formula(paste("yw ~", crit_ps("xw"))), dw)
  c(edf_g = summary(fp$g)$edf, edf_r = summary(fp$r)$edf,
    rmse_g = sqrt(mean((fitted(fp$g) - f(dw$xw))^2)),
    rmse_r = sqrt(mean((fitted(fp$r) - f(dw$xw))^2)))
}
crit_straight <- function() {               # x straight, z no effect at all
  ds <- data.frame(xs = runif(n), zs = runif(n))
  ds$ys <- crit_slope * ds$xs + rnorm(n, 0, sigma)
  fp <- crit_pair(as.formula(paste("ys ~", crit_ps("xs"), "+", crit_ps("zs"))), ds)
  sg <- summary(fp$g)$s.table; sr <- summary(fp$r)$s.table
  c(lin_g = sg[1, "edf"], lin_r = sr[1, "edf"], nul_g = sg[2, "edf"],
    nul_r = sr[2, "edf"], p_g = sg[2, "p-value"], p_r = sr[2, "p-value"])
}
set.seed(14302)
crit_w <- as.data.frame(t(replicate(crit_reps, crit_wiggly())))
crit_s <- as.data.frame(t(replicate(crit_reps, crit_straight())))
## share of samples, with its Monte Carlo standard error
crit_share <- function(v) c(p = mean(v), se = sqrt(mean(v) * (1 - mean(v)) / length(v)))
crit_lin <- sapply(crit_s[c("lin_g", "lin_r")], function(e) crit_share(e > 2))
crit_nul <- sapply(crit_s[c("nul_g", "nul_r")], function(e) crit_share(e > 2))
crit_fp  <- sapply(crit_s[c("p_g", "p_r")], function(p) crit_share(p < 0.05))
crit_q95 <- sapply(crit_s[c("nul_g", "nul_r")], quantile, probs = 0.95)
crit_more <- crit_share(crit_w$edf_r > crit_w$edf_g)
crit_qw <- sapply(crit_w[c("edf_g", "edf_r")], quantile, probs = c(0.05, 0.5, 0.95))
crit_ratio <- crit_w$rmse_g / crit_w$rmse_r
crit_bad_g <- crit_share(crit_ratio > 1.25); crit_bad_r <- crit_share(crit_ratio < 0.8)
crit_bad_up <- (crit_w$edf_g > crit_w$edf_r)[crit_ratio > 1.25]

On the straight term GCV spends more than two effective degrees of freedom in 0.330 of samples and REML in 0.107; on the null term the shares are 0.315 and 0.107 (Monte Carlo standard errors between 0.015 and 0.024), so GCV bends a term that is straight or empty about three times as often as REML (the ratios of the two shares are 3.1 and 2.9). Its tail is long: the 95th percentile, over samples, of the null term’s effective degrees of freedom is 10.8 under GCV and 2.4 under REML. The p-values barely notice. The null term is called significant at 5 per cent in 0.062 of samples under GCV and 0.060 under REML, each with a standard error of about 0.012.

On this post’s wiggly curve the usual sample goes the other way. REML spends more degrees of freedom than GCV in 0.863 of samples (standard error 0.017), so the sample above, where REML smoothed a little less, was the usual case. The median smooth-term edf is 5.8 under GCV and 7.2 under REML; GCV’s 5th to 95th percentile runs from 5.3 to 10.3, REML’s from 6.8 to 7.7. Accuracy against the truth is almost level: mean root mean squared error 0.0823 for GCV and 0.0811 for REML, and a median ratio of GCV error to REML error of 0.99. The price of GCV’s spread is a tail, and what sets it apart is how far the tail reaches more than how often it is hit. In 0.058 of samples (standard error 0.012) its error is more than 1.25 times REML’s, and in 22 of those 23 samples GCV had spent more degrees of freedom than REML: the kind of failure Wood (2011) describes, a small number of severe undersmoothing fits under GCV. REML’s error exceeds 1.25 times GCV’s in 0.028 of samples (standard error 0.008), but never by much: its worst sample is 1.33 times GCV’s error, while GCV’s error is more than 1.5 times REML’s in 9 of the 400 samples and reaches 4.5 times.

Figure 2: Effective degrees of freedom of the smooth term under GCV against REML, one point per simulated sample and both criteria fitted to the same data. The dotted line marks equal edf. In the right panel, orange marks samples in which the GCV fit’s error against the truth is more than 1.25 times the REML fit’s.

GCV’s fault here is spread, not a consistent direction: it adds curvature to straight and empty terms far more often than REML, and on a curve that really bends it usually smooths more than REML, in exchange for a handful of fits much worse than REML’s on the same data, while REML is never much worse than GCV. Neither criterion smooths more on every curve. On a narrow seasonal pulse in Harmonics or cyclic splines for a narrow seasonal peak the typical fit runs the other way, REML smoothing more than UBRE (what method = "GCV.Cp" uses for Poisson counts, where the scale is known) and widening the crest, while UBRE’s occasional runaway spike is the same kind of tail. What method = "REML" buys is a steadier choice of penalty, not a smoother or rougher curve; since it is not the default, it has to be written in the call. This is one sample size, one noise level, one basis (a 20-function P-spline) and a Gaussian response.

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)))

## how much basis is left at and beyond the edge of the data
p6_row_in  <- sum(splineDesign(knots, max(x), ord = deg + 1, outer.ok = TRUE))
p6_row_mid <- sum(splineDesign(knots, 1.05,   ord = deg + 1, outer.ok = TRUE))
p6_row_out <- sum(Be)

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 and returns 0.00. Bounded is not the same as right, since the truth there is 1.48, and the next section says where that zero comes from.

Figure 3: 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 here means dead. A B-spline basis sums to one only across the interval its knots cover; past the last data knot the bumps run out one by one, so the design row that sums to 1.00 at the edge of the data sums to 0.87 at an x of 1.05 and 0.01 at 1.15. The fit is dragged to zero, which is why it returned 0.00 against a true 1.48. That is the basis ending rather than any statement about the process. The penalty controls behaviour inside the data, where there is information; outside, the fit follows the basis. 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 Generalized Additive Models: An Introduction with R, 2nd edn, CRC Press, ISBN 978-1-4987-2833-1.

Related tutorials

• Choosing the basis dimension k in mgcv
• Concurvity in additive models
• Checking a generalised additive model
• Tensor product smooths in mgcv

Reproducibility

• Code last run: 2026-09-30
• R version 4.5.3 (2026-03-11), aarch64-apple-darwin20
• Packages: ggplot2 4.0.3, mgcv 1.9-4