Gaussian random fields on metric graphs

Introduction

In this vignette we will introduce how to work with Gaussian random fields on metric graphs. The main models are the Whittle–Matérn fields introduced in Bolin, Simas, and Wallin (2024) and Bolin, Simas, and Wallin (2023).

The package also has support for isotropic Gaussian processes, and in particular Gaussian processes with isotropic exponential covariance functions as introduced by Anderes, Møller, and Rasmussen (2020). Finally, Gaussian models based on the graph Laplacian, as introduced by Borovitskiy et al. (2021) are also supported, even though these do not defined Gaussian processes on the metric graph, but only at the vertices.

As an example throughout the vignette, we consider the following metric graph:

edge1 <- rbind(c(0,0),c(1,0))
edge2 <- rbind(c(0,0),c(0,1))
edge3 <- rbind(c(0,1),c(-1,1))
theta <- seq(from=pi,to=3*pi/2,length.out = 20)
edge4 <- cbind(sin(theta),1+ cos(theta))
edges = list(edge1, edge2, edge3, edge4)
graph <- metric_graph$new(edges = edges)
graph$plot()

For further details on the construction of metric graphs, see Working with metric graphs

Whittle–Matérn fields

The Whittle–Matérn fields are specified as solutions to the stochastic differential equation $(\kappa^2 - \Delta)^{\alpha/2} \tau u = \mathcal{W}$ on the metric graph $\Gamma$ . We can work with these models without and approximations if the smoothness parameter $\alpha$ is an integer, and this is what we focus on in this vignette. For details on the case of a general smoothness parameter, see Whittle–Matérn fields with general smoothness.

Sampling

As an example, let us simulate the field $u$ on the graph using $\alpha = 1$ . To do so, we first need to specify where to sample it. As a first example, let us specify some locations manually:

PtE <- cbind(rep(1:4, each = 4),
             rep(c(0.2, 0.4, 0.6, 0.8), times = 4))

sigma <- 1.3
alpha <- 1
range <- 0.2
u <- sample_spde(kappa = kappa, sigma = sigma, 
                 range = range,
                 graph = graph, PtE = PtE)
graph$plot(X = u, X_loc = PtE)

In many cases, one wants to sample the field at evenly spaced locations over the graph. To avoid having to specify such locations manually, we can first create a mesh on the graph

graph$build_mesh(h = 0.1)
graph$plot(mesh=TRUE)

In the command build_mesh, the argument h decides the largest spacing between nodes in the mesh. We can now sample the field on this mesh and plot the result as a function as follows:

u <- sample_spde(range = range, sigma = sigma, alpha = alpha,
                 graph = graph, type = "mesh")
graph$plot_function(X = u)

Let us construct a finer mesh, simulate the field, and visualize the simulation in 3D by specifying the plotly argument in the plot function:

graph$build_mesh(h = 0.01)


u <- sample_spde(range = range, sigma = sigma, alpha = alpha,
                 graph = graph, type = "mesh")
graph$plot_function(X = u, type = "plotly")

## Loading required namespace: plotly

Since $\alpha=1$ , these sample paths are continuous but not differentiable. To visualize the correlation structure of the field, we can compute and plot the covariances between some point and all other points in the graph as follows:

C <- spde_covariance(c(2, 0.2), range = range, sigma = sigma, alpha = 1,
                            graph = graph)
graph$plot_function(X = C, type = "plotly")

To obtain a field with differentiable sample paths, we can change to $\alpha=2$ . The corresponding covariance function then looks as follows:

C <- spde_covariance(c(2, 0.2), range = range, sigma = sigma, alpha = 2,
                            graph = graph)
graph$plot_function(X = C, type = "plotly")

Let us simulate a process with $\alpha=2$ as well:

u <- sample_spde(range = range, sigma = sigma, alpha = 2,
                 graph = graph, type = "mesh")
graph$plot_function(X = u, type = "plotly")

Inference

The MetricGraph package contains implementations of a linear mixed effects models, where the random effect can be Whittle–Matérn fields, and the model is observed under Gaussian measurement noise. We also implemented random effects being SPDE fields obtained from the graph Laplacian as well as models with isotropic covariances. In this section we will illustrate these methods. For the use of the Whittle–Matérn fields in more complicated hierarchical models, we recommend using the interfaces to the INLA and inlabru packages. See INLA interface of Whittle–Matérn fields and inlabru interface of Whittle–Matérn fields for further details on these.

Suppose that we want to estimate the model parameters of a Whittle–Matérn field $u(s)$ observed under Gaussian measurement noise. That is, we assume that we are given observations $y_i = u(s_i) + \varepsilon_i, \quad i=1,\ldots,n$ where $s_i\in \Gamma$ are observation locations and $\varepsilon_i$ are independent centered Gaussian variables $N(0,\sigma_e^2)$ representing measurement noise.

Let us start by generating some data like this and adding it to the metric graph. For further details on data manipulation on metric graphs, see Data manipulation on metric graphs.

range <- 0.2
sigma <- 1.3
sigma_e <- 0.1
alpha <- 1

n.obs.per.edge <- 75
PtE <- NULL
for(i in 1:graph$nE){
  #add locations sampled at random to each edge
  PtE <- rbind(PtE, cbind(rep(i, n.obs.per.edge), runif(n.obs.per.edge)))
}

u <- sample_spde(range = range, sigma = sigma, alpha = alpha,
                 graph = graph, PtE = PtE)

y <- u + sigma_e*rnorm(n.obs.per.edge * graph$nE)

df_data <- data.frame(y = y, edge_number = PtE[,1],
                        distance_on_edge = PtE[,2])

graph$clear_observations() # Removing previous observations
graph$add_observations(data = df_data, normalized = TRUE)

## Adding observations...

## list()

graph$plot(data = "y")

We can now use the graph_lme() function to fit the above model. By default a linear regression model is chosen. Since we want to fit a model with a latent model given by a Whittle-Matérn field with $\alpha=1$ , we should set the model argument to either 'alpha1' or to list(model = 'WhittleMatern', alpha = 1) (the first is more convenient but less decriptive). We will choose parameterization_latent as "spde" to obtain the estimated values for kappa and sigma. By default it provides estimated values for the range parameter and sigma.

res <- graph_lme(y ~ -1, graph = graph, model = 'WM1')

summary(res)

## 
## Latent model - Whittle-Matern with alpha = 1
## 
## Call:
## graph_lme(formula = y ~ -1, graph = graph, model = "WM1")

## 
## No fixed effects.

## 
## Random effects:
##       Estimate Std.error z-value
## tau   0.170788  0.009253  18.458
## kappa 6.469880  1.830364   3.535
## 
## Random effects (Matern parameterization):
##       Estimate Std.error z-value
## sigma   1.6277    0.2146   7.586
## range   0.3091    0.0860   3.594
## 
## Measurement error:
##          Estimate Std.error z-value
## std. dev  0.08635   0.02791   3.094
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Log-Likelihood:  -241.4271 
## Number of function calls by 'optim' = 16
## Optimization method used in 'optim' = L-BFGS-B
## 
## Time used to:     fit the model =  0.84931 secs

We can also take a glance at res:

glance(res)

## # A tibble: 1 × 9
##    nobs  sigma logLik   AIC   BIC deviance df.residual model         alpha
##   <int>  <dbl>  <dbl> <dbl> <dbl>    <dbl>       <dbl> <chr>         <dbl>
## 1   300 0.0864  -241.  489.  500.     483.         297 WhittleMatern     1

Let us now compare with the true values:

sigma_e_est <- res$coeff$measurement_error
sigma_est <- res$matern_coeff$random_effects[1]
range_est <- res$matern_coeff$random_effects[2]
results <- data.frame(sigma_e = c(sigma_e, sigma_e_est),
                      sigma = c(sigma, sigma_est),
                      range = c(range, range_est),
                      row.names = c("Truth", "Estimate"))
print(results)

##             sigma_e    sigma     range
## Truth    0.10000000 1.300000 0.2000000
## Estimate 0.08635469 1.627717 0.3091247

Given these estimated parameters, we can now do kriging to estimate the field at locations in the graph. As an example, we now estimate the field on the regular mesh that we previously constructed.

u_est <- predict(res, data.frame(edge_number = graph$mesh$VtE[,1],
                        distance_on_edge = graph$mesh$VtE[,2]), normalized = TRUE)

graph$plot_function(X = u_est$mean, type = "plotly")

David Bolin, Alexandre B. Simas, and Jonas Wallin

Created: 2022-11-23. Last modified: 2024-12-01.

Introduction

Whittle–Matérn fields

Sampling

Inference