# Gaussian random fields on metric graphs

#### David Bolin, Alexandre B. Simas, and Jonas Wallin

#### Created: 2022-11-23. Last modified: 2024-07-19.

Source:`vignettes/random_fields.Rmd`

`random_fields.Rmd`

## 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, plotly = TRUE)
```

`## 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, plotly = TRUE)
```

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, plotly = TRUE)
```

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, plotly = TRUE)
```

### 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...`

`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`

.

```
##
## 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.170764 0.009255 18.452
## kappa 6.535432 1.829112 3.573
##
## Random effects (Matern parameterization):
## Estimate Std.error z-value
## sigma 1.61976 0.21118 7.670
## range 0.30602 0.08422 3.633
##
## Measurement error:
## Estimate Std.error z-value
## std. dev 0.08637 0.02792 3.093
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -241.609
## Number of function calls by 'optim' = 15
## Optimization method used in 'optim' = L-BFGS-B
##
## Time used to: fit the model = 1.07966 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 -242. 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.08636943 1.619762 0.3060242
```

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, plotly = TRUE)
```

We can, alternatively, use the `augment()`

function:

```
pred_aug <- augment(res, data.frame(edge_number = graph$mesh$VtE[,1],
distance_on_edge = graph$mesh$VtE[,2]), normalized = TRUE)
p <- graph$plot_function(X = pred_aug[[".fitted"]], plotly = TRUE)
graph$plot(data = "y", p=p, plotly=TRUE)
```