Devel: Customised model components with the bru_mapper system

Mapper system introduction

Each inlabru latent model component generates an effect, given the latent state vector. The purpose of the mapper system is to define the link from state vectors to effect vectors. For an ordinary model component, named \(c\), this link can be represented as a matrix-vector product \[ \eta_c(u_c) = A_c(\text{input}) u_c, \] where \(u_c\) is the latent state vector, \(\eta_c(u_c)\) is the resulting component effect vector, and \(A_c(\text{input})\) is a component model matrix. This matrix depends on the component inputs (covariates, index information, etc) from main, group, replicate, and weights in the component definition.

The wider scope of component definitions is discussed in the component vignette. Here, we will focus on the mapper system itself, how the basic building blocks are used to construct the built-in mappers, and finally how to construct new mappers.

The in-package documentation for the mapper methods is contained in four parts:

?bru_mapper # Mapper constructors
?bru_mapper_generics # Generic and default methods
?bru_mapper_methods # Specialised mapper methods
?bru_get_mapper # Mapper extraction methods

Regular users normally at most need some of the methods from ?bru_mapper and sometimes ?bru_mapper_generics. The methods in ?bru_mapper_methods provides more details and allows the user to query mappers and to use them outside the context of inlabru model definitions. The bru_mapper_define() and bru_get_mapper() methods are needed for those implementing their own mapper class.

Mappers

The main purpose of each mapper class is to allow evaluating a component effect, from given input and latent state, by calling

ibm_eval(mapper, input, state)

for the mapper associated with the component definition.

mapper <- bru_mapper_scale()
ibm_eval(mapper,
  input = ...,
  state = ...
)

mapper <- bru_mapper_marginal(qfun = ..., pfun = ..., dfun = ..., ..., inverse = ...)
ibm_eval(mapper,
  input = NULL, # If a list is given, it overrides the parameter specification
  state = ...
)

# Examples:
mapper <- bru_mapper_marginal(qfun = qexp, rate = 1 / 8)
mapper <- bru_mapper_marginal(qfun = qexp, pfun = pexp, dfun = dexp, rate = 1 / 8)

mapper <- bru_mapper_aggregate(rescale = ...)
ibm_eval(mapper,
  input = list(block = ..., weights = ...),
  state = ...
)

mapper <- bru_mapper_logsumexp(rescale = ...)
ibm_eval(mapper,
  input = list(block = ..., weights = ...),
  state = ...
)

mapper <- bru_mapper_collect(list(name1 = ..., name2 = ..., ...),
  hidden = FALSE
)
ibm_eval(mapper,
  input = list(name1 = ..., name2 = ..., ...),
  state = ...
)
# If hidden = TRUE, the inla_f=TRUE argument "hides" all but the first mapper:
ibm_eval(mapper,
  input = name1_input,
  state = ...,
  inla_f = TRUE
)

mapper <- bru_mapper_multi(list(name1 = ..., name2 = ..., ...))
ibm_eval(mapper,
  input = list(name1 = ..., name2 = ..., ...),
  state = ...
)

mapper <- bru_mapper_pipe(list(name1 = ..., name2 = ..., ...))
ibm_eval(mapper,
  input = list(name1 = ..., name2 = ..., ...),
  state = ...
)

mapper <-
  bru_mapper_pipe(
    list(
      mapper = bru_mapper_multi(list(main = ..., group = ..., replicate = ...)),
      scale = bru_mapper_scale()
    )
  )
ibm_eval(mapper,
  input = list(
    mapper = list(main = ..., group = ..., replicate = ...),
    scale = ...
  ),
  state = ...
)

bru_get_mapper.my_unique_model_class <- function(model, ...) {
  ...
}

Mapper methods

ibm_n(mapper, inla_f, ...)
ibm_n_output(mapper, input, ...)
ibm_values(mapper, inla_f, ...)
ibm_jacobian(mapper, input, state, ...)
ibm_eval(mapper, input, state, ...)
ibm_names(mapper, ...)
ibm_inla_subset(mapper, ...)

An example where the inla_f argument matters is the bru_mapper_collect class, when the hidden=TRUE argument is used to indicate that only the first mapper should be used for the INLA::f() inputs, e.g. for "bym2" models. For inla_f=FALSE (the default), the ibm_n and ibm_values methods return the total number of latent variables, but with inla_f=TRUE we get the values needed by INLA::f() instead:

mapper <- bru_mapper_collect(
  list(
    a = bru_mapper_index(3),
    b = bru_mapper_index(2)
  ),
  hidden = TRUE
)
ibm_n(mapper)
#> [1] 5

ibm_values(mapper)
#> [1] 1 2 3 4 5

ibm_n(mapper, inla_f = TRUE)
#> [1] 3

ibm_values(mapper, inla_f = TRUE)
#> [1] 1 2 3

ibm_n(mapper, multi = TRUE)
#> $a
#> [1] 3
#> 
#> $b
#> [1] 2

ibm_values(mapper, multi = TRUE)
#> $a
#> [1] 1 2 3
#> 
#> $b
#> [1] 1 2

ibm_n(mapper, inla_f = TRUE, multi = TRUE)
#> $a
#> [1] 3
#> 
#> $b
#> [1] 2

ibm_values(mapper, inla_f = TRUE, multi = TRUE)
#> $a
#> [1] 1 2 3
#> 
#> $b
#> [1] 1 2

For the bru_mapper_multi class, the multi argument provides access to the inner layers of the multi-mapper:

mapper <- bru_mapper_multi(list(
  a = bru_mapper_index(3),
  b = bru_mapper_index(2)
))
ibm_n(mapper)
#> [1] 6

ibm_n(mapper, multi = TRUE)
#> $a
#> [1] 3
#> 
#> $b
#> [1] 2

ibm_values(mapper)
#> [1] 1 2 3 4 5 6

ibm_values(mapper, multi = TRUE)
#>   a b
#> 1 1 1
#> 2 2 1
#> 3 3 1
#> 4 1 2
#> 5 2 2
#> 6 3 2

ibm_n(mapper, inla_f = TRUE)
#> [1] 6

ibm_n(mapper, multi = TRUE, inla_f = TRUE)
#> $a
#> [1] 3
#> 
#> $b
#> [1] 2

ibm_values(mapper, inla_f = TRUE)
#> [1] 1 2 3 4 5 6

ibm_values(mapper, multi = TRUE, inla_f = TRUE)
#>   a b
#> 1 1 1
#> 2 2 1
#> 3 3 1
#> 4 1 2
#> 5 2 2
#> 6 3 2

The default ibm_inla_subset method determines the inla subset by comparing the full values information from ibm_values(mapper, inla_f = FALSE) to the inla specific values information from ibm_values(mapper, inla_f = TRUE), to determine the logical vector identifying the inla values subset. Custom mappers should normally not need to specialise this method.

Mappers for inla.mesh objects

For component models referenced by a character label (e.g. "iid", "bym2", "rw2" etc), inlabru will construct default mappers, that in most cases replicate the default INLA::f() behaviour.

For the fm_mesh_1d, fm_mesh_2d, inla.mesh and inla.mesh.1d classes, default mappers can be constructed by the pre-defined bru_mapper S3 methods. For 2d meshes (fm_mesh_2d, inla.mesh),

bru_mapper(mesh)

For 1d meshes (fm_mesh_1d, inla.mesh.1d),

# If ibm_values() should return mesh$loc (e.g. for "rw2" models
# with degree=1 meshes)
bru_mapper(mesh, indexed = FALSE)
# If ibm_values() should return seq_along(mesh$loc) (e.g. for
# inla.spde2.pcmatern() models)
bru_mapper(mesh, indexed = TRUE)

Customised mappers

A mapper object should store enough information in order for the ibm_* methods to work. The simplest case of a customised mapper is to just attached a new class label to the front of the S3 class() information of an existing mapper, to obtain a class to override some of the standard ibm_* method implementations. More commonly, one would instead start with a basic list(), that might contain an existing mapper object, and then add methods that know how to use that information. In all these cases, the bru_mapper_define() method should be used to properly set the class information:

bru_mapper_define(mapper, new_class)

Implementations can avoid having to define ibm_n and ibm_values methods, by instead computing and storing n, n_inla, values, and values_inla in the mapper object during construction:

bru_mapper_define(
  mapper = list(n = 10, values = 1:10),
  new_class = "my_bru_mapper_class_name"
)

The default ibm_n and ibm_values methods checks if these values are available, and return the appropriate values depending on the inla_f argument. When *_inla values are requested but not available, the methods fall back to the non-inla versions. If the needed information is not found, the default methods give an error message.

Note: Before version 2.6.0, ibm_amatrix was used instead of ibm_jacobian. Version 2.6.0 still supports this, by having the default ibm_jacobian method call ibm_amatrix, but will give a deprecation message later, and it may be removed in a future version.

Example

Let’s build a mapper class bru_mapper_p_quadratic for a model component that takes input covariates with values \(a_{ij}\), \(i=1,\dots,m\) and \(j=1,\dots,p\), in a matrix or data.frame and evaluates a full quadratic expression \[ \eta_i = x_0 + \sum_{j=1}^p a_{ij} x_j + \frac{1}{2}\sum_{j=1}^p\sum_{k=1}^j \gamma_{j,k} a_{ij}a_{ik} x_{j,k}, \] where the latent component vector is \(\boldsymbol{x}=[x_0,x_1,\dots,x_p,x_{1,1},\dots,x_{p,1},x_{2,2}\dots,x_{p,p}]\), and \[ \gamma_{j,k} = \begin{cases} 1, & j = k,\\ 2, & j > k. \end{cases} \]

We start with the constructor method. Like the bru_mapper_matrix, we require the user to supply a vector of covariate labels, and store them as a character vector. We also include a min_degree parameter to control if an intercept (min_degree <= 0) and linear terms (min_degree <= 1) should be included in the model.

bru_mapper_p_quadratic <- function(labels, min_degree = 0, ...) {
  if (is.factor(labels)) {
    mapper <- list(
      labels = levels(labels),
      min_degree = min_degree
    )
  } else {
    mapper <- list(
      labels = as.character(labels),
      min_degree = min_degree
    )
  }
  bru_mapper_define(mapper, new_class = "bru_mapper_p_quadratic")
}

The ibm_n method can compute the value of \(n\) from \(n=1+p+\frac{p(p+1)}{2}\):

ibm_n.bru_mapper_p_quadratic <- function(mapper, ...) {
  p <- length(mapper$labels)
  (mapper$min_degree <= 0) + (mapper$min_degree <= 1) * p + p * (p + 1) / 2
}

For ibm_values, the default method is sufficient, returning the vector \([1,\dots,n]\) with \(n\) obtained by ibm_n(mapper). However, for clearer result naming, we can use a character vector, at least for INLA f() models that allow it (we could have an option argument to the bru_mapper_p_quadratic constructor to control this):

ibm_values.bru_mapper_p_quadratic <- function(mapper, ...) {
  p <- length(mapper$labels)
  n <- ibm_n(mapper)
  jk <- expand.grid(seq_len(p), seq_len(p))
  jk <- jk[jk[, 2] <= jk[, 1], , drop = FALSE]
  c(
    if (mapper$min_degree <= 0) "Intercept" else NULL,
    if (mapper$min_degree <= 1) mapper$labels else NULL,
    paste0(mapper$labels[jk[, 1]], ":", mapper$labels[jk[, 2]])
  )
}

While the approach to ibm_n and ibm_values above works, it is wasteful to recompute the n and values information for each call. Instead we can use that the default method will check if n or values are stored in the mapper object, and then return them. If we change the . in the method definitions to _, we can end the constructor method like this, making subsequent method calls faster:

bru_mapper_p_quadratic <- function(labels, min_degree = 0, ...) {
  ...
  mapper <- bru_mapper_define(mapper, new_class = "bru_mapper_p_quadratic")
  mapper$n <- ibm_n_bru_mapper_p_quadratic(mapper)
  mapper$values <- ibm_values_bru_mapper_p_quadratic(mapper)
  mapper
}

Note that the order matters, since our ibm_values_bru_mapper_p_quadratic function calls ibm_n(mapper).

We can now define the main part of the mapper interface that computes the model matrix linking the latent variables to the component effect. It’s required that NULL input should return a \(0\)-by-\(n\) matrix.

ibm_jacobian.bru_mapper_p_quadratic <- function(mapper, input, ...) {
  if (is.null(input)) {
    return(Matrix::Matrix(0, 0, ibm_n(mapper)))
  }
  p <- length(mapper$labels)
  n <- ibm_n(mapper)
  N <- NROW(input)
  A <- list()
  in_ <- as(input, "Matrix")
  idx <- 0
  if (mapper$min_degree <= 0) {
    idx <- idx + 1
    A[[idx]] <- Matrix::Matrix(1, N)
  }
  if (mapper$min_degree <= 1) {
    idx <- idx + 1
    A[[idx]] <- in_
  }
  for (k in seq_len(p)) {
    idx <- idx + 1
    A[[idx]] <- in_[, seq(k, p, by = 1), drop = FALSE] * in_[, k]
    A[[idx]][, k] <- A[[idx]][, k] / 2
  }
  A <- do.call(cbind, A)
  colnames(A) <- as.character(ibm_values(mapper))
  A
}

If the output of a mapper is of a different size than NROW(input), the mapper should define a ibm_n_output(mapper, input, ...) method to return the output size for a given input.