Bayesian Joint Modelling Approach

• The model's main contribution is integrating multiple submodels into a unified Bayesian framework that rigorously propagates uncertainty.
• The approach employs hierarchical structures and joint priors to mitigate spatial confounding and capture fine- versus coarse-scale effects.
• It facilitates effective inference in settings with interconnected outcomes and latent variables across spatial and statistical domains.

A Bayesian joint modelling approach integrates multiple stochastic and hierarchical processes within a unified probabilistic framework, enabling rigorous sharing and propagation of information—particularly uncertainty—across different components or submodels. Bayesian joint models are widely employed across spatial statistics, biostatistics, machine learning, and other domains to address settings where outcomes are interconnected or exhibit mutual dependence, and where fully accounting for latent variables or nuisance processes is required for valid inference.

1. Formal Specification: Structure and Priors

Central to Bayesian joint modelling is the explicit factorization of the observed data likelihood and the specification of joint priors for all latent quantities and parameters. In the context of spatial statistics, Marques & Wiemann (2023) proposed a Bayesian spatial+ framework that addresses spatial confounding by specifying a full hierarchical model for both a spatially-indexed covariate $x(\mathbf{s})$ and a spatial response $y(\mathbf{s})$. The model adopts a spline-based representation:

• Observations at locations $\mathbf{s}_1, \dots, \mathbf{s}_n$.
• Spline basis $Z = [U \, V]$ (with unpenalized and penalized components via eigendecomposition of the penalty matrix $K$).
• Covariate and response submodels:

$$x = \beta_0^x \mathbf{1}_n + V \gamma^{x}_\mathrm{pen} + U \gamma^{x}_\mathrm{unpen} + \epsilon^x, \quad y = \beta_0 \mathbf{1}_n + \beta \epsilon^x + V \gamma_\mathrm{pen} + U \gamma_\mathrm{unpen} + \epsilon.$$

Here, $\epsilon^x \sim N(0, \sigma_x^2 I_n)$ and $\epsilon \sim N(0, \sigma^2 I_n)$.

The Bayesian hierarchy includes:

• Flat priors on intercepts and regression coefficients,
• Penalized spline coefficients with normal priors parameterized by smoothness parameters $\tau^2$,
• Uniform or weakly-informative priors on variance parameters,
• Critically, a joint prior linking the smoothness parameters of the spatial effects: either $\tau_x = (\tau \sigma_x / \sigma) + \xi$ or $\tau_x^2 = (\tau^2 \sigma_x^2 / \sigma^2) + \xi$, with $\xi$ assigned a non-negative, light-tailed prior (e.g., Uniform or inverse-Gamma).

This joint prior enforces a smoothness ordering that prevents the response spatial effect from capturing finer-scale variability than the covariate effect, thereby mitigating spatial confounding at high frequencies (Marques et al., 2023).

2. Comparison to Two-Stage and Frequentist Methods

A defining technical feature of the Bayesian joint modelling paradigm is the full and automatic propagation of uncertainty across model components. The frequentist spatial+ approach (Dupont et al., 2022) estimates the covariate's spatial structure in a first stage, then plugs these estimates into a second-stage spatial regression. Such "plug-in" estimators do not propagate error from Stage 1 to Stage 2, and can reintroduce bias, especially if the high-frequency structure is not adequately separated.

In contrast, the