<figure><img src="image_1.png" /></figure>

# Addressing spatial confounding in geostatistical regression models: An R-INLA approach

Jérémy Lamouroux<fnref n="1" /> ID | Alizée Geffroy<fnref n="1" /> | Sébastien Leblond<fnref n="2" /> | Caroline Meyer<fnref n="2" /> |
Isabelle Albert<fnref n="1" />

## Abstract

1. Spatial confounding, which has been studied extensively in recent years, can explain inconsistencies between results obtained by regression models with and without spatial modelling. While the most common solutions for spatial confounding mitigation often target areal data or geostatistical data modelling by splines, we aim to extend some resolution methods for geostatistical data modelling by Gaussian Markov Random Fields (GMRF) using R-INLA methodology, aided by the R-inlabru package.

2. We focus on four sought-after approaches for alleviating spatial confounding: Restricted spatial regression (RSR), geoadditive structural equation model (gSEM), Spatial+, and its simplified version, called here Spatial+2.0. We show how each can be implemented from geostatistical data in a GMRF framework using R-inlabru.

3. Next, a simulation study that produces a spatial confounding phenomenon in a geostatistical context is carried out to assess the coherence of the proposed extensions with the expectations of these methods. Then we apply the expanded approaches to a case study, linking cadmium (Cd) concentration in terrestrial mosses to Cd concentration in air.

4. Our findings support the feasibility of our extended approach of spatial confounding resolution methods to geostatistical data using R-inlabru in keeping with the previous contexts, although certain precautions and limitations must be considered.

## KEYWORDS

cadmium, environmental pollution, GMRF, gSEM, inlabru, restricted spatial regression, spatial+, terrestrial mosses

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹Université Paris-Saclay, AgroParisTech, INRAE, UMR MIA Paris-Saclay, Palaiseau, France</footnote>

<footnote n="2">²PatriNat (OFB-MNHN), Paris, France</footnote>

<footnote>**Correspondence**
Jérémy Lamouroux
Email: jeremy.lamouroux@hotmail.fr</footnote>

<footnote>Handling Editor: Pietro Milanesi</footnote>

<footnote>This is an open access article under the terms of the [Creative Commons Attribution-NonCommercial-NoDerivs License](http://creativecommons.org/licenses/by-nc-nd/4.0/), which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.
© 2025 The Author(s). Methods in Ecology and Evolution published by John Wiley & Sons Ltd on behalf of British Ecological Society.</footnote>
<figure><img src="image_3.png" /></figure>

# 1 | INTRODUCTION

Spatial confounding, which has garnered particular interest in recent years (e.g. in 2024 and 2025, Gilbert et al., 2025; Hui et al., 2024; Urdangarin et al., 2024), is evoked when very different results are obtained on the effects of the covariates in a spatial model compared with those of the same model without spatial information. Also, the prediction error can be reduced in the presence of spatial confounding (Page et al., 2017) affecting the model's predictive performance and spatial interpolation.

## 1.1 | Overview of spatial confounding problem

Clayton et al. (1993) describe the possibility of this phenomenon very early, stating that when ‘the pattern of variation of the covariate is similar to the disease risk, the location may act as a confounder’. Gilbert et al. (2024) resume four possible phenomena referred to spatial confounding in literature. ‘(i) Omitted confounder bias: the existence of an unmeasured confounding variable with a spatial structure that influences both the exposure and the outcome; (ii) regularisation bias: the finite-sample bias of methods that use flexible regression functions like splines or Gaussian processes to control for an unknown function of space; (iii) random effect collinearity: the change in fixed effect estimates that may come about when spatially-dependent random effects added to a regression model are collinear with the covariates; (iv) concurvity: the difficulty of assessing the effect of an exposure, which is or is close to a smooth function of space if an arbitrary smooth function of space is also included in a regression model’. The omitted confounder bias can emerge in any statistical model if it is not in the causal inference framework. The regularisation bias is inherent to statistics which work with finite samples, spatial statistics being more particularly affected as they work with complex modelling such as splines or Gaussian fields. The random effect collinearity and concurvity are of particular interest to us because they describe observed superposed spatial variabilities, which may be due to spatial collinearity or ‘concurvity’ (i.e. a space-varying co-smoothness) between covariates and the outcome variable, and, modelling and separating this superimposition in spatial non-causal models can be an objective for producing relevant conclusions from finite-sample case studies.

## 1.2 | Methods to mitigate spatial confounding

Many authors position their work within specific contexts, focusing on particular sources of spatial confounding as the ones listed above. They aim to develop methods that address these complexities while acknowledging the uncertainty that surrounds the phenomenon in reality. In this article, we explore four methods that have been proposed to address more particularly the last two issues of spatial confounding listed above. We first look at a widely used restricted

spatial regression (RSR) method proposed by Reich et al. (2006) and developed later by Hanks et al. (2015). This method addresses collinearity between covariates and the spatial random effects introduced into the regression model while preserving the fixed effect estimates obtained without spatial random effects. We then investigate the geoadditive structural equation model (gSEM), developed by Thaden and Kneib (2018), in their non-unified version including three successive regression steps to address spatial confounding by removing spatial dependence from both the response and the covariate before estimating their relationship. Thereafter we examine the spatial+ method proposed by Dupont et al. (2022), which is a simpler two-step regression method developed for models using splines to model spatial effects. Lastly, we explore the recent simplified spatial+ approach, called here spatial+2.0, developed by Urdangarin et al. (2024), which combines the two regression steps into a single regression and then addresses the residuals' uncertainty from the first regression in spatial+.

These four well-received resolution methods have been developed to handle data types such as areal data or geostatistical data by using partial thin plate spline models to approximate Gaussian random fields (GRF). Areal data describe situations when space is divided into regions associated with a single aggregated value, such as a mean. This data type is very valuable in epidemiology for disease mapping to highlight regions with high risk and as a first step to discovering potential risk factors related to the response of interest. Many formulations have been proposed to capture unobserved spatial variability as the conditional autoregressive (CAR) prior (Besag, 1974) or the BYM model (Besag et al., 1991; Riebler et al., 2016). The RSR model was first proposed for this type of data (Hanks et al., 2015; Reich et al., 2006).

Another type of data to which these resolution methods have been applied is geostatistical data, where $Y(s)$ is a random outcome at a specific location and the spatial index $s$ can vary continuously in the fixed subset of $\mathcal{R}^d$. The location $s$ is typically a two-dimensional vector with latitude and longitude but may also include altitude (Blangiardo & Cameletti, 2015). This data type is very valuable in spatial modelling for interpolation mapping to highlight geographic variations of an outcome and as a first step to discover potential covariates influencing the response of interest. Tobler's First Law of Geography (Tobler, 1970), which states that 'everything is related to everything else, but near things are more related than distant things' is fundamental to the spatial approaches of geostatistical data. This law highlights the spatial autocorrelation often present in geostatistical data, where proximity is critical in the relationships between data points. In such cases, spline models offer a flexible way to capture the space-varying relationships between covariates and the outcome variable in a regression model. Splines are smooth and flexible functions employed to model non-linear relationships in data, particularly to capture continuous variations in space or time (Wahba, 1990). Spline models are used to capture continuous spatial variability from geostatistical data in the spatial+ method. Its extension spatial+2.0 is derived from spatial+ methodology with thin plate splines but it only has been applied to the ICAR model for areal data.
## 1.3 Integrated nested Laplace approximation methodology

To our knowledge, no study has yet explored applications of these four attractive spatial confounding resolution methods to geostatistical data modelled by Gaussian Markov Random Fields (GMRF) in an R-INLA (Integrated Nested Laplace Approximation) methodology (Lindgren et al., 2011), which is currently the most flexible toolbox for analysing data in continuous space. Thus, we propose extensions of these four resolution methods to mitigate spatial confounding from geostatistical data using the R-inlabru package (https://inlabru-org.github.io/inlabru/) (Bachl et al., 2019), the extension of the R-INLA package (http://www.r-inla.org), that addresses Bayesian spatial modelling from ecological survey data. The INLA methodology focuses on models that can be expressed as latent GMRF. Specifying the spatial regression models as Bayesian hierarchical models using latent Gaussian variables, Martinez-Beneito (2013) proposed a unified formulation for the different specifications of the spatial dependence between latent variables, all falling under the umbrella of multivariate GMRF. A GMRF is a GRF represented by an undirected graph and defined through its precision matrix (inverse of the covariance matrix). The key advantage of a GMRF is its efficiency, which comes from the sparsity of the precision matrix. The (i, j)-th entry of the precision matrix is non-zero only if nodes i and j are connected in the graph, meaning they are neighbours (Rue & Held, 2005). Using this sparsity, R-INLA allows efficient estimations of complex spatial models (Lindgren & Rue, 2015). It has also demonstrated high accuracy for these estimations while being computationally much faster than Markov Chain Monte Carlo (MCMC) methods (Bakka et al., 2018), due to Laplace approximation of the posterior parameters' distribution (Rue et al., 2009). With the Stochastic Partial Differential Equation approach (SPDE) (Lindgren et al., 2011), R-INLA accommodates the spatial continuous space and geostatistical data.

## 1.4 The case study

The extensions we proposed are motivated by a case study about the potential link between cadmium (Cd) concentrations in territorial mosses collected in France by the Biosurveillance des Retombées Atmosphériques Métalliques par les Mousses (BRAMM)<fnref n="1" /> and concentrations of this same element in the air, as predicted by the European system EMEP.<fnref n="2" /> These two concentrations in 2016 comprise this study's dataset on French territory. Several studies show the link between the two concentrations with other datasets (Bates, 1992; Travnikov & Ilyin, 2005; Schröder & Nickel, 2019). Indeed, mosses' unique morphology, absence of roots and developed vascular system allow them to function as effective sensors of atmospheric contaminants, as they rely directly on atmospheric deposition for nutrients (Bates, 1992). From our dataset, a linear regression model with a

random spatial effect does not find the link between the Cd concentration in mosses and the Cd concentration in the air predicted by the EMEP model, but the same model with no spatial random effect estimates a highly significant relationship between both concentrations. We have found ourselves in the same context as the examples of Reich et al. (2006), cited by Dupont et al. (2022), evoking the possibility of spatial confounding in this link problem. As the data are geostatistical, we have decided to extend the four previously mentioned resolution methods in this context using the R-INLA tool with R-inlabru.

## 1.5 Objectives and paper structure

The objective of this study is to implement and evaluate four methods for mitigating spatial confounding from geostatistical data within the R-INLA framework. By focusing on practical implementation, we aim to provide researchers and practitioners with effective guidelines to address this challenge within a widely used Bayesian modelling framework. We extend four notable methodological approaches and assess our extensions through a simulation study and a real-world case study, evaluating the covariate's effect and the model predictive performance with Bayesian goodness-of-fit criteria and spatial interpolation in the case study. This paper is structured as follows: Section 2 presents the four approaches we have studied and describes the adaptations developed for their implementation in R-inlabru for geostatistical data. Section 3 is dedicated to the simulation study, covering the data, the fitted models and results to compare the methods in a GMRF/geostatistical data context. This study allows us to examine the expected properties of these resolution methods in the proposed extensions. Section 4 focuses on the case study of Cd concentrations in terrestrial mosses, following a similar structure. Section 5 discusses all the results, highlighting the precautions and limitations needed to use these methods for geostatistical data processed by GMRF using the R-inlabru.

# 2 | MATERIALS AND METHODS

In the following, for the sake of simplicity, the methods and models are presented using a vector X. This vector contains a single covariate corresponding to our simulation and case studies, but note that these methods and models have been developed to handle a matrix X accommodating multiple covariates.

## 2.1 RSR

### 2.1.1 | Model description

It is appropriate to speak of RSR methods in the plural because of the variety of possible approaches (Hodges & Reich, 2010; Hughes & Haran, 2013; Paciorek, 2010; Reich et al., 2006), but all encapsulated by a general regression (Hanks et al., 2015):

<footnote n="1">¹https://bramm.mnhn.fr/.</footnote>

<footnote n="2">²https://www.emep.int/.</footnote>
$$ 
\begin{aligned}
Y(s) &= \beta_{\text{RSR}}X(s) + \tilde{U}(s) + \varepsilon(s), \\
\varepsilon(s) &\stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma_{\varepsilon}^2 I), \quad \tilde{U}(s) = (I - P_X)U(s), \quad U(s) \sim \mathcal{N}(0, \Sigma),
\end{aligned} 
\quad (1) $$

where $Y(s)$ is the outcome variable, where $s$ represents the data location $s = (s_1, ..., s_n)$, $X(s)$ is an observed covariate with unknown fixed effect $\beta_{RSR}$, $\varepsilon(s)$ is an independent and identically distributed errors' noise with standard deviation $\sigma_{\varepsilon}$, $I$ the identity matrix, $P_X = X(X^tX)^{-1}X^t$ projects onto the space spanned by the columns of $X = (X(s_1), ..., X(s_n))^t$, $U(s)$ is a spatial random effect with mean zero and a covariance matrix $\Sigma$, then $\tilde{U}(s)$ is a spatial random effect orthogonal to $X$ (Hanks et al., 2015). The method consists of fitting the regression (1) rather than a spatial regression with $U(s)$. The idea is to restrict the spatial effects to the space orthogonal to the fixed effect to prevent collinearity between the two, which is seen as the cause of spatial confounding in this model type.

Two important properties of the RSR model are (Hodges & Reich, 2010; Khan & Calder, 2022): (i) the posterior mean of the estimator of $\beta_{RSR}$ is equal to that of the null model (simple regression model with the covariate and without spatial random effect):

$$ E[\beta_{\text{RSR}}|Y] = E[\beta_{\text{Null}}|Y]; \quad (2) $$

(ii) the posterior variance of the estimator of $\beta_{RSR}$ is lower or equal to the posterior variance of the estimator $\beta_{Spatial}$ of the spatial model, including the covariate and the random spatial effect:

$$ \mathrm{Var}(\beta_{\text{RSR}}|Y) \le \mathrm{Var}(\beta_{\text{Spatial}}|Y). \quad (3) $$

Khan and Calder (2022) discuss in depth these two important properties of RSR methods as well as others showing, through their results, that some are not achieved, implying counterintuitive consequences for inference. Nevertheless, Dupont et al. (2022) and Urdangarin et al. (2023) find that RSR performs well on predictive criteria (MSE, WAIC).

RSR removes collinearity and all the spatial variability in the direction of the fixed effect is attributed to the observed covariate. The RSR model has been mostly applied to areal data. In this context, the intrinsic conditional autoregressive (ICAR) model is commonly used to model the spatial random effect $U(s)$ (Hanks et al., 2015). Considering an ICAR prior, $U(s)$ is supported on the nodes $V = 1, 2, ..., n$ of a graph $G = (V, A)$ with A the adjacency matrix where $A_{ij} = 1$ indicates an edge between the nodes $i$ and $j$, while $A_{ij} = 0$ indicates no edge. In this case, the spatial precision matrix (the inverse of the covariance matrix) of $U(s)$ is given by $Q = \Sigma^{-1}$ where $Q = \frac{1}{\sigma^2}(\mathrm{diag}(B) - A)$ where $B$ is a matrix with the number of neighbours of each area and $\sigma$ is the standard deviation. Also, RSR has recently been applied to geostatistical data (Hanks et al., 2015). The spatial random effect, $U(s)$, is modelled using a GRF such as the covariance matrix $\Sigma = \sigma^2\rho(s, s')$, where $\Sigma$ is the Matèrn covariance function, $\mathrm{cov}(U(s), U(s')) = \sigma^2\rho(s, s')$, with:

$$ \sigma^2 \rho(s, s') = \frac{\sigma^2}{\Gamma(\nu)2^{\nu-1}} \left( \frac{\sqrt{2\nu d(s, s')}}{\rho} \right)^{\nu} K_{\nu} \left( \frac{\sqrt{2\nu d(s, s')}}{\rho} \right), \quad (4) $$

where $\sigma^2$ is the marginal variance, $d(s, s')$ is a distance between the locations $s$ and $s'$, $\nu$ is the smoothness parameter, $\rho$ is the range parameter, $\Gamma(\nu)$ is the Gamma function evaluated at $\nu$, $K_{\nu}$ is the modified Bessel function of the second kind (Banerjee et al., 2003). Hanks et al. (2015) constructed an MCMC algorithm to obtain samples from the posterior distributions of the RSR model parameters. This approach considers that the constraint $X^tU(s) = 0$ is first sampled from the unconstrained distribution of $U(s)$. As far as we know, this has not been applied to the R-INLA methodology.

## 2.1.2 | Model extension to geostatistical data using inlabru

The R-inlabru methodology (Lindgren et al., 2011) focuses on models that can be expressed as latent GMRF. A spatial regression model can be inferred in inlabru package, the recent extension of R-inlabru package. This spatial regression model is specified as a Bayesian hierarchical model with latent Gaussian variables where the spatial process $U(s)$ is modelled as a GRF with a Matèrn covariance function, as defined in (4). A GRF with a Matèrn covariance can be expressed as the solution to Stochastic Partial Differential Equation (SPDE) and an approximate solution to the SPDE is obtained via the finite element method that discretises the spatial domain $\mathcal{R}^d$ into non-intersecting triangles, forming a triangulated mesh. The continuously indexed GRF is then approximated by a discretely indexed GMRF, using finite basis functions defined over this mesh (Lindgren et al., 2011). To implement RSR in the context, we must find a $X(s)$ projection onto $U(s)$ to mitigate collinearity between these variables. The idea is to ensure that $(X(s), U(s)) = 0$, where $(\cdot, \cdot)$ represents the scalar product, enforcing orthogonality and eliminating spatial redundancy. The projection uses $P_X$ in (1). Implementing this procedure requires specific adjustments to the inference process to avoid concurrently estimating the spatial random and $X$-related effect within the model. A first approach to implementing RSR in R-INLA is to use the argument `extraconstr=list(A=A, e=zero)` in the spatial function, where A is a matrix of centered covariates weighted by the mesh areas in the SPDE framework. This method ensures the computation of the scalar product, facilitating the removal of orthogonality constraints (Rue & Held, 2005). In what follows, we refer to this implementation as 'RSR extra'. Another approach is to construct the model using a formula-like syntax, explicitly defining the latent effect (Bachl et al., 2019). Unlike the syntax typically used in R-INLA, R-inlabru components provide enhanced functionality, allowing more flexible and complex model specifications, including non-linear predictors. For example, when supplying a formula to R-INLA, the resulting object may refer to the fixed effect $\beta_X$, denoted as $X$. This can lead to ambiguity as it becomes unclear whether $X$ refers to the covariate or its associated effect. To address this, R-inlabru provides a shorthand for constructing basic additive models with fixed and random effects. Here it simplified the implementation of $(I - P_X)U(s)$ by
$U(s) - P_X U(s)$. However, this approach is computationally expensive. To reduce the computational cost, we use the numerical formula:

$$ (I - P_X)U = U - \frac{\sum_{s=s_1}^{s_n} X(s)U(s)}{\sum_{s=s_1}^{s_n} X^2(s)}. $$

This numerical form avoids matrix operations, significantly improving efficiency while maintaining the projection's key properties. This approach balances computational speed and accuracy, which is crucial for large spatial datasets (in what follows, we refer to this implementation as 'RSR formula').

Additionally, as inference in R-INLA is inherently Bayesian, it is important to consider the choice of priors for all the model parameters. This requires specifying prior distributions for the parameters that govern the regression coefficient ($\beta_{RSR}$) and the residual standard deviation parameter ($\sigma_e$), as well as the Matérn $U(s)$'s hyperparameters $\rho$ and $\sigma$ (see (4); $\nu$ is fixed to 1, default value in R-inlabru). For instance, a $N(0, 10^3)$ can be used as a vague prior for the $\beta_{RSR}$ parameter and a Gamma (0.01, 0.01) as a quite non-informative and vague prior for the precision, $\frac{1}{\sigma_e^2}$ (default priors in R-inlabru). For the hyperparameters $\rho$ and $\sigma$, two different priors' specifications were employed for the two RSR implementations.

* **RSR extra:** Following Fuglstad et al. (2019) and the use of penalised complexity (PC) priors in R-inlabru, a prior low probability for the range parameter can be set, for example $P(\rho < \rho_0) = 0.05$ and a prior upper probability for the marginal variance parameter can be set, for example, $P(\sigma > \sigma_0) = 0.05$. Recommended choices for $\sigma_0$ range from 2.5 to 40 times the presumed marginal standard deviation, while $\rho_0$ can be chosen between $\frac{1}{10}$ and $\frac{2}{5}$ of a plausible range. These prior choices are supported by our simulation study and case study (see Sections 3 and 4 for more details).

* **RSR formula:** Initially, we have opted for the default priors in R-inlabru as in 'RSR extra'. However, during the simulation study, we observed that the priors for $U(s)$ needed to be constrained in order to be in accordance with the propriety (2) of the RSR model, if not the estimates of the $\beta_{RSR}$ parameter were too unstable. We then impose the constraints for the PC priors: $P(\rho < \rho_0)$ very close to one and $P(\sigma > \sigma_0)$ very close to zero for $U(s)$ to shrink the spatial field towards a base model with $\rho = 0$ and $\sigma = 0$ rather $\rho = \infty$ and $\sigma = 0$ as proposed by default in R-inlabru.

## 2.2 | gSEM

### 2.2.1 | Model description

The geoadditive structural equation model (gSEM) was developed by Thaden and Kneib (2018) as a method to address spatial confounding by removing spatial dependence from both the response and the covariate before estimating their relationship. Spatial dependence is regressed away from both the response and the covariates, and a

regression involving the residuals only is used to identify the original covariate effect.

The gSEM methodology proceeds in three fitting steps. First, spatial models are fitted separately to both the covariate and the response:

$$ X(s) = f_X(s) + \varepsilon_X(s), \quad \varepsilon_X(s) \stackrel{\text{iid}}{\sim} N(0, \sigma_X^2 I) \quad (5) $$

$$ Y(s) = f_Y(s) + \varepsilon_Y(s), \quad \varepsilon_Y(s) \stackrel{\text{iid}}{\sim} N(0, \sigma_Y^2 I) \quad (6) $$

where $f_X(s)$ and $f_Y(s)$ represent smooth spatial functions estimated using thin plate splines, and $\varepsilon_X(s)$ and $\varepsilon_Y(s)$ are independent Gaussian noise terms with residual standard deviation $\sigma_X$ and $\sigma_Y$, respectively. This approach has been shown to broadly remove bias in effect estimates that arises in traditional spatial models (Thaden & Kneib, 2018).

Then, residuals are extracted from both models to remove spatial structure:

$$ r^X(s) = X(s) - \hat{f}_X(s), \quad r^Y(s) = Y(s) - \hat{f}_Y(s). \quad (7) $$

In a third fitting step, the effect of the covariate on the response is estimated using a linear regression of the residuals:

$$ r^Y(s) = \beta_{gSEM} r^X(s) + \varepsilon(s), \quad \varepsilon(s) \stackrel{\text{iid}}{\sim} N(0, \sigma_e^2 I). \quad (8) $$

where $\beta_{gSEM}$ is the unknown effect of the resulting residual $r^X(s)$ and $\varepsilon(s)$ is an iid Gaussian errors' $\sigma_e$.

In simulation studies, the three regression steps of the gSEM model, unified as a structural equation model, have effectively reduced bias, even under strong spatial confounding in the case of areal data, where standard geostatistical models struggle (Thaden & Kneib, 2018). The framework has been primarily applied to areal data, leveraging spatial indicators and Markov random fields to model regional effects. Its extension to geostatistical data, which captures continuous spatial variation, has been implemented using splines (Dupont et al., 2022). To our knowledge, this has not been extended to R-INLA methodology.

### 2.2.2 | Model extension to geostatistical data using R-inlabru

The two first steps of the gSEM methodology require estimating smooth spatial functions for both the covariate and the response. Traditionally, these functions, $f_X(s)$ and $f_Y(s)$, are estimated using thin plate splines. In a geostatistical framework, $f_X(s)$ and $f_Y(s)$ can be modelled as spatial random fields using R-inlabru's implementation of SPDE, which provides a principled way to capture spatial dependence while maintaining computational efficiency. This avoids the need to explicitly choose smoothing spline penalties. Then extracting the residuals by subtracting the estimated spatial effects from both the covariate and response, remains unchanged. Finally, regressing the residuals $r^X(s)$ on $r^Y(s)$ is conducted in the same manner
as in a traditional linear regression model, but within a Bayesian in-ference framework.

When applying R-inlabru, prior distributions must be specified. Non-informative priors can be assigned to βgSEM and residual standard deviations σχ and σγ (as suggested in the RSR approach), unless domain knowledge suggests otherwise. For the hyperparameters of the Matérn GMRFs used to model fχ(s) and fγ(s), penalised complexity (PC) priors should be employed in R-inlabru. Similarly to the approach outlined for the ‘RSR extra’, the same process is applied for specifying prior distributions of hyperparameters in the gSEM method.

## 2.3 | Spatial+

### 2.3.1 | Model description

The spatial+ method is a two-step regression approach developed by Dupont et al. (2022), which, unlike the gSEM approach, does not remove spatial dependency from the response to retain this useful part of the data. The first step involves fitting a spatial model to the covariate to eliminate its spatial dependence:

$$X(s) = f_X(s) + \varepsilon_X(s), \quad \varepsilon_X(s) \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma_X^2 I), \quad (9)$$

where fχ(s) is an unknown bounded function defined on an open bounded domain ℝd which includes the data locations s = (s1, ..., sₙ) and εχ(s) is an iid errors' noise with residual standard deviation σχ. Note that fχ is a smooth function estimated using thin plate splines.

The resulting residual rX(s), which is the difference between the response covariate X(s) and the estimated spatial residual ŷχ(s), is introduced in a second regression replacing the covariate in the spatial model:

$$\begin{align}
Y(s) &= \beta_{\text{Spatial+}} r^X(s) + f(s) + \varepsilon(s), \nonumber \\
\varepsilon(s) &\stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma_e^2 I), \quad r^X(s) = X(s) - \hat{f}_X(s), \tag{10}
\end{align}$$

where Y(s) is the outcome variable, βSpatial+ is the unknown effect of the resulting residual rX(s), f is a smooth function to be estimated using thin plate splines and ε(s) is an iid errors' noise with residual standard deviation σε. Therefore, the method replaces the covariate with its residuals after accounting for spatial dependence, which has been eliminated by regressing X. Spatial+ preserves the spatial random effect f(s) while addressing collinearity by adjusting the covariate X(s) instead of altering the spatial component of the model. In this model rX(s) has a reduced or completely absent spatial structure and still contains the covariate information.

Dupont et al. (2022) show the good proprieties of their approach for alleviating spatial confounding when traditional spatial models can lead to biased effect estimate due to excessive smoothing of spatial terms. Spatial+ model addresses this issue by separating the spatial dependence from the covariate in the second regression, avoiding the bias caused by smoothing. Notably, in their simulation

study (Dupont et al., 2022), ŷβSpatial+ is, on average, closer to the true β value than the estimate ŷβSpatial of the spatial model, including the covariate and the spatial effects.

The methodology of spatial+ has been developed for geostatistical data as it is formulated through the thin plate spline context. However, Dupont et al. (2022) mention that the spatial+ methodology can be applied to other spatial modelling, including the ICAR model for areal data and GMRF models, but they do not implement them.

### 2.3.2 | Model extension to geostatistical data using R-inlabru

The approach spatial+ requires two regression steps. The first step could be done without R-inlabru using fχ(s) = fχ(s¹, s²), with (s¹, s²) as the coordinates (longitude, latitude) of the data localisation s with fχ(s¹, s²) = (fχ(s₁, s₂), ..., fχ(sₙ, s₂))t as a smooth function to be estimated using thin plate splines. A more complex generalised additive model (GAM) could also be performed for this first regression. Then, the regression of interest (10) would be the unique one fitted using R-inlabru. Then the use of R-inlabru requiring Bayesian inference necessitates just one delicate step to fix the priors of all parameters. After it stays to infer a regression model from geostatistical data for which the R-inlabru package is adapted. The R-INLA methodology that captures the continuous space by SPDE and finite element method avoids the choice of the smoothing spline penalties. Then the R-inlabru could be a suitable option preventing the difficult smoothness selection using splines (Dupont et al., 2022).

Concerning the priors, as for gSEM method, default priors in R-inlabru as in ‘RSR extra’ can be chosen for the covariate parameter, here βSpatial+, and the two residual standard deviations σχ and σε of the two regression models, unless outside information is known as expert knowledge. For the hyperparameters of the Matérn GMRFs of the two regression models, we adopt again the approach detailed in the ‘RSR extra’ Section 2.1. As for the gSEM approach, our findings support these priors in the simulation and case studies (see Sections 3 and 4).

## 2.4 | Spatial+ 2.0

### 2.4.1 | Model description

Recently, Urdangarin et al. (2023, 2024) proposed a simplified version of the spatial+ approach, which avoids the two regression steps, that is, one does not fit a spatial regression model to the covariate of interest. This alternative to removing the spatial dependence in the covariate is achieved by decomposing the covariate X(s) into a linear combination of the eigenvectors of the spatial random effect precision matrix Q. The decomposition allows a unique regression that retains short-scale dependence between
<figure><img src="image_4.png" /></figure>

the covariate and the outcome. Here we define $Q = S\Delta S'$ considering the spectral decomposition of Q, where S is the orthogonal matrix of Q's eigenvectors and Δ is a $n \times n$ diagonal matrix with the eigenvalues of Q in the main diagonal. We then arrive at

$$X(s) = a_1 S_1 + \dots + a_s S_s. \quad (11)$$

Assuming that the collinearity between the fixed and random effects primarily arises from the eigenvectors associated with the lowest non-null eigenvalues, the covariate is split into two parts:

$$X(s) = Z(s) + Z^*(s), \quad (12)$$

where $Z^*(s) = a_{s-k}S_{s-k} + \dots + a_sS_s$ comprises large-scale eigenvectors associated with the lowest eigenvalues and $Z(s) = a_1S_1 + \dots + a_{s-(k+1)}S_{s-(k+1)}$ represents the remaining eigenvectors. Then, a spatial regression model is fitted to the outcome variable $Y(s)$:

$$Y(s) = \beta_{\text{Spatial}+2.0} Z(s) + U(s) + \varepsilon(s), \quad (13)$$

where $Z(s)$ is the spatially decorrelated part of $X$, $U(s) \sim N(0, \Sigma)$ is the spatial random effect with mean zero and covariance matrix $\Sigma$ and $\varepsilon(s) \sim N(0, \sigma^2I)$ is the iid errors' noise with $\sigma_\varepsilon$ as standard deviation.

Urdangarin et al. (2024) show in a simulation study that the spatial+2.0 model is effective in accurately recovering the true value of the fixed effect regarding the spatial and RSR models. However, the choice of the number of eigenvectors to keep (or remove) is delicate to have accurate $\beta_{\text{Spatial}+2.0}$'s estimate proprieties. From the simulation and case study, they recommend removing 7% and 20% of the eigenvectors if the spatial covariate pattern is smooth and between 35% and 43% if the covariate pattern is less smooth from their simulation and case studies.

In Urdangarin et al. (2024), spatial+2.0 is applied to areal data, where model fitting and inference are carried out from a full Bayesian perspective using two main techniques: MCMC and INLA. The spatial effect $U(s)$ is modelled using a multivariate intrinsic conditional autoregressive (ICAR) prior with a fixed correlation parameter, or proper conditional autoregressive (PCAR) (Jin et al., 2007) or BYM2 prior (Riebler et al., 2016). As in the spatial+ framework, $U(s)$ could be modelled using spline as described in Section 2.3. However, its extension spatial+2.0 has not yet been applied in a geostatistical data context.

### 2.4.2 | Model extension to geostatistical data using R-inlabru

The spatial+2.0 approach fits a linear model where the covariate $X(s)$ has been previously split into two parts and then uses k-selected eigenvectors of the spatial precision matrix Q as regressors. In the geostatistical data context, as in the areal data case, selecting the optimal value of k is difficult, especially as, in the continuous spatial case, the dimension of X is generally larger (being the number of measurements over the space). As proposed by Urdangarin et al. (2023), it involves balancing the number of included eigenvectors with

minimising the smallest eigenvalue with the effectiveness of the Bayesian WAIC criterion using for model selection.

For the parameter $\beta_{\text{Spatial}+2.0}$ and the residual standard deviations $\sigma_\varepsilon$, default priors in R-inlabru as described in 'RSR extra'. Regarding the hyperparameters of the Matèrn GMRF in the regression model, PC priors are used and we again fix them following the recommendations of Fuglstad et al. (2019). Our findings support these priors in the simulation and case studies (see Sections 3 and 4).

## 3 | SIMULATION

In this section, we aim to assess the behaviours of the four resolution methods extended to geostatistical data using R-inlabru in a scenario where a spatial confounding is simulated. We hope the two extensions of the RSR approach offer inference very similar to the null model as expected and for the three others, to estimate the true value of the covariate parameter with a good empirical coverage and we compare our results to those obtained by these methods implemented with splines if it has been proposed in the literature to see if our extensions do as well. In the simulated scenario, the true relationship between the covariate and the response variable ($\beta$) is known and the random effect $U(s)$ follows a GMRF with a Matèrn covariance matrix as obtained by generating the data by R-inlabru. We compare the models primarily on the basis of the accuracy of $\beta$'s estimates and the Bayesian model selection criteria, DIC and WAIC (Gelman et al., 2014) which inform about the model predictive performance.

### 3.1 | Data

We generate 50 independent datasets, constructing each with an observed covariate $\mathbf{x} = (x(1), \dots, x(n))'$ and the response dataset $\mathbf{y} = (y(1), \dots, y(n))'$ such that:

$$x(s) = 0.35z_x(s) + \varepsilon_x(s),$$

$$y(s) = \beta x(s) + u(s) + \varepsilon_y(s),$$

where $s = 1, \dots, n$ with $n = 500$ randomly selected locations within a continuous spatial field $[0, 10] \times [0, 10]$, $z_x = (z_x(1), \dots, z_x(n))^t$ is the spatial residual of the covariate collected at the $n = 500$ locations with a Matèrn short-scale spatial correlation, $\varepsilon_x(s) \sim N(0, \sigma_x^2I)$ with $\sigma_x^2 = 0.1, \beta = 3$ is the true covariate effect, $u(s) = z_u(s) - z_x(s)$ are the spatial effects where $z_u = (z_u(1), \dots, z_u(n))^t$ is the spatial random effect collected at the $n = 500$ locations with a Matèrn large-scale spatial correlation, $\varepsilon_y(s) \sim N(0, \sigma_y^2I)$ with $\sigma_y^2 = 1$. This approach to creating spatial confounding from geostatistical data is close to Thaden and Kneib (2018), and in turn inspired the method proposed by Dupont et al. (2022). Here, the continuous spatial random fields have a Matèrn spatial covariance matrix constructed using an R-INLA approach. This involves creating for each field a
triangulated mesh which is the support of the finite element method's solution, a GMRF approximative solution, of the SPDE coupled with the parameters of the Matérn covariance function. We employ two meshes of 1283 nodes, the ratio between the number of observations and the mesh's size being 0.39. We see with other mesh sizes that this ratio is a acceptable trade-off between the accuracy of the GMRF representation and the computational cost. We put $\theta_1 = 2$ and $\theta_2 = 0.4$ for $z_x$, and $\theta_1 = 0$ and $\theta_2 = 1$ for $z_u$ in R-INLA ($\theta_1 = \log(\sigma) - \log(\sigma_0)$ and $\theta_2 = \log(\rho) - \log(\rho_0)$) to control our two GRFs. Note that vague PC priors have been fixed creating the datasets: $P(\rho < 0.5) = 0.5$ and $P(\sigma > 0.5) = 0.5$. Alternative simulation setups could be suggested to study the influence of the data generation process on the performance of the different resolution methods as suggested by Marques and Kneib (2022) and Schmidt (2022): the correlation between covariate and spatial random effect, the marginal variances and scales for spatial processes and the spatial model fits if multiple steps (spatial+ and gSEM) could be investigated.

## 3.2 | Models

We fit six regression models to each of the 50 datasets comprised of the response variable Y and the covariate of interest X using R-inlabru. By default, vague priors are given to the parameter $\beta$ and the residual precision $1/\sigma_e$ following a $N(0, 10^3)$ and a Gamma (0.01, 0.01) respectively. The models are as follows:

1. A linear regression model (null model) with no spatial effect given by

$$Y(s) = \beta_{\text{Null}}X(s) + \varepsilon(s), \quad \varepsilon(s) \sim N(0, \sigma_e^2 I), \quad (14)$$

where $\beta_{\text{Null}}$ is the unknown fixed effect to estimate and $\varepsilon(s)$ the iid errors' noise with the standard deviation $\sigma_e$.

2. A spatial model with a spatial random effect given by

$$Y(s) = \beta_{\text{Spatial}}X(s) + U(s) + \varepsilon(s), \\ \varepsilon(s) \sim N(0, \sigma_e^2 I), \quad U(s) \sim N(0, \sigma_u^2 \rho(s, s')), \quad (15)$$

where $\beta_{\text{Spatial}}$ is the covariate fixed effect to estimate, $\varepsilon(s)$ is an iid errors' noise with the standard deviation $\sigma_e$, $U(s)$ is a spatial random effect modelled as a GRF with mean zero and Matérn covariance matrix. PC priors for the Matérn hyperparameters are taken as $P(\rho < 0.05) = 0.05$ and $P(\sigma_u > 3) = 0.05$.

3. An RSR model given in (1) with $\Sigma = \sigma^2\rho(s, s')$ given in (4), where $v$ has been fixed to 1. For the 'RSR extra' approach, we adopted the same vague priors as in the spatial model, following the methodology of Fuglstad et al. (2019). In the 'RSR formula' approach, PC priors for the Matérn hyperparameters are taken

as $P(\rho < 15) = 0.9999$ and $P(\sigma > 1.5) = 0.0001$ as mentioned in Section 1.

4. A gSEM model given in (5), (6) and (8)

$$f_X(s), f_Y(s) \sim N(0, \Sigma), \quad \Sigma \text{ given in (4) with } v = 1$$

PC priors for the Matérn hyperparameters of $f_X(s)$ and $f_Y(s)$, for X and Y spatial regression, are taken as in the spatial model.

5. A spatial+ model given in (9) and (10)

$$f_X(s), f_Y(s) \sim N(0, \Sigma), \quad \Sigma \text{ given in (4) with } v = 1$$

PC priors for the Matérn hyperparameters of $f_X(s)$, for X regression, are taken as $P(\rho^f < 0.01) = 0.01$ and $P(\sigma^f > 0.15) = 0.01$. For the second step of Y regression, the hyperparameters priors associated with $f(s)$ are the same as in the spatial model.

6. A spatial+2.0 model given by (13) with

$$\Sigma = \sigma^2 \rho(s, s'), \quad v = 1$$

PC priors for the Matérn hyperparameters are taken as for the spatial model. We retain 400 of the 500 available eigenvectors based on minimising the WAIC criterion calculated for each of the 500 possible models including from 1 to 500 eigenvectors of the spatial random effect (see Section 2.4). This procedure was replicated 50 times, consistently converging around 400 eigenvectors for each dataset.

For the methods implemented with splines, our data, unlike those of Dupont et al. (2022) which guarantee that the fields lie within the span of the spatial basis vectors used in the models, were not spanned (Marques & Kneib, 2022). The thin plate spline function, referred to as a GAM function hereafter, was then applied with a smoothing penalty to be estimated.

## 3.3 | Results

Figure 1 presents box plots of the posterior means of the estimated $\hat{\beta}$ for the null model, the spatial model and the five approaches fitted to the 50 data replicates within R-inlabru. The results within R-inlabru are shown on the left, while, after the vertical dashed line, the estimated means of $\hat{\beta}$ (frequentist approach) obtained from the three models using a GAM approach (as proposed in the literature). Table 1 summarises, for each approach, the mean of the posterior means or estimated means ($\hat{\beta}$), the empirical standard deviation (ESD) of the posterior means or estimated means, the mean estimated standard errors (SE), the deviance information criterion (DIC) values when applicable and the information criterion (WAIC) values when applicable. Notably, WAIC and AIC are not reported for the gSEM model, as the final regression is performed on the
<figure><img src="image_5.png" /><figcaption>FIGURE 1 Simulation results: Estimated covariate effect $\hat{\beta}$ (posterior mean or estimated mean) for each model and data replicate under different spatial modelling frameworks. The first approaches (on the left of the dashed vertical line) use R-inlabru, while the second (on the right of the dashed vertical line) employ GAM with a nonzero smoothness parameter (estimated at 0.5). The true covariate effect ($\hat{\beta}$ = 3) is represented by the red horizontal dashed line.</figcaption></figure>

<table><caption>TABLE 1 Table 1 summarises the mean of the posterior means or estimated means, $\hat{\beta}$, the empirical standard deviations of the posterior mean, ESD, the mean estimated standard errors, $\overline{SE}$, the mean of the DIC values, the mean of the WAIC values, and the empirical coverage of 95% credibility intervals, CI, over the 50 replicated data, for each model when applicable.</caption><thead><tr><th>Model</th><th>β̂</th><th>ESD</th><th>SĒ</th><th>DIC</th><th>WAIC</th><th>CI(%)</th></tr></thead><tbody><tr><td>Null</td><td>1.031</td><td>0.136</td><td>0.149</td><td>1582.51</td><td>1582.529</td><td>0</td></tr><tr><td>Spatial</td><td>2.649</td><td>0.800</td><td>0.303</td><td>1498.956</td><td>1500.482</td><td>76</td></tr><tr><td>RSR extra</td><td>1.031</td><td>0.136</td><td>0.149</td><td>1580.524</td><td>1580.525</td><td>0</td></tr><tr><td>RSR formula</td><td>1.031</td><td>0.136</td><td>0.149</td><td>1580.524</td><td>1580.525</td><td>0</td></tr><tr><td>gSEM</td><td>4.114</td><td>0.660</td><td>0.474</td><td></td><td></td><td>34</td></tr><tr><td>Spatial+</td><td>2.781</td><td>1.063</td><td>0.303</td><td>1501.736</td><td>1503.327</td><td>74</td></tr><tr><td>Spatial+2</td><td>2.999</td><td>0.375</td><td>0.412</td><td>1503.197</td><td>1507.389</td><td>96</td></tr><tr><td>GAM spatial</td><td>3.335</td><td>0.298</td><td>0.301</td><td></td><td></td><td>80</td></tr><tr><td>GAM gSEM</td><td>4.424</td><td>0.571</td><td>0.480</td><td></td><td></td><td>14</td></tr><tr><td>GAM spatial+</td><td>4.722</td><td>0.495</td><td>0.540</td><td></td><td></td><td>8</td></tr></tbody></table>

residuals of Y rather than on Y itself, making these criteria incomparable those of the other models (Dupont et al., 2022; Thaden & Kneib, 2018). Additionally, Table 1 reports the empirical coverage of 95% credibility intervals (CI) over the 50 replicated datasets for all approaches.

If one first examines the null and spatial models, one observes very different results as expected due to the superposed variability created in each dataset between the covariate and the spatial random effect (spatial confounding). Here, one observes that the spatial

models are closer to the true value with a low coverage at 76% (the estimates are in a large range (ESD = 0.8 and see Figure 1)).

In this context of geostatistical data and R-INLA methodology, the both RSR models (extra and formula) preserve the theoretical properties outlined in (2) and (3) with the PC prior we have taken for the spatial random effect: all the posterior means of $\beta_{RSR}$ and $\beta_{Null}$ are identical (see Figure 1 and Table 1) and the both RSR models have smaller estimated standard errors, $\overline{SE}$, than the spatial models, indicating reduced estimated variance for $\beta_{RSR}$ compared to $\beta_{Spatial}$. The both RSR
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 2 Concentration maps of Cd across France over 445 sites (BRAMM sites) in 2016: Cd concentration in mosses, measurements in µg/g (left); Cd concentrations in air, µg/m<fnref n="3" />, from the physical model EMEP values (right).</figcaption></figure>

remains overly restrictive for addressing spatial confounding, even if its WAIC values are slightly lower than those of the null models. The RSR extra computation takes three times longer than the RSR formula to compute. The gSEM models, which are fitted in R-inlabru through three steps, yield a mean $\hat{\beta}$ estimate of 4.114. This estimate deviates more from the true $\beta$ compared to $\hat{\beta}_{Spatial}$ and exhibits a low empirical coverage of just 34%. The spatial+ models, the two regression steps are fitted here in R-inlabru, produce a mean $\hat{\beta}$ estimate of 2.781 which is on average closer to the true $\beta$ than $\hat{\beta}_{Spatial}$ with a low empirical coverage of 74%. In Figure 1 this estimate range is the largest (ESD=1.063 in Table 1). The spatial+2.0 models, where 400 eigenvectors were retained as explained in Section 3.2, yield a mean estimated $\hat{\beta}$ of 2.999, with an empirical coverage of 96%, and the range of their estimates is narrower than those of the spatial and spatial+ models (ESD=0.375 in Table 1 and Figure 1). Comparing all the DIC and WAIC values, the spatial, spatial+ and spatial+2.0 models give almost similar results.

In this simulation study, the R-inlabru approach outperforms the GAM approach for gSEM and spatial+ models. The GAM approach is not favoured because the data are not spanned, unlike in Dupont et al. (2022).

# 4 | CASE STUDY: APPLICATION ON TERRESTRIAL MOSSES

## 4.1 | Data

We apply the proposed extensions in a spatially confusing case we encountered studying the link between Cd concentrations measured in terrestrial mosses throughout France and Cd concentrations predicted in air in France by an EMEP<fnref n="3" /> model. This is an atmospheric transport model used to predict metal concentrations across Europe.

We analyse a dataset of mosses collected in 2016 as part of the BRAMM<fnref n="4" /> scheme, which covered 445 sites across France (Leblond

et al., 2011). At each of the 445 sites, s, a sample mixture of mosses of a unique species over s (4km²) is collected and Y(s) its Cd concentration, is measured in the laboratory. Using a linear regression model, we relate this concentration to an air Cd concentration value X(s) obtained at location s, from the values predicted by the EMEP model for 2016 on a 0.1°×0.1° longitude-latitude grid. To obtain the value at the mosses' location, X(s), we used the nearest EMEP value from the site s. Figure 2 shows the values Y(s) and X(s) at each location s.

We initially fitted a null model, and the parameter $\beta_{Null}$ was highly positively significant ($p < 5.07e^{-7}$), but a residual spatial autocorrelation of the response values was still there (indicating from Moran's index z-scores, $p < 2.2e^{-16}$). We then fitted a spatial model to the dataset, and the parameter $\beta_{Spatial}$ was no longer significant. These findings were unexpected, so then we turned to models that could resolve potential spatial confounding. In the next section, we provide details of the fitted models for mitigating spatial confounding.

## 4.2 | Models

As for the simulation study, we used the default R-INLA priors for the parameters $\beta_0$ and $\beta$, and the residual standard deviations $\sigma_e$.

1. A null model as defined in (14) with a fixed intercept $\beta_{0Null}$ where Y(s), s ∈ {1, ..., 445}, are a logarithmic transformation of Cd concentrations of mosses to reduce the skewness of their distribution, X(s), are the Cd concentrations in air obtained from the EMEP's values as described above.

2. A spatial model as defined in (15) with a fixed intercept $\beta_{0Spatial}$ with Y(s), X(s)) defined as in the null model above. PC priors for the Matérn hyperparameters of the spatial random effect are taken as: $P(\rho < 6e^4) = 0.05$ and $P(\sigma > 5) = 0.05$.

3. A RSR model given by (1) with a fixed intercept $\beta_{0RSR}$, $\Sigma = \sigma^2\rho(s, s')$ given in (4), where $\nu$ has been fixed to 1, with Y(s), X(s)) defined as in the null model above. For the RSR extra approach, the same PC prior as in the spatial model is used for the Matérn hyperparameters. For the RSR formula approaches, PC priors for the

<footnote n="3">³https://www.emep.int/</footnote>

<footnote n="4">⁴https://bramm.mnhn.fr/</footnote>
<figure><img src="image_8.png" /></figure>

<table><caption>TABLE 2 The posterior mean of the covariate effect $\beta$, quantiles at 2.5% and 97.5% of $\beta$, DIC value, WAIC value and estimated standard error of the posterior mean (SE) for each model fitted in R-inlabru for the case study.</caption><thead><tr><th>Model</th><th>β̂</th><th>0.025quant</th><th>0.975quant</th><th>DIC</th><th>WAIC</th><th>SE</th></tr></thead><tbody><tr><td>Null</td><td>2.306</td><td>1.422</td><td>3.191</td><td>769.07</td><td>775.49</td><td>0.451</td></tr><tr><td>Spatial</td><td>0.230</td><td>-0.820</td><td>1.280</td><td>628.45</td><td>629.28</td><td>0.536</td></tr><tr><td>RSR extra</td><td>2.306</td><td>1.422</td><td>3.191</td><td>769.07</td><td>775.49</td><td>0.451</td></tr><tr><td>RSR formula</td><td>2.302</td><td>1.415</td><td>3.189</td><td>773.21</td><td>778.21</td><td>0.451</td></tr><tr><td>gSEM</td><td>0.118</td><td>-0.569</td><td>0.806</td><td>—</td><td>—</td><td>0.351</td></tr><tr><td>Spatial+</td><td>0.213</td><td>-0.837</td><td>1.258</td><td>628.84</td><td>629.57</td><td>0.535</td></tr><tr><td>Spatial+2.0</td><td>0.748</td><td>-0.845</td><td>2.341</td><td>628.41</td><td>628.25</td><td>0.813</td></tr><tr><td>GAM spatial</td><td>0.166</td><td>-0.957</td><td>1.288</td><td></td><td></td><td>0.573</td></tr><tr><td>GAM gSEM</td><td>1.273</td><td>-3.889</td><td>6.435</td><td></td><td></td><td>2.634</td></tr><tr><td>GAM spatial+</td><td>1.269</td><td>-3.615</td><td>6.152</td><td></td><td></td><td>2.492</td></tr></tbody></table>

Matérn hyperparameters are set as $P(\rho < 18e^6) = 0.9999$ and $P(\sigma > 2.5) = 0.0001$. Note that in the projection $P_X$, X includes the ones of the intercept.

4. A gSEM model given by (5), (6) and (8), with a fixed intercept $\beta_{0gSEM}$ in the third regression model, $f_X(s)$ and $f_Y(s) \sim N^*(0, \Sigma)$, with $\Sigma$ given in (4) with $\nu = 1$, and $(Y(s), X(s))$ defined as in the null model above. PC priors for the Matérn hyperparameters of $f_X(s)$ and $f_Y(s)$, for X and Y spatial regressions, are taken as the spatial model.

5. A spatial+ model given by (9) and (10), with a fixed intercept $\beta_{0Spatial+}$ in the second regression model, $f_X(s)$ and $f_Y(s) \sim N^*(0, \Sigma)$, with $\Sigma$ given in (4) with $\nu = 1$ and $(Y(s), X(s))$ defined as in the null model above. PC priors for the Matérn hyperparameters of $f_X(s)$, for the first regression, are taken as $P(\rho^{fX} < 13e^6) = 0.05$ and $P(\sigma^{fX} > 13) = 0.05$. For the second regression, the same PC prior as in the spatial model are used for the Matérn hyperparameters. 6. A spatial+2.0 model given by (13) with a fixed intercept $\beta_{0Spatial+2.0}$, $\Sigma = \sigma^2\rho(s, s')$ given in (4), where $\nu$ has been fixed to 1, with $(Y(s), X(s))$ defined as in the null model above. PC priors for the Matérn hyperparameters are taken as in the spatial model. We retain 310 of the 445 available eigenvectors (35% removed) based on minimising the WAIC criterion.

As with the simulation part, some methods in this section are implemented using splines. Unlike the data in Dupont et al. (2022), which lie within the span of the spatial basis vectors, our data were not spanned (Marques & Kneib, 2022). We then applied the thin plate spline function, referred to as a GAM function, with an estimated smoothing penalty.

## 4.3 | Results

Table 2 gives the posterior mean of the parameter $\beta$, the 95% CI, DIC and WAIC for each model fitted using R-inlabru. As stated in the simulation results Section 3.3, for the gSEM model we do not

display the results concerning the DIC and WAIC. As mentioned, the covariate effect is significant with the null model and no longer with the spatial model. The RSR model has results very close to those of the null model. Spatial+ and gSEM results are close to the spatial model. The spatial+2.0 model has a higher estimate value for the covariate effect but is still non-significant. DIC and WAIC criteria give similar results, preferring spatial, spatial+ and spatial+2.0 models, with the spatial+2.0 model having the smallest value.

For the spatial+2.0 model, we selected a subset of eigenvectors based on the WAIC values of the 445 different models, splitting $X(s)$ into two parts. Figure 3 gives these different values for all models. The model with the smallest WAIC is chosen, that is, 310 eigenvectors retained (the smallest WAIC values with low eigenvectors [0, 20] retained in the models have been excluded as suggested by Urdangarin et al. (2023) because these models could be too conservative inflating parameters standard errors; Figure 4).

Figure 5 displays the posterior medians prediction (ex- $p(0.50quant(\hat{Y}))$) maps for Cd concentrations in mosses, obtained from the null, spatial, RSR and spatial+ models. A grid consisting of 139,673 pixels of $4 \text{ km}^2$ was constructed to cover the map of France. In the gSEM model, the resulting residual-based analysis does not retain the spatial structure of the response variable Y, making interpolation mapping unfeasible. Also, the spatial+2.0 method necessitates the spectral decomposition of the spatial precision matrix $Q$, a $139,673 \times 139,673$ matrix for which this decomposition cannot be obtained numerically. Given that we cannot produce the prediction map with the method at such a resolution, we attempted a grid of 30,000 pixels, but it was still too large. Decreasing the map resolution again would render the maps too approximate. As expected from the model estimates, the null and RSR models produce identical prediction maps, while the spatial and spatial+ maps are very similar. The latter two maps more accurately highlight a northeastern zone with relatively high values, consistent with observed data (see Figure 2 (left)). In contrast, the null and RSR maps predict elevated values around the city of Marseille in the southeast of France, which are found in the $X(s)$ values (see Figure 2 (right)). By incorporating a
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 3 WAIC value of the spatial+2.0 models versus the number *k* of eigenvectors retained in the model.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 4 Beta estimated value box plots for the spatial+2.0 models versus the number *k* of eigenvectors retained in the model.</figcaption></figure>

larger spatial structure, the spatial and spatial+ models capture more
spatial heterogeneity (Figure 5 (right)).

Figure 6 presents the posterior 95% uncertainty width (exp(0.975quant(Y)) - exp.(0.025quant(Y))) maps of Cd predicted concentrations. These maps use the same grid as those in Figure 5 with the same constraints to generate the spatial+2.0 map. The 95% uncertainty width maps produced by the null and RSR models are very close, and the spatial and spatial+ maps display moderately similar patterns. Indeed the spatial+ map reveals more areas with high uncertainty widths.

5 | DISCUSSION

This paper explores how recent resolution methods for spatial con-
founding (RSR, gSEM, spatial+, spatial+2.0) can be extended to the

context of geostatistical data modelled by latent GMRF in an R-INLA
methodology. We implement them using R-inlabru, the recent R-
package which facilitates the use of R-INLA for modelling ecologi-
cal data. We have generated simulated datasets in a geostatistical
data context of spatial confounding to compare the results of the
extended resolution methods. Our extensions find the expected
properties of the methods. The both RSR (extra and formula) model
gives results close to the null model, as it assigns spatial variability
to the covariate. To obtain this result, the PC prior to the range pa-
rameter for the spatial random effect of the RSR formula method has
been set to shrink the model to a baseline model with a range of zero,
otherwise the estimation of the effect of the covariate was highly
unstable (results not shown). In our results, the Bayesian criteria
(WAIC, DIC) of the RSR model are similar to those of the null model,
not better contrary to what was found by Urdangarin et al. (2023),
perhaps because the RSR model does not fit well the residual spatial
Null

<figure><img src="image_12.png" /></figure>

Spatial

<figure><img src="image_13.png" /></figure>

RSR

<figure><img src="image_14.png" /></figure>

Spatial +

<figure><img src="image_15.png" /></figure>

FIGURE 5 Cd concentration mosses' prediction maps with a common scale between 0 and 1.1 µg/g; null model (range of the Cd prediction values (medians) is [0.12; 0.98] (µg/g)); spatial model (range of the Cd prediction values (medians) is [0.04; 0.41] (µg/g)); RSR model (range of the Cd prediction values (medians) is [0.12; 1.1] (µg/g)); spatial+ model (range of the Cd prediction values (medians) is [0.04; 0.39] (µg/g)).

dependence of Y, an important point of attention to have in the context of RSR models as mentioned by Khan and Calder (2022). The gSEM model tends to overestimate the beta coefficient. This poor behaviour may be due to the fact that the non-unified version of gSEM is implemented here, not gSEM as a structural equation model. More generally regarding the simulation study, the results were fairly dependent on the PC priors on the spatial random effects, showing how difficult it is to separate spatial patterns. Finally the spatial+2.0 has the best 95% empirical coverage (96%, a little bit narrow) and the nearest estimate of the true fixed effect $\beta$. Our simulation study confirms that the gSEM and spatial+ model provides estimates closer to the true beta for the R-inlabru than GAM.

We were interested in a real case study of terrestrial mosses to explore the relationship between Cd concentration in mosses and air Cd concentration given by the EMEP physical model. This subject

can be an important issue to have reliable proxy sources of air pollution for population-based epidemiological studies investigating the health effects of individual metals or groups of metals (Lequy et al., 2023). In this context, we have explored the four methods cited above, but were unable to find a significant covariate effect; however, we observed that spatial+2.0 increases the effect of the covariate (but it does not become significant). Maybe there would be the possibility of a better eigenvector separation criterion than the WAIC criterion, based on expert opinion on the spatial scale of the covariate effect. From the bryologist's perspective, such a scale could be difficult to define a priori, especially because of localised effects that could change the effect of the covariate on the territory. It would, therefore, be interesting to consider also the inclusion of other covariates in the model, such as the wind speed and/or the tree canopy on the sampling site.
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 6 Cd concentration mosses' posterior range at 95% maps with a common scale between 0 and 1.5 µg/g: null model (range of the Cd uncertainty widths is [0.012; 1.461] (µg/g)); spatial model (range of the Cd uncertainty widths is [0.029; 0.420] (µg/g)); RSR model (range of the Cd uncertainty widths is [0.013; 1.144] (µg/g)); spatial+ model (range of the Cd uncertainty widths is [0.041; 0.372] (µg/g)).</figcaption></figure>

Our data consist of geostatistical data with Gaussian responses, but one might also want to consider count data over a continuous domain (or other type of data considered in the GLM framework). It could be handled using the INLA framework because the R-inlabru package easily handles generalised linear models (GLM). However, gSEM and spatial+ implicitly assume that residuals are readily available, typically defined as the difference between observed and predicted values. While this assumption holds naturally for continuous data with additive noise, it becomes less straightforward for non-Gaussian responses. In such cases, residuals must be defined differently. As discussed by McCullagh (1989), various types of residuals, such as Pearson or deviance residuals, can be used to extend the applicability within a GLM framework. In addition, by multiplying successive regression steps, the gSEM and spatial+ approaches omit consideration of the additional uncertainty arising from the first-step estimations determining the residuals. This can be dealt with by joint modelling, but Schmidt (2022) has raised the issue of possible feedback from Y to X when using a single inferential procedure. Marques

and Wiemann (2023) attempt to overcome the problem in a Bayesian version of spatial+, but the new modified version of spatial+, referred to here as spatial+2.0, proposed by Urdangarin et al. (2024) with a single regression, means that the problem no longer seems as significant.

Other spatial confounding mitigation methods, such as TGMRF (Prates et al., 2015) or other spatial spline modellings than GAM as done here, are proposed in literature (Zaccardi et al., 2025). The implementation of the TGMRF method within the R-INLA framework remains challenging due to its departure from the latent Gaussian model structure and the need to treat model parameters as hyperparameters, requiring significant methodological developments beyond the scope of this study. Also, more complex spatial modelling could be an issue to deal with spatial confounding, such as random-effects GAMs or SVC (spatial varying coefficient) models. These models offer greater flexibility in the relationship between the covariate and the response in the regression framework. As highlighted by Sørbye et al. (2019), mitigating spatial confounding requires explicit assumptions about spatial variation and variance
<figure><img src="image_18.png" /></figure>

partitioning, as modelling choices fundamentally influence inference and the interpretation of spatial relationships.

In conclusion, mitigating spatial confounding is a critical challenge for spatial statistics, because when covariates and response have close spatial patterns, it is difficult to distinguish what belongs to the covariates from what does not. Unlike longitudinal or cluster data, there are no replications (intra-individual or intra-cluster), in general (spatio-temporal data offer some possibility in this respect) to distinguish spatial patterns. Here, everything is determined by a single data point at location s and the dependence between ($s_1, ..., s_n$). The answer lies partly with the expert who knows his/her field and can best specify models (spatial scale of the covariates and parameter priors). Nevertheless the proposed methods offer a contribution to how spatial statistics can approach the challenge of spatial confounding in a geostatistical data context in order to reach more valid conclusions.

## AUTHOR CONTRIBUTIONS

Jérémy Lamouroux, Alizée Geffroy and Isabelle Albert contributed conceptually and formulated the models. Jérémy Lamouroux and Alizée Geffroy provided the practical implementation of the model. All the authors assisted in writing and editing the manuscript.

## ACKNOWLEDGEMENTS

The authors wish to acknowledge all the collectors who contributed to the BRAMM survey and those who took care of the samples in the laboratory. The authors would like to thank the staff of the USRAVE laboratory (INRAE-Centre de Bordeaux), accredited by the French Committee for Accreditation, for their analyses of Cd accumulated by mosses. The authors are also thankful to Ilia Ilyin of MSC-East of EMEP for providing the modelled atmospheric deposition and concentration of Cd.

## FUNDING INFORMATION

Open access funding provided by Université Paris-Saclay. The French Agency for Ecological Transition (conventions ADEME-MNHN N°1562c0027 and 1666C0022) supported Mosses collection and metal concentration analysis. Jérémy Lamouroux received a PhD grant from the National Research Institute for Agriculture, Food, and the Environment (INRAE) in France and The French National Museum of Natural History (MNHN).

## CONFLICT OF INTEREST STATEMENT

All authors have no conflict of interest to declare.

## PEER REVIEW

The peer review history for this article is available at https://www.webofscience.com/api/gateway/wos/peer-review/10.1111/2041-210X.70106.

## DATA AVAILABILITY STATEMENT

The localisation data are available via GBIF: https://doi.org/10.15468/c2bqic (Inventaire National du Patrimoine Naturel, 2017).

Data on Cd concentrations in mosses and in the air, as well as the source code, are available via Zenodo: https://doi.org/10.5281/zenodo.15748177 (Lamouroux et al., 2025).

## ORCID

Jérémy Lamouroux https://orcid.org/0009-0001-4410-4771

## REFERENCES

Bachl, F. E., Lindgren, F., Borchers, D. L., & Illian, J. B. (2019). inlabru: An R package for Bayesian spatial modelling from ecological survey data. *Methods in Ecology and Evolution, 10*(6), 760–766. https://doi.org/10.1111/2041-210X.13168

Bakka, H., Rue, H., Fuglstad, G.-A., Riebler, A., Bolin, D., Illian, J., Krainski, E., Simpson, D., & Lindgren, F. (2018). Spatial modeling with r-inla: A review. *Wiley Interdisciplinary Reviews: Computational Statistics, 10*(6), e1443.

Banerjee, S., Carlin, B. P., & Gelfand, A. E. (2003). *Hierarchical modeling and analysis for spatial data*. Chapman and Hall/CRC.

Bates, J. W. (1992). Mineral nutrient acquisition and retention by bryophytes. *Journal of Bryology, 17*(2), 223–240.

Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems. *Journal of the Royal Statistical Society, Series B: Statistical Methodology, 36*(2), 192–225.

Besag, J., York, J., & Mollié, A. (1991). Bayesian image restoration, with two applications in spatial statistics. *Annals of the Institute of Statistical Mathematics, 43*, 1–20.

Blangiardo, M., & Cameletti, M. (2015). *Spatial and Spatio-temporal Bayesian Models with R-INLA*. John Wiley & Sons, Ltd.

Clayton, D. G., Bernardinelli, L., & Montomoli, C. (1993). Spatial correlation in ecological analysis. *International Journal of Epidemiology, 22*(6), 1193–1202.

Dupont, E., Wood, S. N., & Augustin, N. H. (2022). Spatial+: A novel approach to spatial confounding. *Biometrics, 78*(4), 1279–1290.

Fuglstad, G.-A., Simpson, D., Lindgren, F., & Rue, H. (2019). Constructing priors that penalize the complexity of gaussian random fields. *Journal of the American Statistical Association, 114*(525), 445–452.

Gelman, A., Hwang, J., & Vehtari, A. (2014). Understanding predictive information criteria for bayesian models. *Statistics and Computing, 24*, 997–1016.

Gilbert, B., Datta, A., Casey, J. A., & Ogburn, E. L. (2024). A causal inference framework for spatial confounding. *arXiv:2112.14946v12*.

Gilbert, B., Ogburn, E. L., & Datta, A. (2025). Consistency of common spatial estimators under spatial confounding. *Biometrika, 112*, asae070.

Hanks, E. M., Schliep, E. M., Hooten, M. B., & Hoeting, J. A. (2015). Restricted spatial regression in practice: Geostatistical models, confounding, and robustness under model misspecification. *Environmetrics, 26*(4), 243–254.

Hodges, J. S., & Reich, B. J. (2010). Adding spatially-correlated errors can mess up the fixed effect you love. *The American Statistician, 64*(4), 325–334.

Hughes, J., & Haran, M. (2013). Dimension reduction and alleviation of confounding for spatial generalized linear mixed models. *Journal of the Royal Statistical Society, Series B: Statistical Methodology, 75*(1), 139–159.

Hui, F. K. C., Vu, Q., & Hooten, M. B. (2024). Spatial confounding in joint species distribution models. *Methods in Ecology and Evolution, 15*, 1906–1921.

Inventaire National du Patrimoine Naturel. (2017). *Données issues du dispositif de Biosurveillance des Retombées Atmosphériques Métalliques par les Mousses (BRAMM)*. UMS PatriNat (OFB-CNRS-MNHN). Occurrence dataset. https://doi.org/10.15468/c2bqic

Jin, X., Banerjee, S., & Carlin, B. P. (2007). Order-free co-regionalized areal data models with application to multiple-disease mapping. *Journal of the Royal Statistical Society, Series B: Statistical Methodology, 69*(5), 817–838.
<figure><img src="image_19.png" /></figure>

Khan, K., & Calder, C. A. (2022). Restricted spatial regression methods: Implications for inference. *Journal of the American Statistical Association*, 117(537), 482-494.

Lamouroux, J., Geffroy, A., Leblond, S., Meyer, C., & Albert, I. (2025). Zenodo Repository. https://doi.org/10.5281/zenodo.15748176

Leblond, S., Laffray, X., Galsomies, L., & Gombert-Courvoisier, S. (2011). Le dispositif bramm: un outil de biosurveillance de la qualité de l'air. *Pollution Atmospherique*, 49-53.

Lequy, E., Leblond, S., Siemiatycki, J., Meyer, C., Vienneau, D., de Hoogh, K., Zins, M., Goldberg, M., & Jacquemin, B. (2023). Long-term exposure to airborne metals and risk of cancer in the french cohort gazel. *Environment International*, 177, 107999.

Lindgren, F., & Rue, H. (2015). Bayesian spatial modelling with r-inla. *Journal of Statistical Software*, 63(19), 1-25.

Lindgren, F., Rue, H., & Lindström, J. (2011). An explicit link between gaussian fields and gaussian markov random fields: The stochastic partial differential equation approach. *Journal of the Royal Statistical Society, Series B: Statistical Methodology*, 73(4), 423-498.

Marques, I., & Kneib, T. (2022). Discussion on "spatial+: A novel approach to spatial confounding" by Dupont, Emiko, Wood, Simon N., and Augustin, Nicole H. *Biometrics*, 78(4), 1295-1299.

Marques, I., & Wiemann, P. F. V. (2023). Bayesian spatial+: A joint model perspective. *arXiv preprint*, arXiv:2309.05496. https://arxiv.org/abs/2309.05496

Martinez-Beneito, M. A. (2013). A general modelling framework for multivariate disease mapping. *Biometrika*, 100(3), 539-553.

McCullagh, P. (1989). *Generalized linear models* (2nd ed.). Routledge.

Paciorek, C. J. (2010). The importance of scale for spatial-confounding bias and precision of spatial regression estimators. *Statistical Science: A Review Journal of the Institute of Mathematical Statistics*, 25(1), 107-125.

Page, G. L., Liu, Y., He, Z., & Sun, D. (2017). Estimation and prediction in the presence of spatial confounding for spatial linear models. *Scandinavian Journal of Statistics*, 44(3), 780-797.

Prates, M. O., Dey, D. K., Willig, M. R., & Yan, J. (2015). Transformed gaussian markov random fields and spatial modeling of species abundance. *Spatial Statistics*, 14, 382-399.

Reich, B. J., Hodges, J. S., & Zadnik, V. (2006). Effects of residual smoothing on the posterior of the fixed effects in disease-mapping models. *Biometrics*, 62(4), 1197-1206.

Riebler, A., Sørbye, S. H., Simpson, D., & Rue, H. (2016). An intuitive bayesian spatial model for disease mapping that accounts for scaling. *Statistical Methods in Medical Research*, 25(4), 1145-1165.

Rue, H., & Held, L. (2005). *Gaussian Markov random fields: Theory and applications*. CRC Press.

Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. *Journal of the Royal Statistical Society, Series B: Statistical Methodology*, 71(2), 319-392.

Schmidt, A. M. (2022). Discussion on "Spatial+: A novel approach to spatial confounding" by Emiko Dupont, Simon N. Wood, and Nicole H. Augustin. *Biometrics*, 78(4), 1300-1304.

Schröder, W., & Nickel, S. (2019). Moss species-specific accumulation of atmospheric deposition? *Environmental Sciences Europe*, 31, 1-18.

Sørbye, S. H., Illian, J. B., Simpson, D. P., Burslem, D., & Rue, H. (2019). Careful prior specification avoids incautious inference for log-Gaussian Cox point processes. *Journal of the Royal Statistical Society: Series C: Applied Statistics*, 31, 543-564.

Thaden, H., & Kneib, T. (2018). Structural equation models for dealing with spatial confounding. *The American Statistician*, 72(3), 239-252.

Tobler, W. R. (1970). A computer movie simulating urban growth in the detroit region. *Economic Geography*, 46(Sup1), 234-240.

Travnikov, O., & Ilyin, I. (2005). Regional model MSCE-HM of heavy metal transboundary air pollution in Europe. Moscow, Meteorological Synthesizing Centre - East (EMEP/MSC-E Technical Report 6/2005).

Urdangarin, A., Goicoa, T., Kneib, T., & Ugarte, M. D. (2024). A simplified spatial+ approach to mitigate spatial confounding in multivariate spatial areal models. *Spatial Statistics*, 59, 100804.

Urdangarin, A., Goicoa, T., & Ugarte, M. D. (2023). Evaluating recent methods to overcome spatial confounding. *Revista Matemática Complutense*, 36(2), 333-360.

Wahba, G. (1990). *Spline models for observational data*. SIAM.

Zaccardi, C., Valentini, P., Ippoliti, L., & Schmidt, A. M. (2025). Regularized principal spline functions to mitigate spatial confounding. *Biometrics*, 81(2), ujaf076. https://doi.org/10.1093/biomtc/ujaf076

**How to cite this article:** Lamouroux, J., Geffroy, A., Leblond, S., Meyer, C., & Albert, I. (2025). Addressing spatial confounding in geostatistical regression models: An R-INLA approach. *Methods in Ecology and Evolution*, 16, 2082-2097. https://doi.org/10.1111/2041-210X.70106