Stasinopoulos 2007 - Bayesian P-Splines and Distributional Regression in PyMC

Summary

GAMLSS’s penalized likelihood ${\ell }_p$ is the log-posterior of a model with Gaussian random-walk priors on spline coefficients. In PyMC you give each distribution parameter its own P-spline with a second-order random-walk (RW2) prior, a log link for positive scales, and a hyperprior on the smoothing variance. This note is a concept translation (not in the 2007 paper); the PyMC code is untested pseudo-real code.

Overview

Equation (6) of the paper, ${\ell }_p=\ell -\frac{1}{2}\sum {\lambda }_{jk}{\mathbf{\gamma }}_{jk}^{\top }{\mathbf{G}}_{jk}{\mathbf{\gamma }}_{jk}$, equals $\log p(\mathbf{y}\mid \mathbf{\gamma })+\log p(\mathbf{\gamma }\mid \lambda )$ up to constants with prior $\mathbf{\gamma }\sim N(0,(\lambda \mathbf{K}{)}^{-})$. GAMLSS therefore returns the posterior mode for fixed $\lambda$; a Bayesian analysis samples $(\mathbf{\gamma },\lambda )$ jointly and propagates smoothing uncertainty into the centile curves.

Main Content

Bayesian P-spline (RW2)

For parameter $k$, with B-spline basis $\mathbf{B}(x)\in {\mathbb{R}}^{n\times q}$ and difference matrix ${\mathbf{D}}_2$ ($(q-2)\times q$),

$${\eta }_k(x)={\beta }_{0k}+\mathbf{B}(x){\mathbf{\gamma }}_k,{\mathbf{D}}_2{\mathbf{\gamma }}_k\sim N(0,{\tau }_k^2\mathbf{I}),{\tau }_k\sim \text{HalfNormal}.$$

The smoothing parameter is ${\lambda }_k={\tau }_k^{-2}$; the null space of ${\mathbf{D}}_2^{\prime }{\mathbf{D}}_2$ (constant and linear) is handled by the intercept and a free linear slope (“centred RW2”). Equivalent GAMLSS ps(x, order=2).

Growth-chart model for E41 (PyMC 5 sketch, untested)

import numpy as np, pymc as pm, pytensor.tensor as pt
from patsy import dmatrix
 
# x = age (or age**xi), y = BMI, B = (n, q) cubic B-spline design matrix
B = np.asarray(dmatrix("bs(x, df=20, degree=3, include_intercept=False) - 1",
                       {"x": x}))
q = B.shape[1]
 
def rw2_spline(name, scale_sd=1.0):
    tau = pm.HalfNormal(f"tau_{name}", scale_sd)
    z = pm.Normal(f"z_{name}", 0, 1, shape=q)            # non-centred increments
    g = pt.cumsum(pt.cumsum(z * tau)) / q                # RW2 via double cumsum
    b0 = pm.Normal(f"b0_{name}", 0, 2)
    b1 = pm.Normal(f"b1_{name}", 0, 1)                   # linear part
    return b0 + b1 * x_std + pt.dot(B, g)
 
with pm.Model() as m:
    mu = rw2_spline("mu")                    # identity link (centred on data scale)
    sigma = pt.exp(rw2_spline("sigma", 0.5)) # log link
    nu = rw2_spline("nu", 0.5)               # Box-Cox power, identity
    # LMS: Z = ((y/M)**nu - 1)/(nu*S) ~ Normal(0,1), with Jacobian y**(nu-1)/(M**nu * S)
    M = pt.exp(mu)
    z = ((y / M) ** nu - 1) / (nu * sigma)
    pm.Potential("lms", pt.sum(-0.5 * z**2 + (nu - 1) * pt.log(y) - nu * pt.log(M) - pt.log(sigma)))
    idata = pm.sample(1000, tune=1000, target_accept=0.95)

Notes: (i) replace Normal by Student-$t$ for BCT (add a tau = 2 + HalfNormal-type df); (ii) the Jacobian term is the Box-Cox change of variables; (iii) restrict $|\nu \sigma z|<1$ or the support is violated — use a centred parametrisation and tight priors on $\nu$; (iv) E41 compares constant spread, smooth spread, skewed and heavy-tailed variants by LOO.

Comparing variants. Fit (a) $\sigma$ constant, (b) smooth $\sigma$, (c) add $\nu$ (skew), (d) add $\tau$ (heavy tails) and compare by PSIS-LOO or by held-out log score; the Bayesian counterpart of GAIC/VGD (Stasinopoulos 2007 - Fitting Algorithms and Model Selection). Check calibration by PIT / normalised quantile residuals binned by age (worm-plot analogue; Posterior Predictive Checking).

Pitfalls. BCCG/BCT likelihood has a data-dependent support when $\nu \ne 0$; the paper warns about likelihood singularities for small samples (penalise $\sigma \to 0$). With many P-spline blocks use non-centred increments and target_accept $\ge 0.9$.

Examples

Alternative smoother: a Hilbert Space Gaussian Processes prior on each ${\eta }_k(x)$ gives the same distributional structure with kernel-based smoothness and fewer knot-placement decisions.

Connections

Stasinopoulos 2007 - Centile Curves, LMS and Residual Diagnostics, Hierarchical Models (random-effect terms in $\mathbf{Z}\mathbf{\gamma }$ as partial pooling), Prior Predictive Checking for the smoothing priors.

See Also

Quantile Regression, Generalized Linear Models