MetricGraph: Random Fields on Metric Graphs
David Bolin, Alexandre B. Simas, and Jonas Wallin
Created: 2022-11-26. Last modified: 2024-12-13.
Source:vignettes/MetricGraph.Rmd
MetricGraph.Rmd
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 consists of a set of finitely many vertices and a finite set of edges connecting the vertices. Each edge is defined by a pair of vertices and a finite length . An edge in the graph is a curve parametrized by arc-length, and a location is a position on an edge, and can thus be represented as a touple where .
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")