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{array}{r l} \int_X f(x)dx& =\int_X\exp[\log(f(x))]dx\\ &\approx \int_X \exp\left( \log(f(x_0)) + (x-x_0)\left.\frac{\partial \log(f(x))}{\partial x}\right\vert_{x=x_0} + \frac{(x-x_0)^2}{2}\left.\frac{\partial^2 \log(f(x))}{\partial x}^2\right\vert_{x=x_0}\right)dx, \end{array} $$

where the function $f(x)$ will be evaluated at the mode, $\left.f(x)\right\vert_{x=x_0}$, such that

$$x_0=\{x:\frac{\partial f(x)}{\partial x}=0 \wedge \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

$$\sigma^2=\left.\frac{1}{\partial^2 \log(f(x))/\partial x^2}\right\vert_{x=x_0},$$

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

$$\int_Xf(x)dx\approx f(x_0)\cdot \int_X \exp\left(-\frac{(x-x_0)^2}{2\sigma^2}\right)dx.$$

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 $\mathbf{x}\in \mathbb{R}^n$ is said GMRF with reference to a graph $\mathcal{G}=(\mathcal{V}, \mathcal{E})$ with mean $\boldsymbol\mu$ and precision matrix (symmetric positive difenite) $\mathbf{Q}>0$ if its density has the following structure

$$\pi(\mathbf{x})=(2\pi)^{-n/2}|\mathbf{Q}|^{1/2}\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol\mu)^T\mathbf{Q}(\mathbf{x}-\boldsymbol\mu)\right),$$

and

$$(\forall i\neq j) \left\lbrace Q_{ij}\neq 0 \iff \{i,j\}\in \mathcal{E}\right\rbrace.$$

If the precision matrix $\mathbf{Q}$ is completely dense then the network $\mathcal{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 $\mathbf{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 $\mathbf{x}$ is a GMRF with respect to a graph $\mathcal{G}=(\mathcal{V}, \mathcal{E})$, with mean $\boldsymbol\mu$ and SPD precision matrix $\mathbf{Q}>0$, then

$$\begin{array}{rcl} \text{E}(x_i|\mathbf{x}_{-i}) & = & \mu_i - \frac{1}{Q_{ii}}\sum_{j:j\sim i}Q_{ij}(x_j-\mu_j),\\ \text{Prec}(x_i|\mathbf{x}_{-i}) & = & Q_{ii},\\ \text{Corr}(x_i,x_j|\mathbf{x}_{-ij}) & = & -\frac{Q_{ij}}{\sqrt{Q_{ii}Q_{jj}}},\quad i\neq j.\\ \end{array} $$

The diagonal elements of $\mathbf{Q}$ are the conditional precisions of $x_i$ given $\mathbf{x}_{-i}$, while the off-diagonal elements, with a scaling factor, are the conditional correlation between $x_i$ and $x_j$, given $\mathbf{x}_{-ij}$.

Gaussian Latent Fields

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

$$\begin{array}{rcl} \mathbf{y}|\mathbf{\mathcal{X}},\boldsymbol\theta_1 & \sim & \prod_{i=1}^{n}\pi(y_i|\mathcal{X}_i,\boldsymbol\theta_1),\\ \mathbf{\mathcal{X}}|\boldsymbol\theta_2 & \sim & N(\mathbf{0},\mathbf{Q}_{\mathbf{\mathcal{X}}}^{-1}(\boldsymbol\theta_2)),\\ \boldsymbol\theta=\{\boldsymbol\theta_1,\boldsymbol\theta_2\} & \sim & \pi(\boldsymbol\theta),\\ \end{array} $$

where $\mathbf{y}|\mathbf{\mathcal{X}},\boldsymbol\theta_1$ is the data (or likelihood) level, where $\mathbf{\mathcal{X}}$ are the elements of the latent field and $\boldsymbol\theta_1$ are the hyperparameters of the likelihood. The elements of the latent field are distributed according to $\mathbf{\mathcal{X}}|\boldsymbol\theta_2$, following a GMRF with mean ${0}$ and the latent field structure is integrated into the structure of the precision matrix $\mathbf{Q}_{\mathbf{\mathcal{X}}}^{-1}(\boldsymbol\theta_2)$, where $\boldsymbol\theta_2$ are the hyperparameters of the latent field. Finally, the last level is the one concerning the distribution of the hyperparameters of the model $(\boldsymbol\theta)$, comprising both those of the likelihood $(\boldsymbol\theta_1)$ and those concerning the latent field $(\boldsymbol\theta_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{array}{c} \boldsymbol\eta=\beta_0\mathbf{1} + \boldsymbol\beta\mathbf{X} + \sum_{k=1}^K f_k(\mathbf{u}_k), \end{array} $$

where $(\beta_0, \boldsymbol\beta)$ are the parameters associated with the linear effects, while $\{\mathbf{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 $\boldsymbol\eta=\mathbf{A}_j\mathbf{u}_j$, where the effects $(\mathbf{u}_j)$ are linked to the predictor of the observations $(\boldsymbol\eta)$ through a projection matrix $(\mathbf{A}_j)$ relative to each effect $(\mathbf{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

$$ \boldsymbol\eta=\left( \begin{array}{c} \eta_1\\ \hline \vdots \\ \hline \eta_K \end{array}\right)=\left(\mathbf{A}_1| \mathbf{A}_2| \cdots| \mathbf{A}_J\right) \left( \begin{array}{c} \mathbf{u}_1\\ \hline \mathbf{u}_2 \\ \hline \vdots \\ \hline \mathbf{u}_J \end{array}\right), $$

where each effect $\mathbf{u}_j:j\in (1,...,J)$ is related to its corresponding projection matrix $\mathbf{A}_j$, linking the effects with the linear predictor $\boldsymbol\eta=(\eta_1,...,\eta_K)$.

Key Articles

  1. 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
  2. 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
  3. 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
  4. 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
  5. 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

  1. 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
  2. Gómez-Rubio, V. (2020). Bayesian Inference with INLA. In Bayesian Inference with INLA. Chapman & Hall/CRC Press. https://doi.org/10.1201/9781315175584
  3. Moraga, P. (2019). Geospatial Health Data. Chapman and Hall/CRC. https://doi.org/10.1201/9780429341823
  4. 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
  5. Rue, H., & Held, L. (2005). Gaussian Markov Random Fields. Chapman and Hall/CRC. https://doi.org/10.1201/9780203492024
  6. Xiaofeng Wang, Ryan Yue, & Faraway, J. J. (2018). Bayesian Regression Modeling with INLA. Chapman & Hall.

Subsections of INLA

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.

Subsections of Basic Operations

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).

Subsections of INLA's extensions

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

INLA for Fitting Dirichlet Regression Models (Summary)

References

  1. 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
  2. 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

Subsections of dirinla

Example with fixed effects

Example with iid (unstructured diagonal) effects

Example with spatial structured random effects

Subsections of Spatial Models (SPDE-FEM)

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.

Subsections of Various Functions

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
INLA_mesh_with_non_convex_boundary.png

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.
Mathematical reasoning.

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
Polygons_with_holes
Example with nested polygons.