Lecture 6

Modelling geostatistical data

Sara Martino

sara.martino@ntnu.no

Dept. of Mathematical Science, NTNU

Janine Illian

Janine.Illian@glasgow.ac.uk

University of Glasgow

Jafet Belmont

Jafet.Belmont@glasgow.ac.uk

University of Glasgow

Geostatistical modelling

Geostatistical data

In many ecological and environmental settings, measurements are taken from fixed sampling units, aiming to quantify spatial variation and interpolate values at unobserved sites.

Geostatistical data are the most common form of spatial data found in environmental and ecological settings.
We regularly take measurements of an environmental variable of interest at a set of fixed locations.
This could be data from samples taken across a region (eg., water depth in a lake) or from monitoring stations as part of a network (e.g., air pollution).
In each of these cases, our goal is to estimate the value of our variable across the entire space.

Understanding our region

Let \(D\) be our two-dimensional region of interest.
In principle, there are infinitely many locations within \(D\), each of which can be represented by mathematical coordinates (e.g. latitude and longitude).
We can identify any individual location as \(\mathbf{s}_i = (x_i, y_i)\), where \(x_i\) and \(y_i\) are their coordinates.
We can treat our variable of interest as a random variable, \(Z\) which can be observed at any location as \(Z(\mathbf{s}_i)\).

Geostatistical process

A geostatistical process can therefore be written as: \[\{Z(\mathbf{s}); \mathbf{s} \in D\}\]
In practice, data are observed at a finite number of locations, \(m\), and can be denoted as: \[z = \{z(\mathbf{s}_1), \ldots z(\mathbf{s}_m) \}\]
We have observed our data at \(m\) locations, but often want to predict this process at a set of unknown locations.
For example, what is the value of \(z(\mathbf{s}_0)\), where \(\mathbf{s}_0\) is an unobserved site?

Spatial autocorrelation

The key challenge in modelling geostatistical data is understanding correlation.
Typically observations close together in space will be more similar than those which are further apart.
Spatial correlation is usually driven by some unmeasured confounding variable(s) - for example, air pollution is spatially correlated because nearby areas tend to experience similar traffic levels.
It is important that we account for these correlations in our analysis - failing to do so will lead to poor inference.

Spatial Modelling

What We Observed and What We Did Not

We treat the observed process of interest as being measured with error

\[ (\text{observed value})_i = (\text{true value at location } i) + (\text{error})_i \]

alternatively

\[ y_i = Z(\mathbf{s}_i) + \varepsilon_i \]

When geostatistical data are considered, we can often assume that there is a spatially continuous variable underlying the observations that can be modelled using a random field.

we have a process that is occurring everywhere in space \(\rightarrow\) natural to try to model it using some sort of function (of space)
a random field is a random function that generates smooth surfaces
this is hard
we typically make our lives easier by making everything Gaussian

Gaussian Random Fields

A Gaussian random field (GRF) is a collection of random variables, where observations occur in a continuous domain, and where every finite collection of random variables has a multivariate normal distribution

\[ \mathbf{z} = (z(\mathbf{s}_1),\ldots,z(\mathbf{s}_m)) \sim N(\mu(\mathbf{s}_1),\ldots,\mu(\mathbf{s}_m),\Sigma), \]

where \(\Sigma_{ij} = \mathrm{Cov}(z(\mathbf{s}_i),z(\mathbf{s}_j))\) is a dense \(m \times m\) matrix.

This is actually quite tricky: \(\Sigma\) will need to depend on the set of observation sites and has to behave well (be “positive definite”)
Use a covariance function \(C_z(\cdot,\cdot)\) that depends on the distance (\(\Sigma_{ij} = C_z(\mathbf{s}_i,\mathbf{s}_j)\)) between two points and that
- has no negative variances
- is symmetric
- is decreasing, with maximum at distance = 0
\(C_z(\mathbf{s}_i, \mathbf{s}_j)\) measures the strength of the linear dependence between \(Z(\mathbf{s}_i)\) and \(Z(\mathbf{s}_j)\).
\(C_z(\mathbf{s}_i, \mathbf{s}_j) = \mathrm{Var}(Z(\mathbf{s}_i))\) for \(i = j\).

Building the model

The first step in defining a model for a random field in a hierarchical framework is to identify a probability distribution for the observations available at \(m\) sampled spatial locations and represented by the vector \(\mathbf{y} = y_1,\ldots,y_m\).

For example, if we assume our observations follow a Gaussian distribution then

\[ \begin{aligned} Y_i &\sim N(\mu_i,\tau_e^{-1})\\ \eta_i &=\mu_i = \beta_0 + \ldots + Z(\mathbf{s}_i) \end{aligned} \]

\(\tau_e^{-1} = \sigma^2_e\) represents the variance of the zero-mean measurement error
The response mean \(\mu_i\), which coincides with the linear predictor \(\eta_i\), is defined based on:
- the intercept \(\beta_0\) and any additional covariates
- the realization of the latent (unobservable) GF \(Z(\mathbf{s}) \sim \mathrm{MVN}(0,\Sigma)\), which accounts for the spatial correlation through \(\Sigma = C_z(\cdot,\cdot)\).

The Matérn Field

A commonly used covariance function is the Matérn covariance function. The covariance of two points which are a distance \(h\) apart is:

\[ \Sigma =C_{\nu}(h) = \sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{\sqrt{2\nu} h}{\rho} \right)^{\nu} K_{\nu} \left( \frac{\sqrt{2\nu} h}{\rho} \right) \]

\(\Gamma(\cdot)\) is the gamma function
\(K_{\nu}(\cdot)\) is the modified Bessel function of the second kind.
Parameters \(\sigma^2\), \(\rho\) and \(\nu\) are non-negative values of the covariance function.
- \(\sigma^2\) is the spatially structured variance component
- \(\rho\) is the range of the spatial process
- \(\nu\) controls smoothness of the spatial process.

Big n problem!

The disadvantage of the modelling approach involving the spatial covariance function is known as “big n problem” and concerns the computational costs required for algebra operations with dense covariance matrices (such as \(\Sigma\)).

In particular dense matrix operations scale cubically with the matrix size, given by the number of locations where the process is observed. A computationally effective alternative is given by the stochastic partial differential equation (SPDE) approach

The SPDE approach

a different approach to GRF models:

forget about the covariance function!
express spatial structure through an SPDE (stochastic partial differential equation):
- a simple differential equation describes change in time
- a partial differential equation describes change in space (in more than 1 direction!)
- a stochastic partial differential equation makes this random,i.e. not deterministic
need to approximate spatial structures in a clever way.

The SPDE approach

The SPDE approach in a nutshell

The stochastic partial differential equation

We define a (Matérn) GRF as the solution of a stochastic partial differential equation (SPDE)

\[ (\kappa^2-\Delta)^{\alpha/2}Z(t) = W(t) \]

What is this?

\(W(t)\) is random noise
\(\omega(t)\) is the smooth process we want
\((\kappa^2-\Delta)^{\alpha/2}\) is an operator that “smoothes” the white noise.
\(\kappa\) and \(\alpha\) are parameters

Recall the 1-D SPDE

Imagine a guitar string stretched from left to right.

Now imagine someone randomly taps along it at many locations:
- Each tap is independent
- Some taps are strong, some weak
- There is no coordination between taps
- A tap at one location tells you nothing about a tap nearby

This is pure randomness. But a real string does not behave like this

Recall the 1-D SPDE

Imagine a guitar string stretched from left to right.

The string has tension and stiffness
- the tension spreads each tap to nearby points
- Sharp jumps are softened. \((\kappa^2-\Delta)^{\color{red}{\alpha}/2}Z(t) = W(t)\)

\(\Delta\) measures the local curvature
\(\kappa\) controls how far randomness propagates
stronger tension (large \(\kappa\)) = stronger pull

\[(\color{red}{\kappa^2}-\Delta)Z(t) = W(t), ~~~\text{for } \alpha=2\]

Measure the bending at \(z_i\) by taking the with respect the average of the neighbours (i.e. \(z_i - (z_{i-1}+z_i/2) = z_{i-1} + 2z_i +z_{i+1}\)).

if \(= 0 \rightarrow\) straight
if the diff is large \(\rightarrow\) sharply bent
the sign tells you which way it bends

The curvature need to divide by the spacing \(h\) between points (assuming equal spacing), i.e., \((z_{i-1} + 2z_i +z_{i+1}) \times h^{-2}\).

In a continuous space we measured this by taking the second derivative, i.e. the Laplacian operator

As \(h\to 0\) we measure the rate of change of slope, i.e. curvature as:

\[(z_{i-1} + 2z_i +z_{i+1}) \times h^{-2} \to \frac{\partial z}{\partial t^2} = \Delta z\] Where \(z(x)\) is the displacement of the string at position \(t\).

Tension (\(\kappa\)): resists large displacement If a tap pushes part of the string far up or down:

tension pulls it back
stronger tension = stronger pull

\[ \begin{aligned} (\color{red}{\kappa^2}-\Delta)Z(t) &= W(t) \\ \color{red}{\kappa^2}Z(t) - \Delta Z(t) &= W(t)\\ \text{if } Z(t) \text{ is large then }& \color{tomato}{\kappa^2 Z(t)} \text{ becomes large too} \\ \text{large pull}-\text{back bending}&=\text{random tap} \end{aligned} \] In this case, To balance the equation, \(Z(t)\) must be smaller.

Recall the 1-D SPDE

Imagine a guitar string stretched from left to right.

The rope has tension and stiffness
- the tension spreads each tap to nearby points
- Sharp jumps are softened.

Stiffness = controls how smooth the field becomes (\(\alpha\))
Larger \(\alpha\) smoother the process will be

\[ (\kappa^2-\Delta)^{\color{red}{\alpha}/2}Z(t) = W(t) \]

Solving the SPDE

Ok…but we still need to solve the SPDE to find \(Z(t)\)!

Now we need to discretize the domain into T points (we cannot represent continuous space in the computer!)

We represent our solution as

\[ Z(t) = \sum_{i = 1}^T\psi_i(t)w_i, \]

where

\(\psi_i(t)\) are (known) basis functions for nodes \(i=1,\ldots,T\)
- \(\psi_i(t_i)= 1\)
- \(\psi_i(t_j) = 0 ~~\forall~~i \neq j\)
- Linear between neighbouring nodes

“solving the SPDE” means finding a random function \(Z(t)\) such that the equality \((\kappa^2-\Delta)^{\alpha/2}Z(t) = W(t)\) holds in distribution.
in other words, finding \(Z(s)\) such that when the operator \((\kappa^2-\Delta)^{\alpha/2}\) is applied to it, it produces white noise
OR, The Matérn field is the Gaussian field \(Z(s)\) produced by smoothing white noise \(W(t)\) through the inverse operator \((\kappa^2-\Delta)^{\alpha/2}\).
The Matérn field is infinite-dimensional so we approximate it using a finite basis expansion
- we want a simple local basis
- each basis only affects neighbors
- so computations become sparse
Each weight \(w_i\) is basically the value of the field at node \(i\)
The basis functions are:
- Height = 1 at node
- Zero at neighboring nodes
- Linear in between
A hat function is not the field. It is a building block for the field.
the field between nodes is interpolated linearly so the vector of weights fully defines the approximate field.
Since the SPDE is linear and the white noise is Gaussian, the solution is Gaussian

Solving the SPDE

Ok…but we still need to solve the SPDE to find \(Z(t)\)!

Now we need to discretize the domain into T points (we cannot represent continuous space in the computer!)

We represent our solution as

\[ Z(t) = \sum_{i = 1}^T\psi_i(t)w_i \]

Where

\(\psi_i(t)\) are (known) basis functions for nodes \(i=1,\ldots,T\)
\(w_i\) are (unknown) weights
- the field value \(Z(s)\) is a linear interpolation between the two neighboring weights

“solving the SPDE” means finding a random function \(Z(t)\) such that the equality \((\kappa^2-\Delta)^{\alpha/2}Z(t) = W(t)\) holds in distribution.
in other words, finding \(Z(s)\) such that when the operator \((\kappa^2-\Delta)^{\alpha/2}\) is applied to it, it produces white noise
OR, The Matérn field is the Gaussian field \(Z(s)\) produced by smoothing white noise \(W(t)\) through the inverse operator \((\kappa^2-\Delta)^{\alpha/2}\).
The Matérn field is infinite-dimensional so we approximate it using a finite basis expansion
- we want a simple local basis
- each basis only affects neighbours
- so computations become sparse
Each weight \(w_i\) is basically the value of the field at node \(i\)
The basis functions are:
- Height = 1 at node
- Zero at neighboring nodes
- Linear in between
A hat function is not the field. It is a building block for the field.
the field between nodes is interpolated linearly so the vector of weights fully defines the approximate field.
Since the SPDE is linear and the white noise is Gaussian, the solution is Gaussian

Solving the SPDE

Ok…but we still need to solve the SPDE to find \(Z(t)\)!

Now we need to discretize the domain into T points (we cannot represent continuous space in the computer!)

We represent our solution as

\[ Z(t) = \sum_{i = 1}^T\psi_i(t)w_i \]

Where

\(\psi_i(t)\) are (known) basis functions for nodes \(i=1,\ldots,T\)
\(w_i\) are (unknown) weights
This solution is then approximated using a finite combination of piece-wise linear basis functions.
The solution is completely defined by a Gaussian vector of weights with zero mean and a sparse precision matrix.

The SPDE approach on 2D

Now we approximate the GRF using a triangulated mesh.

The SPDE approach represents the continuous spatial process as a continuously indexed Gaussian Markov Random Field (GMRF)

We construct an appropriate lower-resolution approximation of the surface by sampling it in a set of well designed points and constructing a piece-wise linear interpolant.

The SPDE approach on 2D

Now we approximate the GRF using a triangulated mesh.

The SPDE approach represents the continuous spatial process as a continuously indexed Gaussian Markov Random Field (GMRF)

We construct an appropriate lower-resolution approximation of the surface by sampling it in a set of well designed points and constructing a piece-wise linear interpolant.
Note that \(\nu = \alpha - d/2\). For \(\alpha=2 \Rightarrow \nu= 1\) since \(d=2\) we have that:

Note

\[ \begin{aligned} Z(s) &= \sum_{i = 1}^K\psi_i(s)w_i \\ \mathbf{w} &\sim N(\mathbf{0},Q^{-1}) \leftarrow \text{GMRF}\\ Q^{-1} &= \tau^2(\kappa^4 \mathbf{C} + 2\kappa^2 \mathbf{G}+\mathbf{G}\mathbf{C}^{-1}\mathbf{G}) \end{aligned} \]

\(\mathbf{C}\) is diagonal with entries \(C_{ii} =\int \psi_i(s)\mathrm{d}s\) and measures how much of the domain each basis function covers.
\(G_{ij} = \int \nabla \psi_i(s) \nabla \psi_j(s) \mathrm{d}s\) reflects the connectivity of the mesh nodes.
because each basis function overlaps only with nearby ones, the resulting precision matrix is sparse, meaning each coefficient depends directly only on its neighbors

In summary

The continuous Matérn GRF is the solution of a SPDE and is represented as

\[ Z(s) = \sum_{i = 1}^K\psi_i(s)w_i \]

The weights vector \(\mathbf{w} = (w_1,\dots,w_K)\) is Gaussian with a sparse precision matrix \(\longrightarrow\) Computational convenience
The field has two parameters
- The range \(\rho\)
- The marginal variance \(\sigma^2\)
These parameters are linked to the parameters of the SPDE
We need to assign prior to them

Penalized Complexity (PC) priors

Penalized Complexity (PC) priors proposed by Simpson et al. (2017) allow us to control the amount of spatial smoothing and avoid overfitting.

PC priors shrink the model towards a simpler baseline unless the data provide strong evidence for a more complex structure.
To define the prior for the marginal precision \(\sigma^{-2}\) and the range parameter \(\rho\), we use the probability statements:
- Define the prior for the range \(\text{Prob}(\rho<\rho_0) = p_{\rho}\)
- Define the prior for the range \(\text{Prob}(\sigma>\sigma_0) = p_{\sigma}\)

Learning about the SPDE approach

D. Miller, R. Glennie and A. Seaton. (2020). Understanding the Stochastic Partial Differential Equation Approach to Smoothing. Journal of Agricultural, Biological and Environmental Statistics, 25(1), 1-16.
F. Lindgren, H. Rue, and J. Lindström. An explicit link between Gaussian fields and Gaussian Markov random fields: The SPDE approach (with discussion). In: Journal of the Royal Statistical Society, Series B 73.4 (2011), pp. 423–498.
H. Bakka, H. Rue, G. A. Fuglstad, A. Riebler, D. Bolin, J. Illian, E. Krainski, D. Simpson, and F. Lindgren. Spatial modelling with R-INLA: A review. In: WIREs Computational Statistics 10:e1443.6 (2018). (Invited extended review). DOI: 10.1002/wics.1443.
E.T. Krainski, V. Gómez-Rubio, H. Bakka, A. Lenzi, D. Castro-Camilio, D. Simpson, F. Lindgren, and H. Rue. Advanced Spatial Modeling with Stochastic Partial Differential Equations using R and INLA. Github version . CRC press, Dec. 20

Modelling Pacific Cod Distribution

In the next example, we will explore data on the Pacific Cod (Gadus macrocephalus) from a trawl survey in Queen Charlotte Sound.

The dataset contains the presence/absence of Pacific cod in the area swept for a given survey in 2003 as well as depth covariate information.

The Model

Stage 1 Model for the response \[ y(s)|\eta(s)\sim\text{Binomial}(p(s)) \]
Stage 2 Latent field model \[ \eta(s) = \mathrm{logit}(p(s)) = \beta_0 + \beta_1 \,\log \text{depth}(s)+ Z(s) \]
Stage 3 Hyperparameters

The Model

Stage 1 Model for the response \[ y(s)|\eta(s)\sim\text{Binomial}(p(s)) \]
Stage 2 Latent field model \[ \eta(s) = \mathrm{logit}(p(s)) = \beta_0 +\beta_1 \,\log \text{depth}(s) + Z(s) \]
- A global intercept \(\beta_0\)
- A covariate effect \(\beta_1\)
- A Gaussian field \(Z(s)\)
Stage 3 Hyperparameters

The Model

Stage 1 Model for the response \[ y(s)|\eta(s)\sim\text{Binomial}(p(s)) \]
Stage 2 Latent field model \[ \eta(s) = \mathrm{logit}(p(s)) = \beta_0 + \beta_1 \,\log \text{depth}(s) + Z(s) \]
Stage 3 Hyperparameters
- Range and sd in the Gaussian field \(\sigma_{\omega}, \tau_{\omega}\)

Step 1: Define the SPDE representation: The mesh

First, we need to create the mesh used to approximate the random field.

max.edge for maximum triangle edge lengths
offset for inner and outer extensions (to prevent edge effects)
cutoff to avoid overly small triangles in clustered areas

Step 1: Define the SPDE representation: The mesh

All random field models need to be discretised for practical calculations.
The SPDE models were developed to provide a consistent model definition across a range of discretisations.
We use finite element methods with local, piecewise linear basis functions defined on a triangulation of a region of space containing the domain of interest.
Deviation from stationarity is generated near the boundary of the region.
The choice of region and choice of triangulation affects the numerical accuracy.

Step 1: Define the SPDE representation: The mesh

If the mesh is too fine \(\rightarrow\) heavy computation
If the mesh is to coarse \(\rightarrow\) not accurate enough

Step 1: Define the SPDE representation: The mesh

Some guidelines

Create triangulation meshes with fm_mesh_2d():
edge length should be around a third to a tenth of the spatial range
Move undesired boundary effects away from the domain of interest by extending to a smooth external boundary:
Use a coarser resolution in the extension to reduce computational cost (max.edge=c(inner, outer)), i.e., add extra, larger triangles around the border

Step 1: Define the SPDE representation: The mesh

Use a fine resolution (subject to available computational resources) for the domain of interest (inner correlation range) and avoid small edges ,i.e., filter out small input point clusters (0 \(<\) `cutoff \(<\) inner)
Coastlines and similar can be added to the domain specification in fm_mesh_2d() through the boundary argument.
simplify the border

Step 1: Define the SPDE representation: The SPDE

We use the inla.spde2.pcmatern to define the SPDE model using PC priors through the following probability statements

\(P(\rho < 100) = 0.5\)
\(P(\sigma > 1) = 0.5\)

Step 2: Define the model components

The Model

\[ \begin{aligned} y(s)|\eta(s) & \sim\text{Binom}(1, p(s))\\ \eta(s) & = \color{#FF6B6B}{\boxed{\beta_0}} + \color{#FF6B6B}{\boxed{ \beta_1 \,\log \text{depth}(s) }} + \color{#FF6B6B}{\boxed{ Z(s)}}\\ \end{aligned} \]

depth is a column in our sf data
geometry column allow us to use spatial objects as input

pcod_sf %>% select(depth, present) %>% print(n = 3)

Simple feature collection with 232 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 343.0617 ymin: 5639.857 xmax: 576.9152 ymax: 5837.211
Projected CRS: +proj=utm +zone=9 +datum=WGS84 +no_defs +type=crs +units=km
# A tibble: 232 × 3
  depth present            geometry
  <dbl>   <dbl>        <POINT [km]>
1   201       1 (446.4752 5793.426)
2   212       1 (446.4594 5800.136)
3   220       0 (448.5987 5801.687)
# ℹ 229 more rows

The code

# define model component
cmp = ~ -1 + Intercept(1) +  depth_eff(log(depth), model='linear') + 
  space(geometry, model = spde_model)

# define model predictor
eta = present ~ Intercept + depth_eff + space

# build the observation model
lik = bru_obs(formula = eta,
              data = pcod_sf,
              family = "binomial")

# fit the model
fit = bru(cmp, lik)

Step 2: Define the model components

The Model \[ \begin{aligned} y(s)|\eta(s) & \sim\text{Binom}(1, p(s))\\ \eta(s) & = \beta_0 + \beta_1 \,\log \color{#FF6B6B}{\boxed{\text{depth}(s) } } + Z(s)\\ \end{aligned} \]

here depth_r is a raster and depth a layer (be careful with NA’s).
.data. is a keyword referring to whatever data object is passed to bru_obs() or bru()
eval_spatial() extracts the value of the spatial covariate at the locations in the data

Spatial covariates

~ my_sp_effect(
  main = a_spatial_object,
  model = "linear"
)

~ my_sp_effect(
  main = eval_spatial(a_spatial_object, .data.),
  model = "linear"
)

The code

# define model component
cmp = ~ -1 + Intercept(1) +  depth_eff(log(depth_r$depth), model='linear') + 
  space(geometry, model = spde_model)

# define model predictor
eta = present ~ Intercept + depth_eff + space

# build the observation model
lik = bru_obs(formula = eta,
              data = pcod_sf,
              family = "binomial")

# fit the model
fit = bru(cmp, lik)

Step 3: Define the linear predictor

The Model \[ \begin{aligned} y(s)|\eta(s) & \sim\text{Binom}(1, p(s))\\ \color{#FF6B6B}{\boxed{\eta(s)}} & = \color{#FF6B6B}{\boxed{\beta_0 + \beta_1 \,\log \text{depth}(s) + Z(s)}}\\ \end{aligned} \]

pcod_sf %>% select(depth, present) %>% print(n = 3)

Simple feature collection with 232 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 343.0617 ymin: 5639.857 xmax: 576.9152 ymax: 5837.211
Projected CRS: +proj=utm +zone=9 +datum=WGS84 +no_defs +type=crs +units=km
# A tibble: 232 × 3
  depth present            geometry
  <dbl>   <dbl>        <POINT [km]>
1   201       1 (446.4752 5793.426)
2   212       1 (446.4594 5800.136)
3   220       0 (448.5987 5801.687)
# ℹ 229 more rows

The code

# define model component
cmp = ~ -1 + Intercept(1) +  depth_eff(log(depth), model='linear') + 
  space(geometry, model = spde_model)

# define model predictor
eta = present ~ Intercept + depth_eff + space

# build the observation model
lik = bru_obs(formula = eta,
              data = pcod_sf,
              family = "binomial")

# fit the model
fit = bru(cmp, lik)

Step 4: Define the observational model

The Model

\[ \begin{aligned} \color{#FF6B6B}{\boxed{y(s)|\eta(s)}} & \sim \color{#FF6B6B}{\boxed{\text{Binom}(1, p(s))}}\\ \eta(s) & = \beta_0 + \beta_1 \,\log \text{depth}(s) + Z(s)\\ \end{aligned} \]

pcod_sf %>% select(depth, present) %>% print(n = 3)

Simple feature collection with 232 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 343.0617 ymin: 5639.857 xmax: 576.9152 ymax: 5837.211
Projected CRS: +proj=utm +zone=9 +datum=WGS84 +no_defs +type=crs +units=km
# A tibble: 232 × 3
  depth present            geometry
  <dbl>   <dbl>        <POINT [km]>
1   201       1 (446.4752 5793.426)
2   212       1 (446.4594 5800.136)
3   220       0 (448.5987 5801.687)
# ℹ 229 more rows

The code

# define model component
cmp = ~ -1 + Intercept(1) +  depth_eff(log(depth), model='linear') + 
  space(geometry, model = spde_model)

# define model predictor
eta = present ~ Intercept + depth_eff + space

# build the observation model
lik = bru_obs(formula = eta,
              data = pcod_sf,
              family = "binomial")

# fit the model
fit = bru(cmp, lik)

Step 5: Collect Results

We can plot the posterior density of the model components via plot(model, "model component") and also look at the correlation/covariance function as follows:

Code

Inter_post <- plot(fit, "Intercept")
depth_post <- plot(fit, "depth_eff")
corplot <- plot(spde.posterior(fit, "space", what = "matern.correlation"))
covplot <- plot(spde.posterior(fit, "space", what = "matern.covariance"))

Inter_post + depth_post + corplot + covplot + plot_layout(ncol=2)

Spatial predictions

How do we predict at unsampled locations?

In geostatistical applications, the main interest resides in the spatial prediction of the spatial latent field or of the response variable in new locations

Suppose we observe a spatial process \({Z(s): s \in \mathcal{D}}\) at locations \(s_1,\dots,s_n\).
Our goal: predict the variable of interest at an unobserved location \(s_0 \in \mathcal{D}\).

given the data \(y = (y_1,\dots,y_n)\), what can we say about \(Z(s_0)\)?
Rather than a single guess, we want a full uncertainty-aware prediction.
In a Bayesian setting, prediction is a probabilistic task.

Posterior predictive density

The key lies in the posterior predictive distribution

\[ \pi(\tilde{Y} \mid y) = \int \pi(\tilde{Y} \mid \Theta, y)\, \pi(\Theta \mid y)\, d\Theta, \] where \(\Theta\) denotes all latent components and hyperparameters.

\(\pi(\tilde{Y} \mid \Theta, y)\) depends on the task:
- extrapolation (e.g. forecasting): \(\pi(Y_{n+1} \mid \Theta, y_n)\),
- interporlation: \(\pi(Y_i \mid y_{i-1}, y_{i+1}, \Theta)\),
Spatial prediction fits naturally into this framework:
- \(\tilde{Y}\) may represent \(Z(s_0)\), \(\eta_0\), or the response at \(y(s_0)\),
- conditioning reflects the assumed spatial dependence.
INLA approximates \(\pi(\Theta \mid y)\) efficiently, enabling full uncertainty propagation when predicting over \(s_0 \in \mathcal{D}\).

Posterior predictive density

Define a grid of points where we want to predict. fm_pixel() creates a regular grid of points covering the mesh (we will see in the practical how your the raster can be used instead)
Append the values of the covariate(s) (if any) at the prediction locations
inlabru will create a projector matrix \(A\) linking the latent Gaussian field to the prediction locations via the predict function.
We can do predictions for \(Z(s), \eta(s), g^{-1}(\eta(s))\)

pxl = fm_pixels(mesh)

# 2. extract raster values at the point locations
vals <- terra::extract(depth_r$depth, pxl, ID = FALSE)

# 3. bind the extracted values onto the sf object
pred_data <- cbind(pxl, vals) %>%  filter(!is.na(depth))


preds = predict(fit, pred_data, 
      ~data.frame(spde = space,
                  eta = Intercept + 
                        depth_eff +
                        space,
                  p_s = plogis(Intercept + 
                               depth_eff +
                               space)))

Posterior predictive density

Outputs include posterior means, variances, and credible intervals.

spde_pred <- ggplot() +
   gg(preds$spde, aes(fill = q0.5),
      geom = "tile") + 
   scale_fill_scico(name= "median SPDE" ) 

eta_pred <- ggplot() +
   gg(preds$eta, aes(fill = sd),
      geom = "tile") + 
   scale_fill_scico(name=expression(eta~"sd"),
                    palette = "lapaz") 

p_pred <- ggplot() +
   gg(preds$p_s, aes(fill = mean),
      geom = "tile") + 
   scale_fill_scico(palette = "roma") 

spde_pred +  eta_pred + p_pred + plot_layout(ncol=1)

Take-Home Message

Gaussian random fields provide components that reflect spatial structure in continuous space
the SPDE approach allows us to approximate these efficiently and flexibly
the models are still GMRFs ⇝ computationally efficient representation
pc-priors are also relevant for spatial models, e.g. to avoid overfitting

References

Bakka, Haakon, Jarno Vanhatalo, Janine B Illian, Daniel Simpson, and Håvard Rue. 2019. “Non-Stationary Gaussian Models with Physical Barriers.” Spatial Statistics 29: 268–88.

Bryce, Erin, Luigi Lombardo, Cees Van Westen, Hakan Tanyas, and Daniela Castro-Camilo. 2022. “Unified Landslide Hazard Assessment Using Hurdle Models: A Case Study in the Island of Dominica.” Stochastic Environmental Research and Risk Assessment 36 (8): 2071–84.

Cameletti, Michela, Finn Lindgren, Daniel Simpson, and Håvard Rue. 2013. “Spatio-Temporal Modeling of Particulate Matter Concentration Through the SPDE Approach.” AStA Advances in Statistical Analysis 97 (2): 109–31.

Fichera, Andrew, Rachel King, Jarrod Kath, David Cobon, and Kathryn Reardon-Smith. 2023. “Spatial Modelling of Agro-Ecologically Significant Grassland Species Using the INLA-SPDE Approach.” Scientific Reports 13 (1): 4972.

Lamouroux, Jérémy, Alizée Geffroy, Sébastien Leblond, Caroline Meyer, and Isabelle Albert. 2025. “Addressing Spatial Confounding in Geostatistical Regression Models: An r-INLA Approach.” Methods in Ecology and Evolution 16 (9): 2082–97.

Lindgren, Finn, David Bolin, and Håvard Rue. 2022. “The SPDE Approach for Gaussian and Non-Gaussian Fields: 10 Years and Still Running.” Spatial Statistics 50: 100599.

Lindgren, Finn, Johan Lindström, and Håvard Rue. 2010. An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields; the Spde Approach.

Miller, David L, Richard Glennie, and Andrew E Seaton. 2020. “Understanding the Stochastic Partial Differential Equation Approach to Smoothing: DL Miller Et Al.” Journal of Agricultural, Biological and Environmental Statistics 25 (1): 1–16.

Moraga, Paula, Christopher Dean, Joshua Inoue, Piotr Morawiecki, Shahzeb Raja Noureen, and Fengpei Wang. 2021. “Bayesian Spatial Modelling of Geostatistical Data Using INLA and SPDE Methods: A Case Study Predicting Malaria Risk in Mozambique.” Spatial and Spatio-Temporal Epidemiology 39: 100440.

Simpson, Daniel, Håvard Rue, Andrea Riebler, Thiago G. Martins, and Sigrunn H. Sørbye. 2017. “Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors.” Statistical Science 32 (1): 1–28. https://doi.org/10.1214/16-STS576.