Skip to contents

(Vignette under construction!)

Introduction

In traditional INLA code involving inla.mesh objects and inla.spde models, the inla.spde.make.A() function is used to construct the component design matrix that maps between spatial/spatio-temporal locations and the latent variables associated with mesh basis functions. For 2-manifold meshes, such as flat and spherical meshes, the implementation has one latent variable per mesh node, linked to piecewise linear basis functions on the mesh triangles. For 1-manifolds such as intervals and cyclic domains, both piecewise linear and piecewise quadratic basis functions are supported.

The inla.spde.make.A() interface supports a variety of features, that can be broken down in more simple building blocks. With the inlabru bru_mapper system, these building blocks are more easily customised for specific uses, some of which aren’t necessarily connected to spde models, such as the block feature that is used to aggregate rows of a design matrix, e.g. to construct numerical integration schemes.

If you haven’t already, go read the bru_mapper vignette for information about the various bru_mapper classes and methods. Then come back here to continue.

Converting the basic inla.spde.make.A calls into mappers

The most basic inla.spde.make.A call is to map purely spatial points to a mesh:

mesh <- inla.mesh.2d(cbind(0, 0), offset = 2, max.edge = 10)
loc <- matrix(runif(10) * 2 - 1, 5, 2)
ggplot() +
  gg(mesh) +
  geom_point(aes(loc[, 1], loc[, 2])) +
  coord_equal()

A.loc <- inla.spde.make.A(mesh, loc = loc)
A.loc
#> 5 x 9 sparse Matrix of class "dgCMatrix"
#>                                                                          
#> [1,] . 0.29227293 0.09984048 .         .          .         . . 0.6078866
#> [2,] . 0.02324195 0.36614606 .         .          .         . . 0.6106120
#> [3,] . .          0.11935641 0.1416259 .          .         . . 0.7390177
#> [4,] . .          .          0.3522410 0.08180227 .         . . 0.5659568
#> [5,] . .          .          .         0.19233634 0.2974006 . . 0.5102631

With a bru_mapper, this becomes

mapper <- bru_mapper(mesh)
A.loc <- ibm_jacobian(mapper, input = loc)
A.loc
#> 5 x 9 sparse Matrix of class "dgCMatrix"
#>                                                                          
#> [1,] . 0.29227293 0.09984048 .         .          .         . . 0.6078866
#> [2,] . 0.02324195 0.36614606 .         .          .         . . 0.6106120
#> [3,] . .          0.11935641 0.1416259 .          .         . . 0.7390177
#> [4,] . .          .          0.3522410 0.08180227 .         . . 0.5659568
#> [5,] . .          .          .         0.19233634 0.2974006 . . 0.5102631

Mapping with a precomputed location mapping

index <- c(1, 3, 5, 2, 1, 2)
inla.spde.make.A(A.loc = A.loc, index = index)
#> 6 x 9 sparse Matrix of class "dgCMatrix"
#>                                                                         
#> [1,] . 0.29227293 0.09984048 .         .         .         . . 0.6078866
#> [2,] . .          0.11935641 0.1416259 .         .         . . 0.7390177
#> [3,] . .          .          .         0.1923363 0.2974006 . . 0.5102631
#> [4,] . 0.02324195 0.36614606 .         .         .         . . 0.6106120
#> [5,] . 0.29227293 0.09984048 .         .         .         . . 0.6078866
#> [6,] . 0.02324195 0.36614606 .         .         .         . . 0.6106120

mapper <- bru_mapper_taylor(jacobian = A.loc[index, , drop = FALSE])
ibm_jacobian(mapper)
#> 6 x 9 sparse Matrix of class "dgCMatrix"
#>                                                                         
#> [1,] . 0.29227293 0.09984048 .         .         .         . . 0.6078866
#> [2,] . .          0.11935641 0.1416259 .         .         . . 0.7390177
#> [3,] . .          .          .         0.1923363 0.2974006 . . 0.5102631
#> [4,] . 0.02324195 0.36614606 .         .         .         . . 0.6106120
#> [5,] . 0.29227293 0.09984048 .         .         .         . . 0.6078866
#> [6,] . 0.02324195 0.36614606 .         .         .         . . 0.6106120

# For run-time indexing:
mapper <-
  bru_mapper_pipe(
    list(
      matrix = bru_mapper_taylor(jacobian = A.loc),
      index = bru_mapper_index(nrow(A.loc))
    )
  )
ibm_jacobian(mapper, input = list(index = index))
#> 6 x 9 sparse Matrix of class "dgCMatrix"
#>                                                                         
#> [1,] . 0.29227293 0.09984048 .         .         .         . . 0.6078866
#> [2,] . .          0.11935641 0.1416259 .         .         . . 0.7390177
#> [3,] . .          .          .         0.1923363 0.2974006 . . 0.5102631
#> [4,] . 0.02324195 0.36614606 .         .         .         . . 0.6106120
#> [5,] . 0.29227293 0.09984048 .         .         .         . . 0.6078866
#> [6,] . 0.02324195 0.36614606 .         .         .         . . 0.6106120

Group mapping with a group mesh

inla.spde.make.A(..., group = group.values, group.mesh = group.mesh)

mapper <- bru_mapper_multi(list(
  main = bru_mapper(mesh),
  group = bru_mapper(group.mesh)
))
ibm_jacobian(mapper, input = list(main = loc, group = group.values))

Blockwise aggregation

block_rescale <- # one of "none", "count", "weights", "sum"
  inla.spde.make.A(...,
    weights = weights,
    block = block,
    block.rescale = block_rescale,
    n.block = n_block
  )

mapper <- bru_mapper_pipe(
  list(
    main = bru_mapper_multi(list(main = bru_mapper(mesh), ...)),
    block = bru_mapper_aggregate(
      rescale = (block_rescale != "none"),
      n_block = n_block
    )
  )
)
ibm_jacobian(mapper,
  input = list(
    main = list(main = loc),
    block = list(block = block, weights = weights)
  )
)

inla.spde.make.index

ngroup <- 2
nrepl <- 3

summary(
  as.data.frame(
    inla.spde.make.index("field",
      n.spde = mesh$n,
      n.group = ngroup,
      n.repl = nrepl
    )
  )
)
#>      field    field.group    field.repl
#>  Min.   :1   Min.   :1.0   Min.   :1   
#>  1st Qu.:3   1st Qu.:1.0   1st Qu.:1   
#>  Median :5   Median :1.5   Median :2   
#>  Mean   :5   Mean   :1.5   Mean   :2   
#>  3rd Qu.:7   3rd Qu.:2.0   3rd Qu.:3   
#>  Max.   :9   Max.   :2.0   Max.   :3

mapper <- bru_mapper_multi(list(
  field.main = bru_mapper(mesh),
  field.group = bru_mapper_index(ngroup),
  field.replicate = bru_mapper_index(nrepl)
))
summary(ibm_values(mapper, multi = TRUE, inla_f = TRUE))
#>    field.main  field.group  field.replicate
#>  Min.   :1    Min.   :1.0   Min.   :1      
#>  1st Qu.:3    1st Qu.:1.0   1st Qu.:1      
#>  Median :5    Median :1.5   Median :2      
#>  Mean   :5    Mean   :1.5   Mean   :2      
#>  3rd Qu.:7    3rd Qu.:2.0   3rd Qu.:3      
#>  Max.   :9    Max.   :2.0   Max.   :3

The benefit of the mapper approach here is that it encapsulates all the information, so that only the mapper needs to be carried around to code that needs it, and that it doesn’t restrict the group and replicate mappings to integer indices; the index mappers can be replaced by other mappers, e.g. to allow interpolation between group indices, with a 1d mesh mapper.