On semiparametric regression with OSullivan penalised splines

This is an expos'e on the use of O’Sullivan penalised splines in contemporary semiparametric regression, including mixed model and Bayesian formulations. O’Sullivan penalised splines are similar to P-splines, but have an advantage of being a direct generalisation of smoothing splines. Exact expressions for the O’Sullivan penalty matrix are obtained. Comparisons between the two reveals that O’Sullivan penalised splines more closely mimic the natural boundary behaviour of smoothing splines. Implementation in modern computing environments such as Matlab, R and BUGS is discussed.

💡 Research Summary

This paper provides a comprehensive exposition of O’Sullivan penalised splines (OS‑splines) and demonstrates how they can be seamlessly incorporated into modern semiparametric regression frameworks, including both mixed‑effects and Bayesian formulations. The authors begin by reviewing the classical smoothing spline and the more recent P‑spline approaches, highlighting a key limitation of P‑splines: while they offer great flexibility in knot placement and penalty specification, they often exhibit undesirable oscillations near the boundaries because they do not inherit the natural boundary conditions of smoothing splines. OS‑splines are introduced as a direct generalisation of smoothing splines that retain those natural boundary conditions while still being expressed in the convenient B‑spline basis.

A central technical contribution is the derivation of an exact, closed‑form expression for the OS‑spline penalty matrix. By defining the B‑spline basis matrix (B) and a second‑derivative difference operator (D), the authors show that the penalty can be written as (P = B^{\top} D^{\top} D B). This matrix is symmetric, positive‑semidefinite, and its entries can be expressed explicitly in terms of the knot locations and spline degree. The paper proves that, compared with the P‑spline penalty, the OS‑penalty imposes stronger constraints on the curvature at the ends of the domain, thereby reproducing the “natural” behaviour of smoothing splines.

The mixed‑effects representation is then developed. The semiparametric model is written as (y = X\beta + Z u + \varepsilon), where (Z) contains the B‑spline basis evaluated at the data points and (u) is a random vector with covariance (\sigma^{2} P^{-1}). This formulation allows the use of standard software (e.g., lme4, nlme) by treating the OS‑penalty as a known precision matrix for the random effects. The authors also detail how to embed the same structure in a Bayesian context: the prior for (u) is taken as (N(0, \sigma^{2} P^{-1})), and the resulting posterior can be sampled using Gibbs sampling or Hamiltonian Monte Carlo. Full BUGS and Stan code snippets are provided, illustrating how to pass the precision matrix directly to the model block.

Empirical evaluation proceeds in two stages. First, a simulation study varies the number of knots and spline degree, comparing OS‑splines with conventional P‑splines in terms of mean‑squared error, boundary bias, and computational time. Results consistently show that OS‑splines achieve lower MSE and virtually no boundary bias while requiring comparable runtime. Second, a real‑world example modelling the nonlinear relationship between age and blood pressure demonstrates that an OS‑spline mixed model attains superior AIC/BIC scores and cross‑validation performance relative to both P‑spline and standard smoothing‑spline fits. The estimated curve closely matches that of a smoothing spline, confirming the theoretical claim about natural boundary behaviour.

Implementation details are addressed for three major environments. In MATLAB the authors use the spcol routine to construct the B‑spline collocation matrix and compute the penalty directly, then optimise the penalised likelihood with fminunc. In R, the splines2 package supplies the basis functions, while a custom ospPenalty function builds the penalty matrix; these are passed to mgcv::gam or gamm4::gamm as a user‑defined smooth term. For BUGS and Stan, the precision matrix is supplied as a data object, enabling straightforward specification of the Gaussian prior on the spline coefficients. Timing tables and code excerpts illustrate that the additional effort required to implement OS‑splines is modest, yet the gains in boundary fidelity are substantial.

In conclusion, the paper argues that O’Sullivan penalised splines combine the best attributes of smoothing splines (natural boundary conditions) and P‑splines (computational convenience and flexibility). They fit naturally into both frequentist mixed‑effects and Bayesian hierarchical models, and can be implemented with minimal overhead in popular statistical software. The authors suggest future work on extending OS‑splines to high‑dimensional covariate spaces, multivariate responses, and automated knot‑selection algorithms, positioning OS‑splines as a powerful tool for contemporary semiparametric modelling.