Building a mesh in INLA

Introduction

This particular vignette was written to get users started on how to build a Delaunay triangulation (or a mesh) that will serve as a basis to construct the spatial component of the models through an SPDE approach.

Properly choosing the structure of the mesh is essential for the model estimation to be carried out properly. Detailed explanation on the dos and don’ts of constructing a mesh are presented in section 1.3 of the R-INLA tutorial on SPDE models.

Without repeating all the details presented in the above mentioned tutorial, this vignette is designed to get users started. What is important to consider when constructing the mesh is that the triangles in the mesh be roughly of the same size. In addition, it is a good idea to include a buffer zone around the sampling area to prevent problems related to edge effect that may arise. This may be done explicitly or implicitly depending on the way the mesh is contructed.

For the moment, this is essentially to get access to the mite data that will be used in this illustration. In addition, by loading the mapSpecies R package, the INLA R package is also loaded.

library(INLA)
library(raster)

Data

Species data

# Load species data
load("assets/data/mite.spdf.rda")

Environmental data

# Environmental data
load("assets/data/mite.envRaster.rda")

Building the mesh

For our illustration, let’s first build a SpatialPolygons outlining the sampling area. Note that this step is not necessary to construct the mesh but it often the case that a SpatialPolygons object is available to delineate the study area.

poly <- Polygon(matrix(c(0,0,2.6,2.6,0,0,10,10,0,0), 
                       nrow = 5, ncol = 2))
spacePoly <- SpatialPolygons(list(Polygons(list(poly), ID = 1)),
                             proj4string = crs(mite.envRaster))

Using spacePoly, let’s construct the mesh. Usually, constructing the mesh requires a few trials to get exactly what is needed for the species understudy. So, do not despair!

The function used to construct the mesh comes directly from INLA. If you dig deep into INLA, you will find that there existe a few functions to construct meshes. Here we will used inla.mesh.2d. At this point it becomes important to read the help file of inla.mesh.2d to fully understand the meaning of each argument that can be tweaked to construct the mesh you want to have.

Here are a few notes on some of the arguments of inla.mesh.2d that are worth knowing:

max.edge : It is a length given in the same units as the coordinates of the studied system.

cutoff : It is a length given in the same units as the coordinates of the studied system.

Here are a few example of mesh that can be drawn for our study site.

Using a SpatialPolygons as boundaries

Using only max.edge and offset to define the mesh, it is possible to obtain a reasonably good mesh.

meshPoly1 <- inla.mesh.2d(boundary = spacePoly, max.edge = 2,
                          offset = c(2,2), cutoff = 0.5,
                          crs = crs(mite.envRaster))
meshPoly2 <- inla.mesh.2d(boundary = spacePoly, max.edge = 2,
                          offset = c(1,1), cutoff = 0.5,
                          crs = crs(mite.envRaster))
meshPoly3 <- inla.mesh.2d(boundary = spacePoly, max.edge = 0.5,
                          offset = c(2,2), cutoff = 0.5,
                          crs = crs(mite.envRaster))
meshPoly4 <- inla.mesh.2d(boundary = spacePoly, max.edge = 0.5,
                          offset = c(1,1), cutoff = 0.5,
                          crs = crs(mite.envRaster))

par(mfrow = c(2,2), mar = c(1,1,1,1))

plot(meshPoly1, main = "", asp = 1)
points(coordinates(mite.spdf), pch = 19, col = "red")

plot(meshPoly2, main = "", asp = 1)
points(coordinates(mite.spdf), pch = 19, col = "red")

plot(meshPoly3, main = "", asp = 1)
points(coordinates(mite.spdf), pch = 19, col = "red")

plot(meshPoly4, main = "", asp = 1)
points(coordinates(mite.spdf), pch = 19, col = "red")

If computational time is an issue, it is useful to compare the number of edges in the mesh. The larger the number of edges the longer it will take for model estimation. This can be checked with the following code.

meshPoly1$n

## [1] 82

meshPoly2$n

## [1] 58

meshPoly3$n

## [1] 921

meshPoly4$n

## [1] 590

Note that meshPoly4 is probably a good candidate at this point because it is a fine enough mesh to efficiently capture the spatial variation in the data. In addition, it has a reasonable number of edges, although if it turned out to be too long to run the analysis using such a mesh, it might be worth it to increase the max.edge values and thus reduced the number of edges in the analysis.

Using the sampled locations as a basis

This can be done using the argument loc

meshLoc1 <- inla.mesh.2d(loc = mite.spdf, max.edge = 2,
                         offset = c(1,1), cutoff = 0.5,
                         crs = crs(mite.envRaster))
meshLoc2 <- inla.mesh.2d(loc = mite.spdf, max.edge = 0.5,
                         offset = c(1,1), cutoff = 0.5,
                         crs = crs(mite.envRaster))

par(mfrow = c(1,2), mar = c(1,1,1,1))

plot(meshLoc1, main = "", asp = 1)
points(coordinates(mite.spdf), pch = 19, col = "red")

plot(meshLoc2, main = "", asp = 1)
points(coordinates(mite.spdf), pch = 19, col = "red")

# Number of edges in the mesh
meshLoc1$n

## [1] 98

meshLoc2$n

## [1] 890

or the argument loc.domain

meshDom1 <- inla.mesh.2d(loc.domain = mite.spdf, max.edge = 2,
                         offset = c(1,1), cutoff = 0.5,
                         crs = crs(mite.envRaster))
meshDom2 <- inla.mesh.2d(loc.domain = mite.spdf, max.edge = 0.6,
                         offset = c(1,1), cutoff = 0.5,
                         crs = crs(mite.envRaster))
meshDom3 <- inla.mesh.2d(loc.domain = mite.spdf, max.edge = 0.5,
                         offset = c(0.5,0.5), cutoff = 0.5,
                         crs = crs(mite.envRaster))

par(mfrow = c(1,3), mar = c(1,1,1,1))

plot(meshDom1, main = "", asp = 1)
points(coordinates(mite.spdf), pch = 19, col = "red")

plot(meshDom2, main = "", asp = 1)
points(coordinates(mite.spdf), pch = 19, col = "red")

plot(meshDom3, main = "", asp = 1)
points(coordinates(mite.spdf), pch = 19, col = "red")

# Number of edges in the mesh
meshDom1$n

## [1] 77

meshDom2$n

## [1] 608

meshDom3$n

## [1] 596

The main difference between using loc and loc.domain is that when using loc the samples are at the mesh edges which is not the case when loc.domain is used. In this case, meshDom3 is the best mesh because it has a wide enough boundary and the triangle size is roughly the same across the full area. In this respect, meshLoc1 should never be used because the triangles size differs importantly across the sampling region.