Bayesian Penalised Splines & Additive Models (R)

mgcv — the penalised-spline / GAM engine

The R counterpart to pspline_python.ipynb. Where the Python notebook builds the penalised spline from scratch (a random-walk-prior Gibbs) and cross-checks in PyMC, this notebook uses mgcv — the standard implementation of penalised regression splines and additive models. mgcv uses exactly the same idea (many basis functions, a roughness penalty, smoothing selected automatically by REML/GCV rather than a prior), so it is the direct frequentist counterpart. We fit a single smooth of $\log(\text{ozone})$ on temperature and then a three-term additive model.

options(repr.plot.width=12, repr.plot.height=4.3)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths())); suppressMessages(library(mgcv))
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; GREY<-"#718096"
d<-read.csv("airquality.csv"); d$logO<-log(d$Ozone)
cat("airquality:", nrow(d), "complete days\n")

airquality: 111 complete days

1. A single penalised spline, and the smoothing parameter

mgcv::gam(logO ~ s(Temp)) lays down many basis functions and penalises their roughness, choosing the smoothing strength by REML — the frequentist analogue of inferring $\tau^2$. Forcing the smoothing very low or very high shows the penalty at work: too little chases noise, too much straightens toward a line, REML lands between.

g<-gam(logO~s(Temp, k=25), data=d, method="REML")
xg<-seq(min(d$Temp),max(d$Temp),length=200); pg<-predict(g,data.frame(Temp=xg),se.fit=TRUE)
g_under<-gam(logO~s(Temp,k=25,sp=1e-5), data=d); g_over<-gam(logO~s(Temp,k=25,sp=1e5), data=d)
par(mar=c(4,4,3,1)); plot(d$Temp,d$logO,pch=19,col=GREY,cex=.6,xlab="temperature",ylab="log ozone",main="Penalised spline (mgcv): the smoothing parameter at work")
lines(xg, predict(g_under,data.frame(Temp=xg)), col=RED, lwd=1.5, lty=3)
lines(xg, predict(g_over,data.frame(Temp=xg)), col=GREEN, lwd=1.5, lty=2)
polygon(c(xg,rev(xg)),c(pg$fit-1.96*pg$se.fit,rev(pg$fit+1.96*pg$se.fit)),col=adjustcolor(BLUE,.2),border=NA); lines(xg,pg$fit,col=BLUE,lwd=2.5)
legend("topleft", c("tiny penalty (undersmooth)","huge penalty (oversmooth)","REML-selected"), col=c(RED,GREEN,BLUE), lwd=2, lty=c(3,2,1), bty="n", cex=.8)
cat(sprintf("REML picks effective df = %.1f, between the noisy and near-linear extremes -- and returns a confidence band.\n", sum(g$edf)-1))
cat("mgcv's penalised spline is exactly the frequentist twin of the random-walk-prior Bayesian P-spline (Python nb);\n")
cat("as the GP notebook showed, a penalised spline is also a Gaussian-process posterior mean.\n")

REML picks effective df = 2.1, between the noisy and near-linear extremes -- and returns a confidence band.

mgcv's penalised spline is exactly the frequentist twin of the random-walk-prior Bayesian P-spline (Python nb);

as the GP notebook showed, a penalised spline is also a Gaussian-process posterior mean.

2. An additive model — three smooth effects

$\log(\text{ozone})=\beta_0+f_1(\text{solar})+f_2(\text{wind})+f_3(\text{temp})$. mgcv fits all three penalised smooths jointly and returns each partial effect with a confidence band.

gam3<-gam(logO~s(Solar.R)+s(Wind)+s(Temp), data=d, method="REML")
cat("deviance explained:", round(summary(gam3)$dev.expl*100,1), "%   edf per term:", round(summary(gam3)$edf,1), "\n")
# Draw the partial effects from predict(type="terms") rather than plot.gam(shade=TRUE):
# on this mgcv/R version `shade` leaks into the base graphics calls and emits 27 identical
# '"shade" is not a graphical parameter' warnings into the output. This also lets us set ylim
# explicitly so the credible bands cannot run off the top of the panels.
vars <- c("Solar.R","Wind","Temp"); labs <- c("f(solar radiation)","f(wind)","f(temperature)")
par(mfrow=c(1,3), mar=c(4,4,3,1))
for (k in seq_along(vars)) {
  v  <- vars[k]
  xs <- seq(min(d[[v]]), max(d[[v]]), length=200)
  nd <- d[rep(1,200), c("Solar.R","Wind","Temp")]; nd[[v]] <- xs
  pt <- predict(gam3, nd, type="terms", se.fit=TRUE)
  cn <- paste0("s(", v, ")")
  fitk <- pt$fit[, cn]; sek <- pt$se.fit[, cn]
  ylim_k <- range(c(fitk - 1.96*sek, fitk + 1.96*sek))
  plot(xs, fitk, type="n", ylim=ylim_k, xlab=v, ylab="partial effect", main=labs[k])
  polygon(c(xs, rev(xs)), c(fitk - 1.96*sek, rev(fitk + 1.96*sek)),
          col=adjustcolor(BLUE,.25), border=NA)
  lines(xs, fitk, col=BLUE, lwd=2.3); rug(d[[v]], col=GREY)
}
par(mfrow=c(1,1))
cat("Temperature raises ozone strongly and nonlinearly, wind lowers it, solar radiation raises then flattens --\n")
cat("the same three partial effects the from-scratch backfitting Gibbs recovered in the Python notebook.\n")

deviance explained: 70.7 %   edf per term: 2.2 2.5 2 

Temperature raises ozone strongly and nonlinearly, wind lowers it, solar radiation raises then flattens --

the same three partial effects the from-scratch backfitting Gibbs recovered in the Python notebook.

3. Summary

mgcv — the standard R package for penalised splines and additive models — reproduces the Python results by penalised likelihood with REML-selected smoothing rather than a Bayesian random-walk prior. On the ozone data a single smooth of temperature landed between the under- and over-smoothed extremes with a confidence band, and a three-term GAM separated the nonlinear effects of solar radiation, wind and temperature — matching the from-scratch Gibbs and PyMC fits in pspline_python.ipynb.

This is the penalty paradigm, and mgcv is its industrial engine — the frequentist counterpart to the Bayesian P-spline. It complements the variable-selection arc's free-knot splines (adaptive knot locations by reversible-jump MCMC), and, as the Gaussian-process notebook showed, a penalised spline is a GP in disguise. Next the arc turns to nonparametric survival and hazards.