knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = FALSE )
(Vignette under construction!)
library(inlabru)
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.
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.
Basic mappers take covariate vectors or matrices as input and numeric vectors as state.
Constructors:
bru_mapper_const()
bru_mapper_linear()
bru_mapper_index(n)
bru_mapper_factor(values, factor_mapping)
bru_mapper_matrix()
bru_mapper_harmonics(order, scaling, intercept, interval)
The const
mapper defines a mapping from an empty state vector. This is used to
define offset components, $\eta(u)_j = z_j$, where $z$ is a fixed input vector,
indexed by $j$.
The linear
mapper defines a mapper from a length 1 state vector. This is used to
define ordinary "fixed" covariate effects, linear in both the covariate and in the state
variable; $\eta(u)_j = z_j u$, where $u$ is the state, and $z$ is the
input covariate vector.
The index(n)
mapper defines a direct mapping from a length n
state
variable. This is used for structured and unstructured "random effect" component
effects; $\eta(u)j = u{z_j}$, where $u$ is the state vector, and $z$ is the
input index vector.
The factor(values, factor_mapping)
mapper defines a direct mapping between
a state vector of length equal to the number of factor levels of values
(for factor_mapping = "full"
), or one less (for factor_mapping = "contrast"
).
Each factor level is represented as a dummy 0/1 variable, or equivalently,
the input is used as a label index into the state vector; $\eta(u)j = u{z_j}$.
In the "contrast"
case, the first factor level is removed from the model,
and the corresponding effect defined to be zero.
The matrix
mapper defines a direct matrix multiplication between
a pre-computed model matrix $A_\text{input}$ with the state
vector; $\eta(u) = A_\text{input} u$.
This is used e.g. for linear model formula input, that is converted to a
component model matrix before handing it over to the mapper system.
The harmonics(order, scaling, intercept, interval)
mapper defines a sum of
weighted harmonics on a real interval $(L, U)$, each frequency
additionally scaled by factors determined by scale
;
$\eta(u_j) = \sum_{k=1}^{1+2p} A_{j,k} u_k$, where $A_{j,k}=s_k$ for $k=1$,
$A_{j,2k}=s_{k+1}\cos[2\pi k (z_j - L)/(U-L)]$ for $k=1,2,\dots,p$, and
$A_{j,2k+1}=s_{k+1}\sin[2\pi k (z_j - L)/(U-L)]$ for $k=1,2,\dots,p$,
where $p=$order
. When intercept
is TRUE
, the state vector has length
$1+2p$.
When intercept
is FALSE
, the first, constant, column is
removed, and the state vector has length $2p$.
Transformation mappers are mappers that would normally be combined with other mappers, as steps in a sequence of transformations, or as individual transformation mappers.
bru_mapper_scale()
bru_mapper_marginal(qfun, pfun, ...)
bru_mapper_aggregate(rescale)
bru_mapper_logsumexp(rescale)
The scale
mapper multiplies the input
vector with the state
vector.
This is used for INLA models to implement the weights
component scaling,
as the final stage of a pipe
mapper, see below.
mapper <- bru_mapper_scale() ibm_eval(mapper, input = ..., state = ... )
The marginal
mapper transforms the state
vector from standardised Gaussian marginal distribution
to the distribution defined by a quantile function.
This can be used to define components with latent marginal N(0,1)
distribution but non-Gaussian
effect.
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)
The aggregate
mapper aggregates the state
vector elements by blockwise
summation after scaling the elements by weights.
This can be used for summation, integration, and averaging (for rescale=TRUE
).
mapper <- bru_mapper_aggregate(rescale = ...) ibm_eval(mapper, input = list(block = ..., weights = ...), state = ... )
The logsumexp
mapper aggregates the exp(state)
vector elements by blockwise
summation after scaling the elements by weights, and takes the logarithm. The
implementation takes care to avoid numerical overflow.
This can be used for summation, integration, and averaging (for rescale=TRUE
).
mapper <- bru_mapper_logsumexp(rescale = ...) ibm_eval(mapper, input = list(block = ..., weights = ...), state = ... )
Compound mappers define collections or chains of mappings, and can take various forms of input. The state vector is normally a numeric vector, but can in some cases be a list of vectors.
bru_mapper_collect(mappers, hidden)
bru_mapper_multi(mappers)
bru_mapper_pipe(mappers)
The collect
mapper defines a compound mapper where the Jacobian is
block-diagonal, and each block is defined by a separate sub-mapper.
When hidden = TRUE
, it can be used for INLA models where the user-visible
part of the latent state is only a part of the internal state, such as for
"bym"
models.
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 )
The multi
mapper defines a compound mapper where the Jacobian is
the row-wise Kronecker product between the sub-mapper Jacobians.
This can be used for INLA models to construct the internal mapper
for main
, group
, and replicate
, but also for user-defined models
where a single rgeneric
or cgeneric
model is defined on e.g.
a space-time product domain.
mapper <- bru_mapper_multi(list(name1 = ..., name2 = ..., ...)) ibm_eval(mapper, input = list(name1 = ..., name2 = ..., ...), state = ... )
Pipe mappers chain multiple mappers into a sequence, where the evaluation result of each mapper is given as the state vector for the next mapper.
mapper <- bru_mapper_pipe(list(name1 = ..., name2 = ..., ...)) ibm_eval(mapper, input = list(name1 = ..., name2 = ..., ...), state = ... )
All model component mappers are currently defined as a pipe
mapper containing
a multi
mapper followed by a scale
mapper:
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 = ... )
For object with predefined mapper conversion methods, use bru_mapper(object)
to obtain a suitable mapper object.
Current predefined object mappers are defined for objects of these classes:
fm_mesh_2d
and inla.mesh
; Use bru_mapper(object)
fm_mesh_1d
and inla.mesh.1d
; Use bru_mapper(object, indexed)
More details are given below.
For inla model objects with predefined mappers, use bru_get_mapper(object)
Current predefined model object mappers are defined for models of these classes:
inla.spde
; this includes inla.spde2.matern()
and inla.spde2.pcmatern()
models. The conversion method calls bru_mapper()
in the appropriate way for
the mesh type used for the modelinla.rgeneric
; This relies on the rgeneric definition including a "mapper"
callback argument. A less invasive method, that works for both rgeneric and cgeneric
models, is to define a bru_get_mapper()
method for the class, which needs to
include a unique class identifier, e.g. add
class(object) <- c("my_unique_model_class", class(object))
and definebru_get_mapper.my_unique_model_class <- function(model, ...) { ... }
and register it in the namespace.
bru_mapper_taylor()
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) ibm_values(mapper) ibm_n(mapper, inla_f = TRUE) ibm_values(mapper, inla_f = TRUE)
ibm_n(mapper, multi = TRUE) ibm_values(mapper, multi = TRUE) ibm_n(mapper, inla_f = TRUE, multi = TRUE) ibm_values(mapper, inla_f = TRUE, multi = TRUE)
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) ibm_n(mapper, multi = TRUE) ibm_values(mapper) ibm_values(mapper, multi = TRUE)
ibm_n(mapper, inla_f = TRUE) ibm_n(mapper, multi = TRUE, inla_f = TRUE) ibm_values(mapper, inla_f = TRUE) ibm_values(mapper, multi = TRUE, inla_f = TRUE)
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.
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)
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.
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 $\bm{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.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.