INLA methodology

There are various methodologies for the computational implementation of Bayesian inference, simulation methods such as MCMC (Markov chain Monte Carlo), or approximate methods like VB (variational Bayes); all of them with their own challenges. However, INLA (Integrated Nested Laplace Approximation) is a deterministic approximate approach, developed by (Rue, Martino y Chopin, 2009) and expanded upon in (Lindgren y Rue, 2015; Rue y col., 2017; Bakka y col., 2018). . It allows for Bayesian inference in a set of structured additive models, termed latent Gaussian models (LGMs). The INLA method enables the calculation of joint posterior distributions, the marginal distributions of each parameter and hyperparameter, as well as combinations of these or the posterior predictive distributions.

At the core of INLA is the Laplace approximation applied to the expression of the conditional probability distribution of the latent field. This implies that the latent structure must follow a Gaussian Markov Random Field (GMRF) that can be linked to latent Gaussian models (Lindgren, Rue y Lindström). Although many models can be rewritten in such a way that their structure is similar to an LGM.

Laplace Approximation

The Laplace approximation for a density function $f (x)$ involves transformation using logarithms and carrying out a second-order Taylor series expansion, evaluated at the mode of the function:

\begin{aligned} \int_{X} f (x) d x & = \int_{X} \exp [\log (f (x))] d x \\ \approx \int_{X} \exp (\log (f (x_{0})) + (x - x_{0}) {\frac{\partial \log (f (x))}{\partial x} |}_{x = x_{0}} + \frac{(x - x_{0})^{2}}{2} {{\frac{\partial^{2} \log (f (x))}{\partial x}}^{2} |}_{x = x_{0}}) d x, \end{aligned}

where the function $f (x)$ will be evaluated at the mode, ${f (x) |}_{x = x_{0}}$ , such that

x_{0} = {x : \frac{\partial f (x)}{\partial x} = 0 \land \frac{\partial^{2} f (x)}{\partial x^{2}} \neq 0} .

That is, the function is evaluated when the first derivative is null, so the first-order term in the Taylor series expansion can be simplified. Also, if we express the second-order term as

σ^{2} = {\frac{1}{\partial^{2} \log (f (x)) / \partial x^{2}} |}_{x = x_{0}},

then we can express the Laplace approximation as the kernel of a Gaussian function:

\int_{X} f (x) d x \approx f (x_{0}) \cdot \int_{X} \exp (- \frac{(x - x_{0})^{2}}{2 σ^{2}}) d x .

Gaussian Markov Random Field

A Gaussian MArkov Random Field (GMRF) is a Gaussian Field (GF) with Markov properties. This means that, given a random vector $x \in R^{n}$ is said GMRF with reference to a graph $G = (V, E)$ with mean $μ$ and precision matrix (symmetric positive difenite) $Q > 0$ if its density has the following structure

π (x) = (2 π)^{- n / 2} | Q |^{1 / 2} \exp (- \frac{1}{2} (x - μ)^{T} Q (x - μ)),

and

(\forall i \neq j) {Q_{i j} \neq 0 ⟺ {i, j} \in E} .

If the precision matrix $Q$ is completely dense then the network $G$ is completely connected. This implies that any normal distribution with symmetric positive definite (SPD) covariance matrixes is a GMRF and vice versa.

In the case where $Q$ is sparse then the properties of GMRFs are really useful and we can make use of them. In particular, a very useful property is the interpretation of the conditional distributions of the elements of a GMRF.

Suppose $x$ is a GMRF with respect to a graph $G = (V, E)$ , with mean $μ$ and SPD precision matrix $Q > 0$ , then

\begin{array}{rcl} E (x_{i} | x_{- i}) & = & μ_{i} - \frac{1}{Q_{i i}} \sum_{j : j \sim i} Q_{i j} (x_{j} - μ_{j}), \\ Prec (x_{i} | x_{- i}) & = & Q_{i i}, \\ Corr (x_{i}, x_{j} | x_{- i j}) & = & - \frac{Q_{i j}}{\sqrt{Q_{i i} Q_{j j}}}, i \neq j . \end{array}

The diagonal elements of $Q$ are the conditional precisions of $x_{i}$ given $x_{- i}$ , while the off-diagonal elements, with a scaling factor, are the conditional correlation between $x_{i}$ and $x_{j}$ , given $x_{- i j}$ .

Gaussian Latent Fields

The structure on which INLA is based can be summarised in the following hierarchical model:

\begin{array}{rcl} y | X, θ_{1} & \sim & \prod_{i = 1}^{n} π (y_{i} | X_{i}, θ_{1}), \\ X | θ_{2} & \sim & N (0, Q_{X}^{- 1} (θ_{2})), \\ θ = {θ_{1}, θ_{2}} & \sim & π (θ), \end{array}

where $y | X, θ_{1}$ is the data (or likelihood) level, where $X$ are the elements of the latent field and $θ_{1}$ are the hyperparameters of the likelihood. The elements of the latent field are distributed according to $X | θ_{2}$ , following a GMRF with mean $0$ and the latent field structure is integrated into the structure of the precision matrix $Q_{X}^{- 1} (θ_{2})$ , where $θ_{2}$ are the hyperparameters of the latent field. Finally, the last level is the one concerning the distribution of the hyperparameters of the model $(θ)$ , comprising both those of the likelihood $(θ_{1})$ and those concerning the latent field $(θ_{2})$ .

The second level is the latent Gaussian field, which constitutes a latent Gaussian model (LGM). LGMs are a class of models that follow Gaussian processes, be it for time series, spatial models, iid random effects, cluster random effects, etc. Therefore, the Gaussian field that has the above structure can also be formulated according to the linear predictor of the model as

\begin{matrix} η = β_{0} 1 + β X + \sum_{k = 1}^{K} f_{k} (u_{k}), \end{matrix}

where $(β_{0}, β)$ are the parameters associated with the linear effects, while ${f}$ are the unknown functions of the random effects $\mathbf{U}=\{\mathbf{u}\}, which can have very different structures: iid, random walks, Besag, SPDE (Stochastic Partial Differential Equation), etc.

Based on the above expression, we can reformulate it in matrix terms $η = A_{j} u_{j}$ , where the effects $(u_{j})$ are linked to the predictor of the observations $(η)$ through a projection matrix $(A_{j})$ relative to each effect $(u_{j})$ . This projection matrix integrates the weights associated with the effects, i.e. the values of the explanatory variables for linear effects, smoothing weights or associated weights with SPDEs. That is, we can rewrite it as

η = (\begin{matrix} η_{1} \\ ⋮ \\ η_{K} \end{matrix}) = (A_{1} | A_{2} | \dots | A_{J}) (\begin{matrix} u_{1} \\ u_{2} \\ ⋮ \\ u_{J} \end{matrix}),

where each effect $u_{j} : j \in (1, . . ., J)$ is related to its corresponding projection matrix $A_{j}$ , linking the effects with the linear predictor $η = (η_{1}, . . ., η_{K})$ .

Key Articles

Bakka, H., Rue, H., Fuglstad, G.-A., Riebler, A. I., Bolin, D., Illian, J., Krainski, E., Simpson, D. P., & Lindgren, F. K. (2018). Spatial modelling with INLA: A review. In Wires (Vol. xx, Issue Feb). https://doi.org/10.1002/wics.1443
Lindgren, F., Rue, H., & Lindström, J. (2011). An explicit link between gaussian fields and gaussian markov random fields: The stochastic partial differential equation approach. Journal of the Royal Statistical Society. Series B: Statistical Methodology, 73(4). https://doi.org/10.1111/j.1467-9868.2011.00777.x
Lindgren, F., & Rue, H. (2015). Bayesian spatial modelling with R-INLA. Journal of Statistical Software, 63(19). https://doi.org/10.18637/jss.v063.i19
Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society. Series B: Statistical Methodology, 71(2). https://doi.org/10.1111/j.1467-9868.2008.00700.x
Rue, H., Riebler, A., Sørbye, S., Illian, J., Simpson, D. & Lindgren, F. (2017). Bayesian Computing with INLA: A Review. Annual Review of Statistics and Its Application, 4:1, 395-421. https://doi.org/10.1146/annurev-statistics-060116-054045

Books

Blangiardo, M., & Cameletti, M. (2015). Spatial and Spatio-temporal Bayesian Models with R - INLA. In Spatial and Spatio-temporal Bayesian Models with R - INLA. Wiley. https://doi.org/10.1002/9781118950203
Gómez-Rubio, V. (2020). Bayesian Inference with INLA. In Bayesian Inference with INLA. Chapman & Hall/CRC Press. https://doi.org/10.1201/9781315175584
Moraga, P. (2019). Geospatial Health Data. Chapman and Hall/CRC. https://doi.org/10.1201/9780429341823
Krainski, E., Gómez-Rubio, V., Bakka, H., Lenzi, A., Castro-Camilo, D., Simpson, D., Lindgren, F., & Rue, H. (2018). Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA. In Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA. https://doi.org/10.1201/9780429031892
Rue, H., & Held, L. (2005). Gaussian Markov Random Fields. Chapman and Hall/CRC. https://doi.org/10.1201/9780203492024
Xiaofeng Wang, Ryan Yue, & Faraway, J. J. (2018). Bayesian Regression Modeling with INLA. Chapman & Hall.

Basic Operations

In this part we will show different basic functions and operations, which are the basis to modeling with INLA.

Basic Linear Models
Simple examples of linear models.
Stacks
A short introduction to inla stacks and their uses.
Mesh Construction
Brief introduction to the construction of meshes: applications and operational basis.
Advanced Features
A short introduction to INLA advanced features and their uses.

Basic Linear Models

In this section we will show

Stacks

Mesh Construction

Advanced Features

INLA's extensions

In this section we will show

dirinla
Examples using the dirinla R package for Compositional Data (CoDa).

dirinla

Esta sección está dedicada específicamente al paquete dirinla que está construido sobre el enfoque implementado en R-INLA. Para que los ejemplos resulten comprensibles dedicaremos esta parte introductoria para sintetizar los fundamentos matemáticos y operativos de dirinla. Esto es, realizaremos una síntesis del artículo Martínez-Minaya et al. (2023) en el que se presenta el paquete y también se evaluarán algunas funciones que integran el paquete, modificándolas allí donde se considere pertinente.

Examples

Example with fixed effects
Simulated Dirichlet data fitted with a fixed effects model.
Example with iid (unstructured diagonal) effects
Simulated Dirichlet data fitted with iid (unstructured diagonal) random effects model.
Example with spatial structured random effects
Simulated Dirichlet data fitted with spatial structured random effects model.

INLA for Fitting Dirichlet Regression Models (Summary)

References

Joaquín Martínez-Minaya, Finn Lindgren, Antonio López-Quílez, Daniel Simpson & David Conesa (2023) The Integrated Nested Laplace Approximation for Fitting Dirichlet Regression Models, Journal of Computational and Graphical Statistics, 32:3, 805-823, DOI: https://doi.org/10.1080/10618600.2022.2144330
Depaoli, S., Clifton, J. P., & Cobb, P. R. (2016). Just Another Gibbs Sampler (JAGS): Flexible Software for MCMC Implementation. Journal of Educational and Behavioral Statistics, 41(6), 628–649. http://www.jstor.org/stable/26447820

Example with fixed effects

Example with iid (unstructured diagonal) effects

Example with spatial structured random effects

Spatial Models (SPDE-FEM)

In this section we will show

Geostatistical Model
A simple geostatical model in INLA.
Preferential Model
A joint model for preferentially sampled data.
SVC Model
Spatially Varying Coefficient Model.

Geostatistical Model

In this section we will show

Preferential Model

In this section we will show

SVC Model

In this section we will show

Various Functions

Here we will present different functions that have been developed for implementation in automated processes or that may be of interest.

Transformation to sp/sf
Here we define some functions related to the transformations of inla.mesh objects to sp and sf.

Transformation to sp/sf

Version infromation

In this we expose a list with all the versions of the programms, operative system and the libraries used in this section code, which almost is extracted from sessionInfo().

SO (version): Windows 11 Home (22H2).
R: 4.2.2 (2022-10-31 ucrt).
RStudio: 2022.07.2 (Build 576).
Libraries:
- INLA: INLA_22.05.07
- inlabru: inlabru_2.5.3
- sp: sp_1.5-0
- sf: sf_1.0-8
- ggplot2: ggplot2_3.3.6
- viridis: viridis_0.6.2

There may be cases when you need to work with inla elements as if they were spatial objects (usually sp or sf), this will be particularly appropriate when working with meshes or mesh edges. Whatever the reasons for passing the mesh or mesh edges to spatial objects, there is no built-in function available in library(INLA), so we will proceed with its definition.

Mesh to spatial object

Suppose we start from a set of points given by the XY data frame, on which we apply the function inla.nonconvex.hull to obtain the boundary of our mesh, resulting from the morphological operations of dilation and closure incorporated for the formation of the outer edge.

Get XY data

library(INLA)
library(inlabru)
XY <- readRDS(file="XY.rds")
INLAnonconvex <- inla.nonconvex.hull(XY, convex=1000, resolution=400)
INLAnonconvexMesh <- inla.mesh.2d(boundary=INLAnonconvex, max.edge=5000)

If we plot the mesh ggplot() + inlabru::gg(INLAnonconvexMesh) + theme_void() we obtain the following figure:

Fig. 1: INLA mesh with non convex boundary

Once the mesh is available, to transform the tessellation into a spatial object (sp) we use the information integrated in the mesh object (mesh). Specifically in the rows of the mesh$graph$tv matrix, which store the vertex positions for each i-th triangle. In this way, we simply inscribe the vertex indices for each triangle in the coordinate matrix of the vertices themselves, forming a polygon with each one and finally grouping the result in a SpatialPolygons object, or in a This is synthesized in the following function:

SpatialPolygonDelaunay <- function(mesh){
  return(SpatialPolygons(sapply(1:nrow(mesh$graph$tv), function(i){
    list(Polygons(list(Polygon(mesh$loc[mesh$graph$tv[i,],1:2])),i
    ))
  })
  ))
}

SpatialPolygonDelaunay <- function(mesh){
  return(st_polygon(sapply(1:nrow(mesh$graph$tv), function(i){
    list(mesh$loc[c(mesh$graph$tv[i,], mesh$graph$tv[i,1]),1:2])
  })
  ))
}

Such that the object transformed would be given by:

SpatialPolygonMesh <- SpatialPolygonDelaunay(mesh=INLAnonconvexMesh)

Boundary to spatial object

In this case we intend to transform the edges of the mesh to SpatialPolygons, in such a way that we can account for the gaps present in the mesh. So, the first thing to do is to convert the edges into simply connected polygons (without gaps).

To correctly conform a polygon, the order of succession of the points on the perimeter of the polygon must be indicated, such that the first point is also the last point to close the polygon. This can be easily done because the mesh stores the locations of the vertices that make up the edges as pairs of points in the rows of the mesh$segm$bnd$idx matrix. That is, if the points refer to the same polygon we will have the connection established from the pairs of indexes going through the rows of this matrix, such that it will come to repeat in the first column the position of the initial point to close the polygon.

Therefore, knowing that the sequence of indices for each closed curve implies that the last pair ends at the first point (first pair of indices), we will have that the element of the first column of the next row will be different from the element of the second column of the reference row: mesh$segm$bnd$idx[i+1,1]!=mesh$segm$bnd$idx[i,1]. So, through this condition, we can set the indices of the polygon boundaries that we want to reconstruct as spatial objects (sp or sf).

SpatialPolygonsBoundary <- function(mesh){
  BoundaryIndx <- sapply(1:(nrow(mesh$segm$bnd$idx)-1), function(i){
    return(mesh$segm$bnd$idx[i+1,1]==mesh$segm$bnd$idx[i,2])
  })
  BoundaryIndxLims <- which(!BoundaryIndx)
  BoundaryIndxLims <- c(0, BoundaryIndxLims, nrow(mesh$segm$bnd$idx))
  SpatialPolygonsBoundary <- SpatialPolygons(sapply(2:(length(BoundaryIndxLims)), function(i){
    VertixBoundaryIndex <- matrix(c(mesh$segm$bnd$idx[(BoundaryIndxLims[i-1]+1):BoundaryIndxLims[i],1],
                                    mesh$segm$bnd$idx[BoundaryIndxLims[i],2]), 
                                  ncol=1, byrow=TRUE)
    return(list(Polygons(list(Polygon(mesh$loc[VertixBoundaryIndex,1:2])),i)))
  }))
  return(SpatialPolygonsBoundary)
}

SpatialPolygonsBoundary <- function(mesh){
  BoundaryIndx <- sapply(1:(nrow(mesh$segm$bnd$idx)-1), function(i){
    return(mesh$segm$bnd$idx[i+1,1]==mesh$segm$bnd$idx[i,2])
  })
  BoundaryIndxLims <- which(!BoundaryIndx)
  BoundaryIndxLims <- c(0, BoundaryIndxLims, nrow(mesh$segm$bnd$idx))
  SpatialPolygonsBoundary <- st_polygon(sapply(2:(length(BoundaryIndxLims)), function(i){
    VertixBoundaryIndex <- c(mesh$segm$bnd$idx[(BoundaryIndxLims[i-1]+1),1], 
                             mesh$segm$bnd$idx[(BoundaryIndxLims[i-1]+1):BoundaryIndxLims[i],2])
    return(list(cbind(mesh$loc[VertixBoundaryIndex,1:2])))
  }))
  return(SpatialPolygonsBoundary)
}

Then, once all the polygons have been formed from the edges given by the INLA function, these polygons must be transformed to account for the gaps they may contain and are present. To do this we can either perform a mathematical approach to the properties of the problem, using set theory for ordered sets, or we can approach it operationally to solve it.

Operational reasoning.

To do this according to an operational approach we must consider that we have formed a polygon for each closed curve, for each edge or boundary, so that these curves may contain within them other surfaces given by other curves. In addition, some of these polygons may be associated with the inner edge of the polygons, which would refer to the boundary of the inner hole of the polygon, so that we can have to classify the polygons obtained in two classes or types: full polygons and hollow polygons, depending on whether they refer to the border of a “full” area or a “hollow” area. In addition, it must be taken into account that a filled polygon will contain in the next level of maximum internal intersection, in case of containing in its interior other polygons, a hollow polygon. While for the case of hollow polygons this condition assumes the counterpart; namely, that the next level of maximum internal intersection will be a filled polygon. Such that the level of maximal internal intersectionality with respect to a given polygon we shall understand as the polygon of largest area that is completely contained in the reference polygon. But we can extend this definition according to different sets of mutually independent polygons, where the maximum internal intersectionality for the reference polygon is no longer a single polygon, but the largest polygon of each of these mutually disjoint sets.

On the other hand, if a polygon, bounded by one of these given curves, has a non-zero intersection with others, then this intersection will totally reproduce one of the two polygons.

Subsequently, if analogously to an adjacency or neighborhood matrix, we elaborate an intersectionality matrix, determined according to the intersection operation of a given polygon on the rest of polygons. Such that the matrix will contain the relations between a given polygon with the whole set of polygons according to the intersection operation, operation that will be worth 1 if there is intersection and 0 if there is not. Therefore, for each row or column we will obtain the intersection of a polygon with the rest, where the diagonal of the matrix would manifest the intersection of a polygon with itself, so the main diagonal will be a vector of ones.

From the intersectionality matrix we can conclude that the rows or columns whose sum is unitary (only the element of the main diagonal will be different from 0 for that row or column), will belong to a polygon that will not contain any other polygon, whether it is hollow or full. Therefore, we could extract these polygons as correctly defined polygons that do not require any transformation.

Finally, in case a given row sums to more than one it will imply that it contains at least one hole. So to form the polygons with the holes we can approach the procedure in the following way:

First, for each row we can search for the smallest polygon (area) in the whole set of intersecting polygons, from which a completely ordered set will be constructed according to the intersections of the polygon.
Secondly, the set of polygons can be recomposed as the rearrangement into the smallest number of totally ordered sets. Such that from these totally ordered sets we will know that the largest of the sets refers to a filled polygon, the next to a hollow one, … and so on.
To thirdly, to conform by means of operations of difference and union the different polygons that account for the existing hollows.

Mathematical reasoning.

From a mathematical approach we can use set theory to evaluate the set of operations to be performed to solve the question. To do this, the first thing we must do is to establish the conditions of the problem:

We have a set of $n$ polygons $p_{i}$ , such that the universal set of reference will be $P = {p_{1}, \dots, p_{n}}$ . These polygons ar esubsets of the two dimensional Euclidian space, $p_{i} \subset R^{2}$ .
There are two kinds of polygons: (i) filled polygons $P_{F} = {p_{f i}}$ and (ii) hollow polygons $P_{H} = {p_{h i}}$ , by jointly reproducing the universal set of polygons $P = P_{F} \cup P_{H}$ .
There is a function $f$ that stablish the relation between a polygon $p_{i}$ with its area $a_{i}$ , such that $f : P ⟼ A$ implies $f (p_{i}) = a_{i}$ . However, this function is only surjective if we assume that there can be two or more polygons with the same area, while it would be surjective and injective if for each polygon it has a different area.
If the intersection of two polygons is non-zero, then the intersection results in the polygon of smaller area. That is, $p_{i} \cap p_{j} = {\emptyset ⊻ p_{j}} | f (p_{i}) > f (p_{j})$ . In addition, we have as a condition that there is no polygon intersecting with another that has the same area, $(∄ p_{i}) (p_{i}, p_{j} \in P) : f (p_{i}) = f (p_{j}) ⟺ {p_{i} \cap p_{j} \neq \emptyset \land i \neq j}$ .
For the set of all non-zero intersections of polygons $(p_{i}, p_{j})$ , we will have that the set with the largest area is a filled polygon $f (p_{i}) > f (p_{j}) : p_{i} \in P_{F}$ .
Finally, as a weak condition between the non-zero intersection of two filled polygons, we have that there exists, at least, a hollow polygon with greater area than that of the filled polygon of smaller area and that also intersects with them. In other words, if $p_{j} \subset p_{i} : p_{i}, p_{j} \in P_{F} \Rightarrow (\exists p_{k}) p_{k} \in P_{H} : p_{j} \subset p_{k} \subset p_{i}$ . While the strong condition will be, that if we have two polygons such that one is a subset of another with no filled polygon that can be incorporated as a middle element between the two by the subset operation, then there exists one and only one hollow polygon that lies between the two:

(\forall p_{i}, p_{j} \in P_{F}) p_{j} \subset p_{i} : ∄ p_{k} \in P_{F}, p_{j} \subset p_{k} \land p_{k} \subset p_{i} ⟺ \exists! p_{z} \in P_{H} : p_{j} \subset p_{z} \land p_{z} \subset p_{i} .

Given these conditions it is obvious that we can reorder the set $P$ of all the polygons, in such a way that we have the least number of completely ordered sets of higher cardinality whose union reproduces the set { $P$ } of reference. That is to say, we have to conform all those completely ordered sets $(P_{i}, \subset)$ , which implies that $(\forall p_{i}) p_{i} \in P_{i} (\exists p_{j}) p_{j} \in P_{i} \to p_{j} \subset p_{i} \lor p_{i} \subset p_{j}$ . Therefore, we shall construct these completely ordered sets $P_{i} = {p_{i 1} \subset p_{2 i} \subset \dots \subset p_{n i}}$ , of which we know from conditions (5) and (6) that for any of these sets $P_{i}$ , we will have that the filled polygons will be:

{p_{k i}, p_{(k - 2) i}, . . .,} \in P_{F} \Rightarrow p_{(k - 2 m) i} \in P_{F} : m \in M \subset N,

while the hollow polygons would be:

{p_{(k - 1) i}, p_{(k - 3) i}, . . .,} \in P_{H} \Rightarrow p_{(k - 1 - 2 m) i} \in P_{H} : m^{'} \in M^{'} \subset N .

The code to perform the pertinent transformations and obtain the polygons with the corresponding holes is as follows:

CompleteOrderedSetsFromSpatialPolygonsClass <- function(SPDINLAnonconvexBoundary){
  SPDINLABoundaryList <- list(Polygons=list(), 
                              DFIDArea=data.frame(matrix(ncol=2,nrow=0, dimnames=list(NULL, c("ID", "Area"))))
  )
  for(i in 1:length(SPDINLAnonconvexBoundary)){
    SPDINLABoundaryList$Polygons[[i]] <- SPDINLAnonconvexBoundary[i]
    SPDINLABoundaryList$DFIDArea[i,] <- c(i, SPDINLAnonconvexBoundary[i]@polygons[[1]]@area)
  }
  
  IndexIDV <- SPDINLABoundaryList$DFIDArea$ID
  
  ConnectivityMatrix <- matrix(ncol=length(SPDINLAnonconvexBoundary), nrow=length(SPDINLAnonconvexBoundary))
  for(i in 1:length(SPDINLAnonconvexBoundary)){
    for(j in i:length(SPDINLAnonconvexBoundary)){
      ConnectivityMatrix[i,j] <- ConnectivityMatrix[j,i] <- as.numeric(gIntersects(SPDINLAnonconvexBoundary[i],
                                                                                   SPDINLAnonconvexBoundary[j]))
    }
  }
  
  IndexIsolatedPolygonFill <- which(apply(ConnectivityMatrix, FUN=sum, MARGIN=1)==1) 

  IndexIDV <- if(length(IndexIsolatedPolygonFill)==0){IndexIDV} else{IndexIDV[-which(apply(ConnectivityMatrix, FUN=sum, MARGIN=1)==1)]}
  
  CompleteOrderedSets <- list(CompleteOrderedSets=list(), DFCompleteOrderedSets=list(), 
                              CompleteOrderedIsolatedSets=list(), DFCompleteOrderedIsolatedSets=list())
  
  if(length(IndexIsolatedPolygonFill)!=0){
    for(i in 1:length(IndexIsolatedPolygonFill)){
      CompleteOrderedSets$CompleteOrderedIsolatedSets[[i]] <- SPDINLABoundaryList$Polygons[[IndexIsolatedPolygonFill[i]]]
      CompleteOrderedSets$DFCompleteOrderedIsolatedSets[[i]] <- 
        data.frame(ID=SPDINLABoundaryList$Polygons[[IndexIsolatedPolygonFill[i]]]@polygons[[1]]@ID,
                   Area=SPDINLABoundaryList$Polygons[[IndexIsolatedPolygonFill[i]]]@polygons[[1]]@area,
                   HOLE=FALSE)
    }
  }
  
  if(length(IndexIsolatedPolygonFill)!=length(IndexIDV)){
    n_step <- 1
    while(length(IndexIDV)!=0){
      IDminArea <- SPDINLABoundaryList$DFIDArea[IndexIDV,"ID"][which.min(SPDINLABoundaryList$DFIDArea[IndexIDV,"Area"])] 
      IndxCOS <- which(as.logical(ConnectivityMatrix[IDminArea,]))
      DFIDAreaV <- SPDINLABoundaryList$DFIDArea[IndxCOS,]
      CompleteOrderedSets$CompleteOrderedSets[[n_step]] <- 
        SpatialPolygons(sapply(DFIDAreaV[order(DFIDAreaV$Area, decreasing=TRUE),"ID"],
                               function(i){
                                 list(SPDINLABoundaryList$Polygons[[i]]@polygons[[1]])
                               })
        )
      CompleteOrderedSets$DFCompleteOrderedSets[[n_step]] <- 
        data.frame(ID=DFIDAreaV[order(DFIDAreaV$Area, decreasing=TRUE),"ID"],
                   Area=DFIDAreaV[order(DFIDAreaV$Area, decreasing=TRUE),"Area"],
                   HOLE=floor(1:length(IndxCOS)/2)==1:length(IndxCOS)/2
        )
      
      IndexIDV <- setdiff(IndexIDV,IndxCOS)
      n_step <- n_step+1
    } 
    SpatialPolygonsNames <- c()
    for(i in 1:length(CompleteOrderedSets$CompleteOrderedSets)){
      if(length(CompleteOrderedSets$CompleteOrderedSets[[i]])/2!=floor(length(CompleteOrderedSets$CompleteOrderedSets[[i]])/2)){
        SpatialPolygonsNames <- 
          c(SpatialPolygonsNames, 
            paste0("SPDINLA", 
                   CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[length(CompleteOrderedSets$CompleteOrderedSets[[i]])]))
        
        assign(x=paste0("SPDINLA",
                        CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[length(CompleteOrderedSets$CompleteOrderedSets[[i]])]),
               value=SpatialPolygons(list(
                 CompleteOrderedSets$CompleteOrderedSets[[i]]@polygons[[length(CompleteOrderedSets$CompleteOrderedSets[[i]])]]
               )))
      }
      
      for(j in 1:floor(length(CompleteOrderedSets$CompleteOrderedSets[[i]])/2)){
        indx <- 1:2+2*(j-1)
        if(exists(paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]))){
          assign(x=paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]),
                 value=gDifference(eval(parse(text=paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]))), 
                                   SpatialPolygons(list(CompleteOrderedSets$CompleteOrderedSets[[i]]@polygons[[indx[2]]])),
                                   id=paste0(CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]])
                 )
          )
        } else{
          SpatialPolygonsNames <- c(SpatialPolygonsNames, paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]))
          assign(x=paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]),
                 value=gDifference(SpatialPolygons(list(CompleteOrderedSets$CompleteOrderedSets[[i]]@polygons[[indx[1]]])), 
                                   SpatialPolygons(list(CompleteOrderedSets$CompleteOrderedSets[[i]]@polygons[[indx[2]]])),
                                   id=paste0(CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]])
                 )
          )
        }}
      
      if(length(IndexIsolatedPolygonFill)!=0){
        SpatialPolygonRefined <- SpatialPolygons(c(lapply(1:length(SpatialPolygonsNames), function(i){
          return(eval(parse(text=paste0(SpatialPolygonsNames[i], "@polygons[[1]]"))))
        }), lapply(1:length(CompleteOrderedSets$CompleteOrderedIsolatedSets), function(i){
          return(CompleteOrderedSets$CompleteOrderedIsolatedSets[[i]]@polygons[[1]])
        })))        
      } else{
        SpatialPolygonRefined <- SpatialPolygons(lapply(1:length(SpatialPolygonsNames), function(i){
          return(eval(parse(text=paste0(SpatialPolygonsNames[i], "@polygons[[1]]"))))
        }))
      }

    }
  } else{
    SpatialPolygonRefined <- SpatialPolygons(lapply(1:length(CompleteOrderedSets$CompleteOrderedIsolatedSets), function(i){
      return(CompleteOrderedSets$CompleteOrderedIsolatedSets[[i]]@polygons[[1]])
    }))
  }
  
  SpatialPolygonRefined@proj4string <- CRS("+proj=utm")
  CompleteOrderedSetsAndPolygons <- list(CompleteOrderedSets=CompleteOrderedSets,
                                         SpatialPolygonRefined=SpatialPolygonRefined)
  
  return(CompleteOrderedSetsAndPolygons)
}

Thanks to this function we have all the polygons according to the closed curves that determine them. This allows us to perform operations on them, so that we are taking into account the existing gaps.

Therefore, following the example, we would have to apply the first function SpatialPolygonsBoundary on the result of the mesh INLAnonconvexMesh, and then apply the function just defined.

SPDINLAnonconvexBoundary <- SpatialPolygonsBoundary(mesh = INLAnonconvexMesh)
PolygonsRefined <- CompleteOrderedSetsFromSpatialPolygonsClass(SPDINLAnonconvexBoundary = SPDINLAnonconvexBoundary)

In order to visualize that we now have all the polygons available independent of each other, we can plot them in different colors:

col <- viridis::turbo(length(PolygonsRefined$SpatialPolygonRefined))
ggplotPolygons <- ggplot()
for(i in 1:length(PolygonsRefined$SpatialPolygonRefined)){
    ggplotPolygons <- ggplotPolygons + geom_sf(data=st_as_sf(PolygonsRefined$SpatialPolygonRefined[i]), fill=col[i], alpha=0.5)
}

So, if we represent the polygons as ggplotPolygons + theme_void() we will obtain the following figure:

Fig. 2: Refined polygons

Example with nested polygons.

INLA methodology

Laplace Approximation

Gaussian Markov Random Field

Gaussian Latent Fields

Key Articles

Books

Subsections of INLA

Basic Operations

Subsections of Basic Operations

Basic Linear Models

Stacks

Mesh Construction

Advanced Features

INLA's extensions

Subsections of INLA's extensions

dirinla

Examples

INLA for Fitting Dirichlet Regression Models (Summary)

References

Subsections of dirinla

Example with fixed effects

Example with iid (unstructured diagonal) effects

Example with spatial structured random effects

Spatial Models (SPDE-FEM)

Subsections of Spatial Models (SPDE-FEM)

Geostatistical Model

Preferential Model

SVC Model

Various Functions

Subsections of Various Functions

Transformation to sp/sf

Mesh to spatial object

Boundary to spatial object