In spatial analysis, the starting point of the SPDE approach is the linear fractional stochastic partial differential equation:
\[ (\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s)) = \mathcal{W}(s), \quad s \in \mathbb{R}^d \]
Where:
Traditionally, we define \(\xi(s)\) directly using the Matérn covariance function:
\[ Cov(\xi(s_i), \xi(s_j)) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}(\kappa ||h||)^\nu K_\nu(\kappa ||h||) \]
However, using this covariance matrix for \(n\) locations requires \(O(n^3)\) operations—impossible for large datasets (e.g., satellite data with millions of pixels).
The SPDE equation is a mathematical shortcut:
In plain English: The SPDE disguises a global, dense covariance matrix as a local, sparse “nearest-neighbor” model that computers can solve incredibly fast.
Yes, but with important distinctions:
| Scenario | Mathematical Formulation | What is \(\xi(s)\)? |
|---|---|---|
| No covariates, no error | \(Y(s) = \xi(s)\) | The raw observed data. |
| No covariates, with error | \(Y(s) = \xi(s) + \varepsilon(s)\) | The true, noise-free spatial surface. |
| With covariates, no error | \(Y(s) = \mathbf{X}(s)\boldsymbol{\beta} + \xi(s)\) | The spatial residual after removing covariate effects. |
| With covariates and error (Standard) | \(Y(s) = \mathbf{X}(s)\boldsymbol{\beta} + \xi(s) + \varepsilon(s)\) | The latent (hidden) spatial random effect. |
Where:
The SPDE equation itself NEVER changes when you add covariates:
\[ (\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s)) = \mathcal{W}(s) \]
Covariates are added to the observation equation (the likelihood), not inside the differential operator.
Data level (Likelihood):
\[ Y(s_i) \mid \xi(s_i) \sim \text{Normal}\left(\mathbf{X}(s_i)\boldsymbol{\beta} + \xi(s_i), \sigma^2_e\right) \]
Process level (SPDE):
\[ (\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s)) = \mathcal{W}(s) \]
Parameter level (Priors):
\[ \boldsymbol{\beta} \sim \text{Normal}(0, 1000) \]
\[ \log(\kappa) \sim \text{Normal}(0, 1) \]
\[ \log(\tau) \sim \text{Normal}(0, 1) \]
\[ \sigma^2_e \sim \text{Inverse-Gamma}(1, 0.01) \]
In the simulation example from Krainski and Lindgren (2013), the model is:
\[ y_i \sim \text{Normal}(\eta_i, \sigma^2_e), \quad i = 1, \dots, 200 \]
where:
\(\sigma^2 = 5\) is the variance of the latent spatial field \(\xi(s)\):
\[ \text{Var}(\xi(s)) = \sigma^2 = 5 \]
This means:
In this simulation, the total variance of your observations \(y_i\) comes from two independent sources:
Since \(\xi(s)\) and \(\varepsilon_i\) are independent:
\[ \text{Var}(y_i) = \text{Var}(\xi_i + \varepsilon_i) = \sigma^2 + \sigma^2_e = 5 + 0.3 = 5.3 \]
So the total variance of your simulated observations is 5.3.
This is often called the proportion of spatial variance or the “fraction”:
\[ \frac{\sigma^2}{\sigma^2 + \sigma^2_e} = \frac{5}{5 + 0.3} = \frac{5}{5.3} \approx 0.943 \]
Interpretation: About 94.3% of the total variation in your data is due to spatial correlation, and only 5.7% is pure measurement noise. This means the spatial field is very strong and the data are highly spatially correlated.
Example: - If at location \(s\), \(\xi(s) = 2\) (slightly above average spatial effect) and \(\varepsilon = 0.2\) (small measurement error), then: - \(y = 10 + 2 + 0.2 = 12.2\)
In the SPDE formulation:
\[ (\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s)) = \mathcal{W}(s) \]
The variance \(\sigma^2\) is related to the SPDE parameters \(\tau\), \(\kappa\), and \(\alpha\):
For a Matérn field in \(\mathbb{R}^d\):
\[ \sigma^2 = \frac{\Gamma(\nu)}{\Gamma(\nu + d/2)} \cdot \frac{1}{(4\pi)^{d/2} \kappa^{2\nu} \tau^2} \]
where:
So in the simulation, the parameter values were chosen such that:
| Parameter | Value | Meaning |
|---|---|---|
| \(b_0\) | 10 | Global mean |
| \(\sigma^2\) | 5 | Variance of spatial field \(\xi(s)\) |
| \(\sigma^2_e\) | 0.3 | Variance of measurement error |
| \(\kappa\) | 7 | Controls spatial range (larger = shorter range) |
| \(\lambda\) | 1 | Smoothness parameter (\(\nu = 1\)) |
“The observed rainfall \(Y\) at location \(s\) equals the linear effect of altitude and other covariates, plus a spatially correlated residual \(\xi(s)\), plus independent measurement noise. The SPDE equation forces \(\xi(s)\) to follow a Matérn covariance function with range \(\kappa\) and smoothness \(\alpha\).”
Key takeaway: The SPDE only models the residual spatial dependence after accounting for covariates. This separates:
| Concept | Equation | Role |
|---|---|---|
| SPDE | \((\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s)) = \mathcal{W}(s)\) | Defines the spatial covariance of \(\xi(s)\). |
| Observation Model | \(Y(s) = \mathbf{X}(s)\boldsymbol{\beta} + \xi(s) + \varepsilon(s)\) | Links data \(Y\) to the spatial field and covariates. |
| Covariates | \(\mathbf{X}(s)\boldsymbol{\beta}\) | Explain large-scale trends (fixed effects). |
| Spatial Random Effect | \(\xi(s)\) | Captures residual spatial correlation (random effects). |
| Measurement Error | \(\varepsilon(s) \sim \text{N}(0, \sigma^2_e)\) | Captures independent noise. |
| Spatial Variance | \(\sigma^2 = 5\) | Magnitude of spatial variation in \(\xi(s)\). |
When we fit this model to real data:
Here is a minimal example of fitting this model using the
R-INLA package:
# Load required libraries
library(INLA)
library(sf)
# Assume we have:
# - coords: matrix of spatial coordinates
# - y: vector of observations
# - x1, x2: covariates
# 1. Create the mesh (triangulation of the spatial domain)
mesh <- inla.mesh.2d(
loc = coords,
max.edge = c(0.5, 1),
cutoff = 0.2
)
# 2. Build the SPDE model object
spde <- inla.spde2.matern(
mesh = mesh,
alpha = 2 # for 2D space, alpha = 2 gives smoothness nu = 1
)
# 3. Create the projection matrix (A) to map mesh to data locations
A <- inla.spde.make.A(mesh, loc = coords)
# 4. Create the stack (combine data, covariates, and SPDE)
stack <- inla.stack(
data = list(y = y),
A = list(A, 1),
effects = list(
list(spde = 1:spde$n.spde), # SPDE random effect
list(b0 = 1, x1 = x1, x2 = x2) # Fixed effects
)
)
# 5. Fit the model
formula <- y ~ 0 + b0 + x1 + x2 + f(spde, model = spde)
result <- inla(formula,
data = inla.stack.data(stack),
control.predictor = list(A = inla.stack.A(stack), compute = TRUE)
)
# 6. View results
summary(result)
# 7. Extract key parameters
# Range parameter
range <- inla.spde2.result(result, "spde", spde, do.transf = TRUE)
print(range$summary.range)
# Variance parameter
variance <- inla.spde2.result(result, "spde", spde, do.transf = TRUE)
print(variance$summary.variance)
# 8. Predict on a grid for mapping
# Create prediction coordinates
grid_points <- expand.grid(
x = seq(min(coords[,1]), max(coords[,1]), length = 50),
y = seq(min(coords[,2]), max(coords[,2]), length = 50)
)
# Project predictions to grid
A_pred <- inla.spde.make.A(mesh, loc = as.matrix(grid_points))
# Create prediction stack
stack_pred <- inla.stack(
data = list(y = NA),
A = list(A_pred, 1),
effects = list(
list(spde = 1:spde$n.spde),
list(b0 = 1, x1 = grid_points$x, x2 = grid_points$y)
)
)
# Combine with estimation stack
stack_full <- inla.stack(stack, stack_pred)
# Refit with predictions
result_pred <- inla(formula,
data = inla.stack.data(stack_full),
control.predictor = list(
A = inla.stack.A(stack_full),
compute = TRUE
)
)
# Extract predictions
pred_mean <- result_pred$summary.linear.predictor[
inla.stack.index(stack_full, "pred")$data, "mean"
]Here’s a more general version with multiple covariates:
# Suppose we have covariates in a data frame
df <- data.frame(
y = y,
coords_x = coords[,1],
coords_y = coords[,2],
altitude = altitude,
temperature = temperature,
precipitation = precipitation
)
# Create mesh
mesh <- inla.mesh.2d(
loc = as.matrix(df[, c("coords_x", "coords_y")]),
max.edge = c(0.5, 1),
cutoff = 0.2
)
# Build SPDE model
spde <- inla.spde2.matern(mesh, alpha = 2)
# Projection matrix
A <- inla.spde.make.A(mesh, loc = as.matrix(df[, c("coords_x", "coords_y")]))
# Stack with all covariates
stack <- inla.stack(
data = list(y = df$y),
A = list(A, 1),
effects = list(
list(spde = 1:spde$n.spde),
list(
intercept = 1,
altitude = df$altitude,
temperature = df$temperature,
precipitation = df$precipitation
)
)
)
# Formula with all covariates
formula <- y ~ 0 + intercept + altitude + temperature + precipitation +
f(spde, model = spde)
# Fit model
result <- inla(formula,
data = inla.stack.data(stack),
control.predictor = list(A = inla.stack.A(stack), compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE)
)
# Summary
summary(result)
# Extract spatial variance
variance <- inla.spde2.result(result, "spde", spde, do.transf = TRUE)
print(variance$summary.variance)
# Extract measurement error
print(result$summary.hyperpar)Here’s how we would simulate data with the parameters from the Krainski and Lindgren (2013) example:
# Set seed for reproducibility
set.seed(123)
# Parameters from the simulation
b0 <- 10 # Intercept
sigma2 <- 5 # Spatial variance
sigma2_e <- 0.3 # Measurement error variance
kappa <- 7 # Scale parameter
nu <- 1 # Smoothness parameter (lambda = 1)
# Generate 200 random locations in a unit square
n <- 200
coords <- cbind(runif(n), runif(n))
# Build the Matérn covariance matrix
# This uses the fields package for demonstration
library(fields)
# Compute distances
dist_matrix <- rdist(coords)
# Matérn covariance function
matern_cov <- function(h, kappa, nu, sigma2) {
# h: distance
# kappa: scale parameter
# nu: smoothness
# sigma2: variance
h <- h * kappa
if (nu == 0.5) {
return(sigma2 * exp(-h))
} else {
return(sigma2 * (h^nu) * besselK(h, nu) / (2^(nu-1) * gamma(nu)))
}
}
# Build covariance matrix
Sigma <- matrix(0, n, n)
for (i in 1:n) {
for (j in 1:n) {
Sigma[i,j] <- matern_cov(dist_matrix[i,j], kappa, nu, sigma2)
}
}
# Simulate spatial field
library(MASS)
xi <- mvrnorm(1, mu = rep(0, n), Sigma = Sigma)
# Simulate measurement error
epsilon <- rnorm(n, mean = 0, sd = sqrt(sigma2_e))
# Generate observations
y <- b0 + xi + epsilon
# Check variance decomposition
cat("Variance of xi:", var(xi), "\n")
cat("Variance of epsilon:", var(epsilon), "\n")
cat("Total variance of y:", var(y), "\n")
cat("Proportion spatial:", var(xi) / var(y), "\n")
# Fit SPDE model to the simulated data
mesh <- inla.mesh.2d(loc = coords, max.edge = c(0.2, 0.5))
spde <- inla.spde2.matern(mesh, alpha = 2)
A <- inla.spde.make.A(mesh, loc = coords)
stack <- inla.stack(
data = list(y = y),
A = list(A, 1),
effects = list(
list(spde = 1:spde$n.spde),
list(intercept = rep(1, n))
)
)
formula <- y ~ 0 + intercept + f(spde, model = spde)
result <- inla(formula,
data = inla.stack.data(stack),
control.predictor = list(A = inla.stack.A(stack), compute = TRUE)
)
# Compare estimated parameters
summary(result)The SPDE approach is revolutionary because:
The equation \((\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s)) = \mathcal{W}(s)\) is the engine; covariates are the steering wheel. They work together to model complex spatial phenomena efficiently.
Key understanding: When you see \(\sigma^2 = 5\) in the SPDE context, it’s the variance of the hidden spatial field \(\xi(s)\), not the observation noise! This represents the magnitude of spatial variation in your process after accounting for fixed effects.