Space is inherent to all ecological processes, influencing dynamics such as migration, dispersal, and species interactions. A primary goal in ecology is to understand how these processes shape species distributions and dynamics across space.
often inaccessible to practitioners (literature aimed at statisticians)
methodology not always linked to applications
models too simple to reflect real-life data
difficult to apply without expertise in spatial stats and computational statistics
we shall see that inlabru can help with this.
Spatial data structures
We can distinguish three types of spatial data structures
Areal data
Map of bird conservation regions (BCRs) showing the proportion of bird species within each region showing a declining trend
Geostatistical data
Scotland river temperature monitoring network
Point pattern data
Occurrence records of four ungulate species in Tibet,
Spatial data structures
We can distinguish three types of spatial data structures
Areal data
In areal data our measurements are summarised across a set of discrete, non-overlapping spatial units.
Map of bird conservation regions (BCRs) showing the proportion of bird species within each region showing a declining trend
Geostatistical data
Scotland river temperature monitoring network
Point pattern data
Occurrence records of four ungulate species in Tibet
Spatial data structures
We can distinguish three types of spatial data structures
Areal data
Map of bird conservation regions (BCRs) showing the proportion of bird species within each region showing a declining trend
Geostatistical data
In geostatistical data, measurements of a continuous process are taken at a set of fixed locations.
Scotland river temperature monitoring network
Point pattern data
Occurrence records of four ungulate species in Tibet,
Spatial data structures
We can distinguish three types of spatial data structures
Areal data
Map of bird conservation regions (BCRs) showing the proportion of bird species within each region showing a declining trend
Geostatistical data
Scotland river temperature monitoring network
Point pattern data
In point pattern data we record the locations where events occur (e.g. trees in a forest, earthquakes) and the coordinates of such occurrences are our data.
Occurrence records of four ungulate species in Tibet,
Types of spatial data
Discrete space: - Data on a spatial grid (areal data)
Continuous space: - Geostatistical (geo-referenced) data - Spatial point pattern data
Model components are used to reflect spatial dependence structures in discrete and continuous space.
Areal Data
How do we model this?
Imagine we have animal counts in each region. We can model them as Poisson
\[
y_i \sim \mathrm{Poisson}(e^{\eta_i})
\] How do we model the linear predictor \(\eta_i\)?
We could model the number of animals in each region independently
Where \(Q\) denotes the precision matrix defined as
\[
Q_{i,j} = \begin{cases}
d_i, & i = j \\
-1, & i \sim j \\
0, &\text{otherwise}
\end{cases}
\]
This structure matrix directly defines the neighbourhood structure and is sparse.
Modelling spatial similarity
The ICAR model accounts only for spatially structured variability and does not include a limiting case where no spatial structure is present.
We typically add an unstructured random effect \(z_i|\tau_v \sim N(0,\tau_{z}^{-1})\)
The resulting model \(v_i = u_i + z_i\) is known as the Besag-York-Mollié model (BYM)
The structured spatial effect is controlled by \(\tau_u\) which control the degree of smoothing:
Higher \(\tau_u\) values lead to stronger smoothing (less spatial variability).
Lower \(\tau_u\) values allow for greater local variation.
In the next example we will illustrate how to fit this model using inlabru.
Example: Respiratory hospitalisation in Glasgow
In this example we model the number of respiratory hospitalisations across Intermediate Zones (IZ) that make up the Greater Glasgow and Clyde health board in Scotland.
In epidemiology, disease risk is assessed using Standardized Mortality Ratios (SMR):
\[ SMR_i = \dfrac{Y_i}{E_i} \]
A value \(SMR > 1\) indicates a high risk area.
A value \(SMR<1\) suggests a low risk area.
Example: Respiratory hospitalisation in Glasgow
Stage 1: We assume the responses are Poisson distributed: \[
\begin{aligned}y_i|\eta_i & \sim \text{Poisson}(E_i\lambda_i)\\\text{log}(\lambda_i) = \color{#FF6B6B}{\boxed{\eta_i}} & = \color{#FF6B6B}{\boxed{\beta_0 + \beta_1 \mathrm{pm10} + u_i + z_i} }\end{aligned}
\]
Stage 2:\(\eta_i\) is a linear function of four components: an intercept, pm10 effect , a spatially structured effect \(u\) and an unstructured iid random effect \(z\): \[
\eta_i = \beta_0 + u_i + z_i
\]
Stage 3:\(\{\tau_{z},\tau_u\}\): Precision parameters for the random effects
The latent field is \(\mathbf{x}= (\beta_0, \beta_1, u_1, u_2,\ldots, u_n,z_1,...)\), the hyperparameters are \(\boldsymbol{\theta} = (\tau_u,\tau_z)\), and must be given a prior.
# define model componentcmp =~Intercept(1) +pm10(pm10, model ="linear") +space(space, model ="besag", graph = Q) +iid(space, model ="iid") # define model predictoreta = observed ~ Intercept + space + iid# build the observation modellik =bru_obs(formula = formula, family ="poisson",E = expected,data = resp_cases)# fit the modelfit =bru(cmp, lik)
# define model componentcmp =~Intercept(1) +pm10(pm10, model ="linear") space(space, model ="besag", graph = Q) +iid(space, model ="iid") # define model predictoreta = observed ~ Intercept + space + iid# build the observation modellik =bru_obs(formula = formula, family ="poisson",E = expected,data = resp_cases)# fit the modelfit =bru(cmp, lik)
# define model componentcmp =~Intercept(1) +pm10(pm10, model ="linear")space(space, model ="besag", graph = Q) +iid(space, model ="iid") # define model predictoreta = observed ~ Intercept + space + iid# build the observation modellik =bru_obs(formula = formula, family ="poisson",E = expected,data = resp_cases)# fit the modelfit =bru(cmp, lik)
# define model componentcmp =~Intercept(1) +pm10(pm10, model ="linear")space(space, model ="besag", graph = Q) +iid(space, model ="iid") # define model predictoreta = observed ~ Intercept + space + iid# build the observation modellik =bru_obs(formula = formula, family ="poisson",E = expected,data = resp_cases)# fit the modelfit =bru(cmp, lik)
neighbourhood structure
library(spdep)W.nb <-poly2nb(GGHB.IZ,queen =TRUE)
Example: Respiratory hospitalisation in Glasgow
Posterior summaries
mean
0.025quant
0.975quant
Intercept
−1.52
−1.91
−1.12
pm10
0.09
0.07
0.12
Precision for space
2.82
2.28
3.45
Precision for iid
21,825.18
1,420.88
85,676.44
Estimated relative risks
Geostatistical modelling
Geostatistical Data
Measurements of an a spatial continuous variable of interest at a set of fixed locations.
This could be data from samples taken in locations across a region of interest (e.g., 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 and/or to model the relationship with covariates.
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?
Gaussian Random Fields
If we have a process that is occurring everywhere in space, it is natural to try to model it using some sort of function.
If \(z\) is a vector of observations of \(z(\mathbf{s})\) at different locations, we want this to be normally distributed:
\[
\mathbf{z} = (z(\mathbf{s}_1),\ldots,z(\mathbf{s}_m)) \sim \mathcal{N}(0,\Sigma),
\] where \(\Sigma_{ij} = \mathrm{Cov}(z(\mathbf{s}_i),z(\mathbf{s}_j))\) is a dense \(m \times m\) matrix.
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
Stationary random fields
A GRF is stationary if:
has mean zero.
the covariance between two points depends only on the distance and direction between those points.
It is isotropic if the covariance only depends on the distance between the points.
The SPDE approach
The goal: approximate the GRF using a triangulated mesh via the so-called SPDE approach.
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
A GF with Matérn covariance \(c_{\nu}(d;\sigma,\rho)\) is a solution to a particular PDE.
This solution is then approximated using a finite combination of piece-wise linear basis functions defined on a triangulation .
The solution is completely defined by a Gaussian vector of weights (defined on the triangulation vertices) with zero mean and a sparse precision matrix.
How do we choose sensible priors for \(\sigma,\rho\)?
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
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
SPDE models
We call spatial Markov models defined on a mesh SPDE models.
SPDE models have 3 parts
A mesh
A range parameter \(\kappa\)
A precision parameter \(\tau\)
Lets see how this is done in practice!
Example: Species Distribution Model
In the next example, we will explore data on the Pacific Cod (Gadus macrocephalus) from a trawl survey in Queen Charlotte Sound. The
The dataset contains presence/absence data in 2003 as well as depth covariate infromation.
Example: Species Distribution Model
Stage 1 Model for the response\[
y(s)|\eta(s)\sim\text{Binom}(1, p(s))
\]
Stage 2 Latent field model \[
\eta(s) = \text{logit}(p(s)) = \beta_0 + f( x(s)) + \omega(s)
\]
Stage 3 Hyperparameters
Bayesian Geostatistics
Stage 1 Model for the response \[
y(s)|\eta(s)\sim\text{Binom}(1, p(s))
\]
A Gaussian field \(\omega(s)\) (will discuss this later..)
Stage 3 Hyperparameters
Example: Species Distribution Model
Stage 1 Model for the response \[
y(s)|\eta(s)\sim\text{Binom}(1, p(s))
\]
Stage 2 Latent field model \[
\eta(s) = \text{logit}(p(s)) = \beta_0 + \beta_1 x(s) + \omega(s)
\]
Stage 3 Hyperparameters
Precision for the smooth function \(f(\cdot)\)
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
# define model componentcmp =~-1+Intercept(1) +depth_smooth(log(depth), model='rw2') +space(geometry, model = spde_model)# define model predictoreta = present ~ Intercept + depth_smooth + space# build the observation modellik =bru_obs(formula = eta,data = pcod_sf,family ="binomial")# fit the modelfit =bru(cmp, lik)
Model Results
Point Processes
Point process data
Many of the ecological and environmental processes of interest can be represented by a spatial point process or can be viewed as an aggregation of one.
Many contemporary data sources collect georeferenced information about the location where an event has occur (e.g., species occurrence, wildfire, flood events).
This point-based information provides valuable insights into ecosystem dynamics.
Defining a Point Process
Consider a fixed geographical region \(A\).
The set of locations at which events occur are denoted by \(\mathbf{s} = (\mathbf{s}_1, \ldots, \mathbf{s}_n)\).
We let \(N(A)\) be a random variable which represents the total number of events in every subset of region \(A\).
Our primary interest is in measuring where events occur, so the locations are our data.
Homogeneous Poisson Process
A simplestpoint process model is the homogeneous Poisson process (HPP).
The likelihood of a point pattern \(\mathbf{y} = \left[ \mathbf{s}_1, \ldots, \mathbf{s}_n \right]^\intercal\) distributed as a HPP with intensity \(\lambda\) and observation window \(\Omega\) is
\(|\Omega|\) is the size of the observation window.
\(\lambda\) is the expected number of points per unit area.
\(|\Omega|\lambda\) the total expected number of points in the observation window.
A key property of a Poisson process is that the number of points within any subset\(B\) of region \(A\) is Poisson distributed with constant rate \(|B|\lambda\).
Inhomogeneous Poisson process
The inhomogeneous Poisson process has a spatially varying intensity \(\lambda(\mathbf{s})\).
If the case of an HPP the integral in the likelihood can easily be computed as \(\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} =|\Omega|\lambda\)
For IPP, integral in the likelihood has to be approximated numerically as a weighted sum.
Inhomogeneous Poisson process
This integral is approximated as \(\sum_{j=1}^J w_j \lambda(\mathbf{s}_j)\)
\(w_j\) are the integration weights
\(\mathbf{s}_j\) are the quadrature locations.
This serves two purposes:
Approximating the integral
re-writing the inhomogeneous Poisson process likelihood as a regular Poisson likelihood.
Inhomogeneous Poisson process
The idea behind the trick to rewrite the approximate likelihood is to introduce a dummy vector \(\mathbf{z}\) and an integration weights vector \(\mathbf{w}\) of length \(J + n\)
This is similar to a product of Poisson distributions with means \(\eta_i\), exposures \(w_i\) and observations \(z_i\).
This is the basis for the implementation of Cox process models in inlabru, which can be specified using family = "cp".
Limitations with IPP
IPP models assume that data points are conditionally independent given the covariates, meaning that any spatial variation is fully explained by environmental and sampling factors.
Unmeasured endogenous and exogenous factors can create spatial dependence.
Ignoring them can lead to bias in our conclusions.
The Log-Gaussian Cox Process
Log-Gaussian Cox processes (LGCP) extend the IPP by allowing the intensity function to vary spatially according to a structured spatial random effect, i.e. the intensity is random
The events are then assumed to be independent given the covariates and \(\xi(s)\) - a GMRF with Matérn covariance.
inlabru has implemented some integration schemes that are especially well suited to integrating the intensity in models with an SPDE effect.
Example: Forest fires in Castilla-La Mancha
In this example we model the location of forest fires in the Castilla-La Mancha region of Spain between 1998 and 2007.
We are now going to use the elevation as a covariate to explain the variability of the intensity \(\lambda(s)\) over the domain of interest and a spatially structured SPDE model.
