Ngme2 Matern SPDE model • ngme2

In this vignette, we will introduce the SPDE Matérn model in ngme2. First we introduce a little about Gaussian process.

Gaussian Matérn Fields

Gaussian processes and random fields provide powerful methods for modeling spatial and spatio-temporal dependence structures. In spatial statistics, particularly in geostatistics, Gaussian fields (GFs) play a fundamental role due to their mathematical tractability and flexibility.

A standard geostatistical model can be expressed as: $Y_i = x(\mathbf{s}_i) + \varepsilon_i, \quad i=1,\ldots,N, \quad \varepsilon_i\sim N(0, \sigma_\varepsilon^2),$ $x(\mathbf{s}) \sim GP\left(\sum_{k=1}^{n_b} b_k(\mathbf{s})w_k, c(\mathbf{s},\mathbf{s}')\right),$ where:

$Y_i$ represents observations at spatial locations $\mathbf{s}_i$
$N$ is the number of spatial observations
$GP(m,c)$ denotes a Gaussian process with mean function $m$ and covariance function $c$
$n_b$ is the number of basis functions in the mean structure
$\{b_k(\cdot)\}_{k=1}^{n_b}$ are basis functions (often polynomials or splines)
$w_k$ are coefficients to be estimated
$\varepsilon_i$ represents measurement error with variance $\sigma_\varepsilon^2$

Among the various covariance functions available for modeling spatial dependence, the Matérn covariance function is particularly valued for its flexibility:

$c(\mathbf{s}, \mathbf{s}') = \frac{\sigma^2}{\Gamma(\nu)2^{\nu-1}}(\kappa \|\mathbf{s}-\mathbf{s}'\|)^\nu K_\nu(\kappa\|\mathbf{s}-\mathbf{s}'\|),$

This function is parameterized by:

$\sigma^2$ : the marginal variance, controlling the overall variability of the process
$\kappa > 0$ : the spatial scale parameter, inversely related to the range of correlation.
$\nu > 0$ : the smoothness parameter, determining the differentiability of the field
$\Gamma(\cdot)$ : the Gamma function
$K_\nu(\cdot)$ : the modified Bessel function of the second kind of order $\nu$

The Matérn family is particularly useful because:

It includes the exponential covariance ( $\nu = 0.5$ ) and the Gaussian covariance (as $\nu \to \infty$ ) as special cases
The parameter $\nu$ directly controls the mean-square differentiability of the process
The effective range (distance at which correlation drops to approximately 0.1) is roughly $\sqrt{8\nu}/\kappa$

Traditional approaches to fitting these models rely on maximum likelihood estimation, which becomes computationally prohibitive for large datasets. The computational complexity typically scales as $\mathcal{O}(N^3)$ , making it impractical for modern spatial datasets with thousands or millions of observations.

The SPDE approach with Gaussian noise

It is well-known (Whittle, 1963) that a Gaussian process $u(\mathbf{s})$ with Matérn covariance function solves the stochastic partial differential equation (SPDE)

$$\begin{equation}\label{spde} (\kappa^2 -\Delta)^{\frac{\alpha}{2}} u = \mathcal{W}\quad \hbox{in } \mathcal{D}, \end{equation}$$ where $\Delta = \sum_{i=1}^d \frac{\partial^2}{\partial_{x_i^2}}$ is the Laplacian operator, $\mathcal{W}$ is the Gaussian spatial white noise on $\mathcal{D}=\mathbb{R}^d$ , and $\alpha = \nu + d/2$ .

Inspired by this relation between Gaussian processes with Matérn covariance functions and solutions of the above SPDE, Lindgren et al. (2011) constructed computationally efficient Gaussian Markov random field approximations of $u(\mathbf{s})$ , where the domain $\mathcal{D}\subsetneq \mathbb{R}^d$ is bounded and $\alpha\in\mathbb{N}$ . The approximate solutions of the SPDE are obtained through a finite element discretization.

In ngme2, we use the matern function to specify the SPDE Matérn model. For example, we can use the following code to specify the SPDE Matérn model in 1d case with Gaussian noise:

f(map = 1:10, 
  model = "matern",
  mesh = fmesher::fm_mesh_1d(1:10),
  kappa = 1, # the spatial scale parameter
  noise = noise_normal(sigma = 1)
)
#> Model type: Matern
#>     kappa = 1
#> Noise type: NORMAL
#> Noise parameters: 
#>     sigma = 1

The SPDE approach with non-Gaussian noise

Then we will describe how to generalize this approach with non-Gaussian noise. Our goal now is to describe the SPDE approach when the noise is non-Gaussian. The motivation for handling non-Gaussian noise comes from the fact that many features cannot not be handled by Gaussian noise. Some of these reasons are:

Skewness;
Heavier tails;
Jumps in the sample paths;
Asymmetries in the sample paths.

Non-Gaussian Matérn fields

The idea is to replace the Gaussian white noise $\mathcal{W}$ in the SPDE by a non-Gaussian white noise $\dot{\mathcal{M}}$ : $(\kappa^2 - \Delta)^\beta u = \dot{\mathcal{M}}.$ The solution $u$ will have Matérn covariance function, but their marginal distributions will be non-Gaussian.

We will consider the same setup. More precisely, we consider $V_n = {\rm span}\{\varphi_1,\ldots,\varphi_n\}$, where $\varphi_i(\mathbf{s}), i=1,\ldots, n$ are piecewise linear basis functions obtained from a triangulation of $\mathcal{D}$ and we approximate the solution $u$ by $u_n$ , where $u_n$ is written in terms of the basis functions as $u_n(\mathbf{s}) = \sum_{i=1}^n w_i \varphi_i(\mathbf{s}).$ In the right-hand side we obtain a random vector $\mathbf{f} = (\dot{\mathcal{M}}(\varphi_1),\ldots, \dot{\mathcal{M}}(\varphi_n)),$ where the functional $\dot{\mathcal{M}}$ is given by $\dot{\mathcal{M}}(\varphi_j) = \int_{\mathcal{D}} \varphi_j(\mathbf{s}) d\mathcal{M}(\mathbf{s}).$ By considering $\mathcal{M}$ to be a type-G Lévy process, we obtain that $\mathbf{f}$ has a joint distribution that is easy to handle.

We say that a Lévy process is of type G if its increments can be represented as location-scale mixtures: $\gamma + \mu V + \sigma \sqrt{V}Z,$ where $\gamma, \mu$ are parameters, $Z\sim N(0,1)$ and is independent of $V$ , and $V$ is a positive infinitely divisible random variable.

Therefore, given a vector $\mathbf{V} = (V_1,\ldots,V_n)$ of independent stochastic variances (in our case, positive infinitely divisible random variables), we obtain that $$\mathbf{f}|\mathbf{V} \sim N(\gamma + \mu\mathbf{V}, \sigma^2{\rm diag}(\mathbf{V})).$$ So, if we consider, for instance, the non-fractional and non-Gaussian SPDE $(\kappa^2 - \Delta) u = \dot{\mathcal{M}},$ we obtain that the FEM weights $\mathbf{w} = (w_1,\ldots,w_n)$ satisfy $$\mathbf{w}|\mathbf{V} \sim N(\mathbf{K}^{-1}(\gamma+\mu\mathbf{V}), \sigma^2\mathbf{K}^{-1}{\rm diag}(\mathbf{V})\mathbf{K}^{-1}),$$ where $\mathbf{K} = \kappa^2\mathbf{C}+\mathbf{G}$ is the discretization of the differential operator.

In ngme2, we use the matern function to specify the SPDE Matérn model. For example, we can use the following code to specify the SPDE Matérn model in 1d case with NIG noise:

f(map = 1:10, 
  model = "matern",
  mesh = fmesher::fm_mesh_1d(1:10),
  kappa = 1, # the spatial scale parameter
  noise = noise_nig(mu = 0, sigma = 1, nu = 0.5)
)
#> Model type: Matern
#>     kappa = 1
#> Noise type: NIG
#> Noise parameters: 
#>     mu = 0
#>     sigma = 1
#>     nu = 0.5

Using SPDE Matérn model in ngme2

Simulation

We first use the simulation function to generate the data.

# define the domain
library(ggplot2)
library(viridis)
#> Loading required package: viridisLite
pl01 <- cbind(c(0, 1, 1, 0, 0) * 10, c(0, 0, 1, 1, 0) * 5)
mesh <- fmesher::fm_mesh_2d(
  loc.domain = pl01, cutoff = 0.1,
  max.edge = c(0.3, 10)
)

n_obs <- 1000
loc <- cbind(runif(n_obs, 0, 10), runif(n_obs, 0, 5))

# define the simulated model
true_noise <- noise_nig(mu=-2, sigma=1, nu=0.5)
true_model <- f(
  map = loc,
  model="matern",
  theta_K = log(4),
  mesh = mesh,
  noise = true_noise
)

# simulate the data
W <- simulate(true_model)[[1]]
Y <- W + rnorm(n_obs, sd=0.5)

plot(mesh)
points(loc, col="red", pch=16)


##### Fit the model
out <- ngme(
  Y ~ 0 + f(loc,
    model="matern",
    name="spde",
    mesh = mesh,
    noise=noise_nig(),
  ),
  data = data.frame(Y = Y),
  control_opt = control_opt(
    iterations = 3000,
    optimizer = precond_sgd(),
    solver = "supernodal",
    rao_blackwellization = TRUE,
    n_parallel_chain = 4,
    std_lim = 0.01
  )
)
#> Starting estimation... 
#> 
#> Starting posterior sampling... 
#> Posterior sampling done! 
#> Note:
#>       1. Use ngme_post_samples(..) to access the posterior samples.
#>       2. Use ngme_result(..) to access different latent models.
out
#> *** Ngme object ***
#> 
#> Fixed effects: 
#>   None
#> 
#> Models: 
#> $spde
#>   Model type: Matern
#>       kappa = 3.94
#>   Noise type: NIG
#>   Noise parameters: 
#>       mu = -1.7
#>       sigma = 1.15
#>       nu = 0.464
#> 
#> Measurement noise: 
#>   Noise type: NORMAL
#>   Noise parameters: 
#>       sigma = 0.492
#> 
#> 
#> Number of replicates is  1
traceplot(out, "spde", hline=c(4, -2, 1, 0.5)) # compare with true parameters

#> Last estimates:
#> $kappa
#> [1] 3.945232
#> 
#> $mu
#> [1] -1.704709
#> 
#> $sigma
#> [1] 1.15019
#> 
#> $nu
#> [1] 0.4643931
traceplot(out, hline=0.5)

#> Last estimates:
#> $sigma
#> [1] 0.493301

# Compare with the true density 
plot(true_noise,
  out$replicates[[1]]$models[[1]]$noise)

Real data example (argo float data)

We use the argo_float data to illustrate how to use the SPDE Matérn model in ngme2.

# 2d example
data(argo_float)
head(argo_float)
#>       lat     lon           sal       temp
#> 1 -64.078 175.821 -0.0699508100  0.4100305
#> 2 -63.760 162.917 -0.0320931260 -0.2588680
#> 3 -63.732 163.294 -0.0008063143 -0.1151362
#> 4 -63.700 162.568 -0.0209534220 -0.2378965
#> 5 -63.269 169.623  0.0409914840  0.3375048
#> 6 -63.113 171.526  0.0269408910  0.2145556
# take longitude and latitude to build the mesh

max.edge    <- 1
bound.outer <- 5
loc_2d <- unique(cbind(argo_float$lon, argo_float$lat))
# nrow(loc) == nrow(dat) no replicate
argo_mesh <- fmesher::fm_mesh_2d(loc = loc_2d,
                    # the inner edge and outer edge
                    max.edge = c(1,5),
                    cutoff = 0.3,
                    # offset extension distance inner and outer extenstion
                    offset = c(max.edge, bound.outer)
)

Let’s use the previous argo_float spatial (2d) example. First we explore the how the data look like:

# tempearture
ggplot(data=argo_float) +
  geom_point(aes(
    x = loc_2d[, 1], y = loc_2d[, 2],
    colour = temp
  ), size = 2, alpha = 1) +
  scale_color_gradientn(colours = viridis(100))


# salinity
ggplot(data=argo_float) +
  geom_point(aes(
    x = loc_2d[, 1], y = loc_2d[, 2],
    colour = sal
  ), size = 2, alpha = 1) +
  scale_color_gradientn(colours = viridis(100))

Next, we specify a model formula, and then fit the model.

formula <- temp ~ sal + f(
  loc_2d, 
  model = "matern", 
  name = "spde",
  mesh=argo_mesh, 
  noise = noise_nig()
)

out <- ngme(
  formula = formula,
  family = "nig",
  data = argo_float,
  control_opt = control_opt(
    n_parallel_chain = 4,
    iterations = 2000,
    seed = 7,
    optimizer = precond_sgd(),
    solver = "supernodal",
    print_check_info = FALSE
  )
)
#> Starting estimation... 
#> 
#> Starting posterior sampling... 
#> Posterior sampling done! 
#> Note:
#>       1. Use ngme_post_samples(..) to access the posterior samples.
#>       2. Use ngme_result(..) to access different latent models.
out
#> *** Ngme object ***
#> 
#> Fixed effects: 
#> (Intercept)         sal 
#>       0.868       8.030 
#> 
#> Models: 
#> $spde
#>   Model type: Matern
#>       kappa = 0.733
#>   Noise type: NIG
#>   Noise parameters: 
#>       mu = -0.851
#>       sigma = 2.09
#>       nu = 0.475
#> 
#> Measurement noise: 
#>   Noise type: NIG
#>   Noise parameters: 
#>       mu = -0.0604
#>       sigma = 0.655
#>       nu = 0.0788
#> 
#> 
#> Number of replicates is  1

traceplot(out, "spde")

#> Last estimates:
#> $kappa
#> [1] 0.9223424
#> 
#> $mu
#> [1] -0.8479764
#> 
#> $sigma
#> [1] 2.136856
#> 
#> $nu
#> [1] 0.4900381
traceplot(out)

#> Last estimates:
#> $mu
#> [1] 0.1331467
#> 
#> $sigma
#> [1] 0.6392965
#> 
#> $nu
#> [1] 0.0790425
#> 
#> $`fixed effect 1`
#> [1] 0.8914478
#> 
#> $`fixed effect 2`
#> [1] 8.000974

Non-stationary SPDE Matérn model

In stationary case, the $K = \kappa^2 C + G$ . However, if we allow $\kappa$ to change over space, we can introduce non-stationary version, where $K = diag(\kappa) C diag(\kappa) + G$ , where $\kappa = \exp(B_K \theta_K)$ , $B_K$ is the basis matrix user provide.