Efficient Bayesian Multivariate Surface Regression

Methods for choosing a fixed set of knot locations in additive spline models are fairly well established in the statistical literature. While most of these methods are in principle directly extendable to non-additive surface models, they are less likely to be successful in that setting because of the curse of dimensionality, especially when there are more than a couple of covariates. We propose a regression model for a multivariate Gaussian response that combines both additive splines and interactive splines, and a highly efficient MCMC algorithm that updates all the knot locations jointly. We use shrinkage priors to avoid overfitting with different estimated shrinkage factors for the additive and surface part of the model, and also different shrinkage parameters for the different response variables. This makes it possible for the model to adapt to varying degrees of nonlinearity in different parts of the data in a parsimonious way. Simulated data and an application to firm leverage data show that the approach is computationally efficient, and that allowing for freely estimated knot locations can offer a substantial improvement in out-of-sample predictive performance.

💡 Research Summary

The paper introduces a Bayesian framework for multivariate Gaussian regression that simultaneously captures linear effects, additive nonlinearities, and high‑dimensional interactions through spline bases. The model is expressed as
  Y = X₀B₀ + Xₐ(ξₐ)Bₐ + Xₛ(ξₛ)Bₛ + E,
where X₀ contains the original covariates (including an intercept), Xₐ(ξₐ) holds univariate spline bases for each covariate with scalar knots ξₐ, and Xₛ(ξₛ) contains multivariate radial‑basis (thin‑plate) spline functions for the interaction surface with vector‑valued knots ξₛ. The error matrix E follows an i.i.d. multivariate normal distribution with covariance Σ.

A central difficulty in spline regression is the selection of knot locations, especially for surface components where knots live in a q‑dimensional space and the curse of dimensionality makes a uniform coverage impossible. Traditional approaches fix knots (often at k‑means centroids) and rely on variable‑selection priors or reversible‑jump MCMC to prune them. Such methods struggle when the number of covariates exceeds a few, because moving a single knot at a time leads to slow mixing and poor exploration of the posterior.

To overcome this, the authors propose three interlocking innovations:

1. Model decomposition – By separating the regression into three distinct parts (linear, additive, surface), the additive component absorbs most of the marginal nonlinearity, allowing the surface component to be represented with a relatively small set of multivariate knots. This reduces the dimensionality of the knot space that must be explored.

2. Joint knot updating – The authors derive a marginal posterior for the knot locations after analytically integrating out the regression coefficients B. This marginal posterior is used to construct a Metropolis‑Hastings proposal that updates all knots simultaneously. The proposal distribution is tuned by the gradient of the log‑posterior with respect to the knots, which is computed efficiently by exploiting sparsity in the design matrices. Consequently, the algorithm achieves rapid mixing even in high‑dimensional covariate spaces.

3. Separate shrinkage priors – For each of the three coefficient blocks (B₀, Bₐ, Bₛ) and for each response variable, a Gaussian prior with a precision matrix λ P is placed, where P is either the identity or an g‑prior (XᵀX). The hyper‑parameters λ are treated as unknown and assigned log‑normal priors, enabling the data to dictate the amount of shrinkage separately for linear, additive, and surface parts. This hierarchical shrinkage mitigates over‑fitting while preserving flexibility where the data demand it.

The number of knots is not sampled within the MCMC; instead, it is selected beforehand by Bayesian cross‑validation using the log‑predictive score, which can be parallelised across candidate knot counts. Although this does not incorporate knot‑count uncertainty as RJMCMC does, it yields a robust and computationally efficient procedure.

The methodology is evaluated on two fronts. In simulated experiments, data are generated with varying degrees of nonlinearity and interaction complexity. The “free‑knot” model consistently outperforms a comparable model with fixed knots in terms of out‑of‑sample predictive error. In a real‑world application to firm leverage (two‑response) data, the same pattern holds: allowing knots to move freely improves predictive performance and produces sensible heat‑maps of knot density that highlight regions of the covariate space where the surface is most active.

Computationally, the joint‑knot MCMC runs for several thousand iterations in a matter of minutes on a standard workstation, demonstrating that the analytical integration of B and the gradient‑based proposal dramatically reduce the per‑iteration cost.

In summary, the paper delivers a coherent solution to the knot‑selection problem in multivariate spline regression by (i) decomposing the model to keep the surface knot set small, (ii) integrating out regression coefficients to obtain a tractable marginal posterior for knots, (iii) using gradient‑adapted joint proposals for efficient exploration, and (iv) employing hierarchical shrinkage priors that adapt to each model component and each response. The resulting algorithm is both statistically powerful—yielding superior predictive accuracy—and computationally scalable, making it a valuable addition to the toolbox of Bayesian non‑parametric regression.