Automatic smoothness selection for arbitrary statistical models with penalised splines, simple random effects, or tensorproducts

Efficiently fits models involving quadratic penalties of the form $$\sum_{i} \lambda_i b^T S_i b,$$ by combining automatic differentiation via RTMB with a custom implementation of the extended Fellner-Schall update.

Users only need to supply a custom penalised log-likelihood function that calls penalty() or penalty2().

Usage

qreml(
  pnll,
  par,
  dat,
  random,
  map = NULL,
  silent = 1,
  spname = NULL,
  tol_edf = 0.01,
  maxiter = 100,
  alpha = 0.1,
  tol = 0.01,
  smoothing = 1,
  method = "BFGS",
  lsp_max = 15,
  lsp_min = -15,
  step_small = 0.05,
  max_halve = 6,
  control = list(),
  spHess = FALSE,
  joint_unc = FALSE,
  saveall = FALSE
)

Arguments

pnll

penalised negative log-likelihood function that is structured as dictated by RTMB and uses the penalty or penalty2 function to compute the penalty

par

named list of initial parameters

dat

initial data list that contains the data used in the likelihood function, hyperparameters, and the initial penalty strength vector, called sp or lambda (see spname)

random

vector of names of the random effects/ penalised parameters in par

Caution: The ordering of random needs to match the order of the random effects passed to penalty.

map

optional map argument, containing factor vectors to indicate parameter sharing or fixing

silent

integer silencing level: 0 corresponds to full printing of inner and outer iterations, 1 to printing of outer iterations only, and 2 to no printing

spname

optional name of the penalty strength vector in dat.

By default, dat$sp is used if present, otherwise dat$lambda. Hence, if dat contains sp, lambda is free to be used for something else, e.g. the mean of a Poisson distribution. In the likelihood function, the penalty strengths are accessed by the same name, e.g. penalty(beta, S, sp), as qreml() replaces them with the current values in each iteration. The returned model object stores them under this name as well, e.g. mod$sp.

tol_edf

convergence tolerance on the effective degrees of freedom.

This is the primary convergence criterion. The iteration stops once no smooth's effective degrees of freedom has changed by more than tol_edf across the last four outer iterations, provided the step in log(lambda) is also small. The effective degrees of freedom are what say whether the fitted smooth is still changing, and unlike a tolerance on the restricted likelihood they mean the same thing on every model, being measured in effective parameters rather than in nats. The penalty strengths themselves are a poor proxy: lambda can slide along a flat ridge for many iterations, moving substantially in relative terms, while neither the criterion nor the fit changes appreciably.

maxiter

maximum number of outer iterations

alpha

smallest factor by which a penalty strength may decrease in one outer iteration, a number in [0, 1). Defaults to 0.1.

Penalty strengths are free to increase as fast as the update proposes, but cannot collapse faster than this per iteration. Reducing a penalty strength too quickly can push the inner optimisation into a local optimum or a numerically awkward region, which matters more here than in a GAM because the likelihood is user-written. Set to zero to remove the floor entirely. Step length is handled separately, by max_halve and the step multiplier.

tol

fallback convergence tolerance, on the restricted log-likelihood. Defaults to 0.01.

Used only when the effective degrees of freedom are not trustworthy, i.e.\ when they fall outside $[0, K_i]$, which happens when the data Hessian is indefinite (see tol_edf and the returned edf_valid). The iteration then stops once the restricted log-likelihood has changed by less than tol over the last four outer iterations. A tolerance in nats is not comparable across models, which is why it is the fallback rather than the primary criterion.

smoothing

optional scaling factor for the final penalty strength parameters. Increasing this beyond one leads to a smoother final model

method

inner optimisation method to be used by optim. Defaults to "BFGS"

lsp_max

largest value allowed for log(lambda). Defaults to 15, as in mgcv, i.e. penalty strengths saturate at roughly 3.3e6

lsp_min

smallest value allowed for log(lambda). Defaults to -15, i.e. penalty strengths stop at roughly 3e-7. A penalty strength on a smooth the data cannot inform otherwise keeps falling towards zero along a flat restricted likelihood, and its step never becomes small enough for convergence. A penalty strength ending at this bound gives a warning, as the penalty is then effectively switched off. If no other penalty acts on the same coefficients, the smooth is effectively unpenalised, which usually means its basis dimension is too small. If another penalty still does, e.g. the ridge of the random effect margin of a tensor product with bs = "re", the penalty at the bound is merely redundant: this is usually harmless, and a larger basis dimension does not help.

step_small

size of a step in log(lambda) below which convergence may be declared. Defaults to 0.05, as in mgcv

max_halve

maximum number of times a step that decreases the restricted likelihood is halved before it is accepted anyway. Defaults to 6.

mgcv never shortens below the full Fellner-Schall step and accepts a worse one instead, which for a user-written likelihood can drift downhill for tens of iterations.

control

list of control parameters for optim to use in the inner optimisation

spHess

logical, if TRUE, the sparse automatic differentiation Hessian is used for evaluation. The factorisation is dense either way

joint_unc

logical, if TRUE, joint RTMB object is returned allowing for joint uncertainty quantification

saveall

logical, if TRUE, then all model objects from each iteration are saved in the final model object

Value

model object of class "qremlModel", inheriting from "LaMaModel" (see lama), carrying the fitted quantities, the smoothness selection diagnostics shown by summary.qremlModel, outer_hessian(), the outer Hessian that sdreport_outer turns into standard errors, and joint_uncertainty(), the joint precision (or corrected covariance) accounting for the uncertainty in the penalty strengths

Details

Step size control follows the extended Fellner-Schall implementation in mgcv while convergence is judged on the effective degrees of freedom (differing from mgcv).

The criterion. The penalty strengths are chosen to maximise the Laplace-approximate restricted log-likelihood $$V(\lambda) = \ell_p(\hat{b}; \lambda) + \tfrac{1}{2} \log |S_\lambda|_+ - \tfrac{1}{2} \log |J|,$$ where $\ell_p$ is the penalised log-likelihood at its mode $\hat{b}$, $S_\lambda = \sum_i \lambda_i S_i$, $|\cdot|_+$ denotes the product of the positive eigenvalues, and $J$ is the Hessian of the penalised negative log-likelihood at $\hat{b}$. Every outer iteration refits the inner problem to obtain $\hat{b}$ and $J$, so one call is a sequence of complete model fits, each warm started at the previous one.

The update. The penalty strengths are updated multiplicatively, $\lambda_i \leftarrow \lambda_i r_i$, by the extended Fellner-Schall ratio $$r_i = \frac{\mathrm{tr}(S_\lambda^- S_i) - \mathrm{tr}(J^{-1} S_i)}{\hat{b}^T S_i \hat{b}},$$ which equals one exactly where $\partial V / \partial \lambda_i = 0$ and otherwise moves $\lambda_i$ in the direction of the gradient, without ever making it negative. For a smooth carrying a single penalty the first trace is $\mathrm{rank}(S_i) / \lambda_i$, obtained once from an eigendecomposition instead of at every iteration. Only the diagonal blocks of $J^{-1}$ enter, and these come from the Cholesky factor by triangular solves rather than from a full inverse.

Step control. The step is taken in $\log \lambda$ and scaled by a multiplier that persists across iterations: it is halved, up to max_halve times, while the criterion does not improve, and doubles back towards one after each improving step. Unlike mgcv, there is no acceleration beyond the full step: each trial at a larger step costs a complete inner fit, which for a user-written likelihood is expensive, and rarely pays off. A step that still worsens the criterion after that many halvings is accepted, which is the only way the iteration can move downhill. Changes in the criterion smaller than $10^{-6} (1 + |V|)$ count as noise here, so that a step of essentially zero at the optimum does not set off halvings. In addition no penalty strength may fall by more than a factor alpha in one iteration, and $\log \lambda$ is kept within [lsp_min, lsp_max]. The fit returned is the best iterate seen, not necessarily the last one.

Convergence. The iteration stops once no smooth's effective degrees of freedom $$\mathrm{edf}_i = K_i - \mathrm{tr}(J^{-1} \lambda_i S_i)$$ has moved by more than tol_edf over the last four iterations, provided the step in $\log \lambda$ is also small; see tol_edf. It also stops once the inner optimiser no longer moves from its warm start (after the first three iterations and at a full step), as the change in the penalty strengths then no longer affects the fit; this ends a penalty strength creeping towards zero along a flat restricted likelihood, whose step never becomes small. Being eigenvalue sums, the effective degrees of freedom must lie in $[0, K_i]$, and they leave that range only if the data Hessian is indefinite, in which case the traces driving the update are not the quantities the method assumes. That is reported rather than repaired, through the returned edf_valid and a warning, and the criterion then falls back to tol on the restricted log-likelihood.

Numerics. $J$ is the exact automatic differentiation Hessian, which matters because the finite differencing error of optimHess is of the same order as the criterion differences the convergence test has to resolve. It is factorised by a pivoted Cholesky decomposition after scaling to a unit diagonal; should it not be positive definite, a ridge is added and reported through the returned hessian_repaired.

Uncertainty. The returned outer_hessian() gives the Hessian of the negative restricted log-likelihood in the log penalty strengths, by central differences of its approximate gradient, which sdreport_outer turns into standard errors. Note that these describe the curvature of $V$ only, and treat the approximation to it as exact. The returned joint_uncertainty() propagates this uncertainty to the parameters by the first-order correction of Wood, Pya and Saefken (2016): it gives the joint precision matrix of the parameters and the log penalty strengths, as sdreport(getJointPrecision = TRUE) does in TMB, which stays sparse if the Hessian is (spHess = TRUE), or, with precision = FALSE, the corrected covariance matrix of the parameters alone. Both take the outer Hessian as argument H_rho, by default computed by outer_hessian().

References

Wood, S. N., & Fasiolo, M. (2017). A generalized Fellner-Schall method for smoothing parameter optimization with application to Tweedie location, scale and shape models. Biometrics, 73(4), 1071-1081.

Koslik, J. O. (2024). Efficient smoothness selection for nonparametric Markov-switching models via quasi restricted maximum likelihood. arXiv preprint arXiv:2411.11498.

Examples

data = trex[1:2000,] # subset of the data

# initial parameter list
par = list(logmu = log(c(0.3, 2.5)), # step mean
           logsigma = log(c(0.3, 1.5)), # step sd
           beta0 = c(-2,-2), # state process intercept
           beta_spline = matrix(rep(0, 18), nrow = 2)) # state process spline coefs
          
# data object with initial penalty strength sp
dat = list(step = data$step, # step length
           tod = data$tod, # time of day covariate
           N = 2, # number of states
           sp = rep(20,2)) # initial penalty strength

# building model matrices
modmat = modmat(~ s(tod, bs = "cp"), 
                       data = data.frame(tod = 1:24), 
                       knots = list(tod = c(0,24))) # wrapping points
dat$Z = modmat$Z # spline design matrix
dat$S = modmat$S # penalty matrix

# penalised negative log-likelihood function
pnll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm(cbind(beta0, beta_spline), Z) # transition probabilities
  delta = stationary_p(Gamma, t = tod[1] %% 24 + 1) # initial distribution
  mu = exp(logmu) # step mean
  sigma = exp(logsigma) # step sd
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step)) # only for non-NA obs.
  for(j in 1:N) allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])
  -forward(delta, Gamma[,,tod], allprobs) +
      penalty(beta_spline, S, sp) # this does all the penalization work
}

# model fitting
mod = qreml(pnll, par, dat, random = "beta_spline", silent = 2)