MetricGraph: Random Fields on Metric Graphs

David Bolin, Alexandre B. Simas, and Jonas Wallin

Created: 2022-11-26. Last modified: 2024-12-13.

Introduction

There has lately been much interest in statistical modeling of data on compact metric graphs such as street or river networks based on Gaussian random fields.

The R package MetricGraph contains functions for working with data and random fields on compact metric graphs. The main functionality is contained in the metric_graph class, which is used for specifying metric graphs, adding data to them, visualization, and other basic functions that are needed for working with data and random fields on metric graphs. The package also implements three types of Gaussian processes on metric graphs: The Whittle–Matérn fields introduced by Bolin, Simas, and Wallin (2024) and Bolin, Simas, and Wallin (2023), Gaussian processes with isotropic covariance functions Anderes, Møller, and Rasmussen (2020), and Gaussian models based on the graph Laplacian Borovitskiy et al. (2021).

Basic statistical tasks such as likelihood evaluation and prediction is implemented for these three types of models in MetricGraph. Further, the package also contains interfaces to R-INLA Lindgren and Rue (2015), package available from http://R-INLA. org/download/ and inlabru Bachl et al. (2019) that facilitates using those packages for full Bayesian inference of general Latent Gaussian Models (LGMs) that includes Whittle–Matérn fields on metric graphs.

The package is available to install via the repository https://github.com/davidbolin/MetricGraph.

The following sections describe the main functionality of the package and summarizes some of the required theory. Section 2 introduces metric graphs and the metric_graph class, Section 3 shows how to work with random fields on metric graphs, and Section 4 introduces the inlabru interface of the package through an application to real data. For a complete introduction to the functionality of the package, we refer to the Vignettes available on the package homepage https://davidbolin.github.io/MetricGraph/. In particular, that contains an introduction to the INLA interface, the implementation of Whittle–Matérn fields with general smoothness, and further details and examples for all methods in the package.

The metric_graph class

A compact metric graph $\Gamma$ consists of a set of finitely many vertices $\mathcal{V}=\{v_i\}$ and a finite set $\mathcal{E}=\{e_j\}$ of edges connecting the vertices. Each edge $e$ is defined by a pair of vertices $(v_i,v_k)$ and a finite length $l_e \in (0,\infty)$. An edge in the graph is a curve parametrized by arc-length, and a location $s\in \Gamma$ is a position on an edge, and can thus be represented as a touple $(e,t)$ where $t\in[0,l_e]$.

Basic constructions

A metric graph is represented in the MetricGraph package through the class metric_graph. An object of this class can be constructed in two ways. The first is to specify the vertex matrix V and the edge matrix E, and it is then assumed that all edges are straight lines. The second, more flexible, option is to specify the object from a SpatialLines object using the sp package (Bivand, Pebesma, and Gomez-Rubio 2013).

To illustrate this, we use the osmdata package to download data from OpenStreetMap. In the following code, we extract the streets on the campus of King Abdullah University of Science and Technology (KAUST) as a SpatialLines object:

call <- opq(bbox = c(39.0884, 22.33, 39.115, 22.3056))
call <- add_osm_feature(call, key = "highway",value=c("motorway",
                                                      "primary","secondary",
                                                      "tertiary",
                                                      "residential"))
data <- osmdata_sf(call)

We can now create the metric graph as follows.

graph <- metric_graph$new(data)
graph$plot(vertex_size = 0.5)

Let us also build an interactive plot with the underlying map:

graph$plot(vertex_size = 0.5, type = "mapview")

Observe that when we create this graph from osmdata, it also includes the edge data as edge weights:

graph$get_edge_weights()

## # A tibble: 368 × 11
##    osm_id   name    bridge highway layer lit   maxspeed maxweight oneway surface
##    <chr>    <chr>   <chr>  <chr>   <chr> <chr> <chr>    <chr>     <chr>  <chr>  
##  1 38085402 Discov… NA     reside… NA    yes   40       NA        NA     NA     
##  2 38085409 King A… NA     primary NA    yes   60       NA        yes    NA     
##  3 39425304 King A… NA     primary NA    NA    20       NA        yes    NA     
##  4 39425484 Peace … NA     reside… NA    NA    40       NA        NA     NA     
##  5 39425550 Island… NA     reside… NA    NA    40       NA        NA     NA     
##  6 39425743 Coral … NA     reside… NA    yes   40       NA        NA     NA     
##  7 39425802 Hammou… NA     reside… NA    yes   40       NA        NA     NA     
##  8 39426595 Najel … NA     reside… NA    yes   40       NA        yes    NA     
##  9 39478848 Discov… NA     reside… NA    NA    40       NA        NA     NA     
## 10 39478849 Explor… NA     reside… NA    yes   40       NA        NA     NA     
## # ℹ 358 more rows
## # ℹ 1 more variable: .weights <dbl>

Now, note that when creating the graph there was a warning saying that the graph is not connected, so let us create a graph_components object that contains all connected components as graphs and then extract the largest connected component to work with

graphs <- graph_components$new(edges = data)
graph <- graphs$get_largest()
graph$plot(vertex_size = 0)

Here is the interactive plot:

graph$plot(vertex_size = 0, type = "mapview")