#' Create an inlabru model object from model components
#'
#' The [inlabru] syntax for model formulae is different from what
#' `INLA::inla` considers a valid.
#' In inla most of the effects are defined by adding an `f(...)` expression to
#' the formula.
#' In [inlabru] the `f` is replaced by an arbitrary (exceptions: `const` and
#' `offset`) string that will determine the label of the effect. See Details for
#' further information.
#'
#' @details
#' For instance
#'
#' `y ~ f(myspde, ...)`
#'
#' in INLA is equivalent to
#'
#' `y ~ myspde(...)`
#'
#' in inlabru.
#'
#' A disadvantage of the inla way is that there is no clear separation between
#' the name of the covariate and the label of the effect. Furthermore, for some
#' models like SPDE it is much more natural to use spatial coordinates as
#' covariates rather than an index into the SPDE vertices. For this purpose
#' [inlabru] provides the `main` argument. For convenience, the `main` argument
#' can be used like the first argument of the f function, e.g., and is the first
#' argument of the component definition.
#' The `INLA` model formula
#'
#' `y ~ f(temperature, model = 'linear')`
#'
#' is equivalent to the `inlabru` component and formula definition
#'
#' `y ~ temperature(temperature, model = 'linear')`
#' and
#' `y ~ temperature(main = temperature, model = 'linear')`
#' as well as
#' `y ~ temperature(model = 'linear')`
#' which sets `main = temperature`.
#'
#' On the other hand, `main` can also be a function mapping, e.g the
#' [sp::coordinates()] function:
#'
#' `y ~ mySPDE(coordinates, ...)`
#'
#' This extracts the coordinates from the data object, and maps it to the latent
#' field via the information given in the `mapper`, which by default is
#' extracted from the `model` object, in the case of `spde` model objects.
#'
#' Morevover, `main` can be any expression that evaluates within your data as an
#' environment.
#' For instance, if your data has columns 'a' and 'b', you can create a fixed
#' effect of 'sin(a+b)' by setting `map` in the following way:
#'
#' `y ~ myEffect(sin(a+b))`
#'
#'
#' @export
#' @param components A [component_list] object
#' @param lhoods A list of one or more `lhood` objects
#' @return A [bru_model] object
#' @keywords internal
bru_model <- function(components, lhoods) {
stopifnot(inherits(components, "component_list"))
# Back up environment
env <- environment(components)
# Create joint formula that will be used by inla
formula <- BRU_response ~ -1
included <- bru_used(lhoods)
included <- union(included[["effect"]], included[["latent"]])
# Detect linearity
linear <- all(vapply(lhoods, function(lh) lh[["linear"]], TRUE))
# TODO: detect pure effect additivity (allowing nonlinear components)
# Complete the used component definitions based on data
components <- component_list(components[included], lhoods)
for (cmp in included) {
if (linear ||
!(components[[cmp]][["main"]][["type"]] %in% c("offset", "const"))) {
formula <- update.formula(formula, components[[cmp]]$inla.formula)
}
}
# Restore environment
environment(components) <- env
environment(formula) <- env
# Make model
mdl <- structure(
list(effects = components, formula = formula),
class = "bru_model"
)
return(mdl)
}
#' @export
#' @method summary bru_model
#' @param object Object to operate on
#' @param \dots Arguments passed on to other methods
#' @rdname bru_model
summary.bru_model <- function(object, ...) {
result <- structure(
list(
components =
summary(object[["effects"]], ...)
),
class = "summary_bru_model"
)
result
}
#' @export
#' @param x A `summary_bru_model` object to be printed
#' @rdname bru_model
print.summary_bru_model <- function(x, ...) {
print(x[["components"]])
invisible(x)
}
#' Evaluate or sample from a posterior result given a model and locations
#'
#' @export
#' @param model A [bru] model
#' @param state list of state lists, as generated by [evaluate_state()]
#' @param data A `list`, `data.frame`, or `Spatial*DataFrame`, with coordinates
#' and covariates needed to evaluate the predictor.
#' @param input Precomputed inputs list for the components
#' @param comp_simple Precomputed `comp_simple_list` for the components
#' @param predictor A formula or an expression to be evaluated given the
#' posterior or for each sample thereof. The default (`NULL`) returns a
#' `data.frame` containing the sampled effects. In case of a formula the right
#' hand side is used for evaluation.
#' @param format character; determines the storage format of predictor output.
#' Available options:
#' * `"auto"` If the first evaluated result is a vector or single-column matrix,
#' the "matrix" format is used, otherwise "list".
#' * `"matrix"` A matrix where each column contains the evaluated predictor
#' expression for a state.
#' * `"list"` A list where each element contains the evaluated predictor
#' expression for a state.
#' @param n Number of samples to draw.
#' @param seed If seed != 0L, the random seed
#' @param num.threads Specification of desired number of threads for parallel
#' computations. Default NULL, leaves it up to INLA.
#' When seed != 0, overridden to "1:1"
#' @param used A [bru_used()] object, or NULL (default)
#' @param \dots Additional arguments passed on to `inla.posterior.sample`
#' @details * `evaluate_model` is a wrapper to evaluate model state, A-matrices,
#' effects, and predictor, all in one call.
#'
#' @keywords internal
#' @rdname evaluate_model
evaluate_model <- function(model,
state,
data = NULL,
input = NULL,
comp_simple = NULL,
predictor = NULL,
format = NULL,
used = NULL,
...) {
used <- bru_used(used, labels = names(model$effects))
if (is.null(state)) {
stop("Not enough information to evaluate model states.")
}
if (is.null(input)) {
input <- input_eval(
components = model$effects[used$effect],
data = data,
inla_f = TRUE
)
}
if (is.null(comp_simple) && !is.null(input)) {
comp_simple <- evaluate_comp_simple(model$effects[used$effect],
input = input,
inla_f = TRUE
)
}
if (is.null(comp_simple)) {
effects <- NULL
} else {
effects <- evaluate_effect_multi_state(
comp_simple,
state = state,
input = input
)
}
if (is.null(predictor)) {
return(effects)
}
values <- evaluate_predictor(
model,
state = state,
data = data,
effects = effects,
predictor = predictor,
format = format,
used = used
)
values
}
#' @details * `evaluate_state` evaluates model state properties or samples
#' @param result A `bru` object from [bru()] or [lgcp()]
#' @param property Property of the model components to obtain value from.
#' Default: "mode". Other options are "mean", "0.025quant", "0.975quant",
#' "sd" and "sample". In case of "sample" you will obtain samples from the
#' posterior (see `n` parameter). If `result` is `NULL`, all-zero vectors are
#' returned for each component.
#' @param internal_hyperpar logical; If `TRUE`, return hyperparameter properties
#' on the internal scale. Currently ignored when `property="sample"`.
#' Default is `FALSE`.
#' @export
#' @rdname evaluate_model
#' @keywords internal
evaluate_state <- function(model,
result,
property = "mode",
n = 1,
seed = 0L,
num.threads = NULL,
internal_hyperpar = FALSE,
...) {
# Evaluate random states, or a single property
if (property == "sample") {
state <- post.sample.structured(result,
n = n, seed = seed,
num.threads = num.threads,
...
)
} else if (is.null(result)) {
state <- list(lapply(
model[["effects"]],
function(x) {
rep(0.0, ibm_n(x[["mapper"]]))
}
))
} else {
state <- list(extract_property(
result = result,
property = property,
internal_hyperpar = internal_hyperpar
))
}
state
}
#' @export
#' @rdname evaluate_effect
evaluate_effect_single_state <- function(...) {
UseMethod("evaluate_effect_single_state")
}
#' @export
#' @rdname evaluate_effect
evaluate_effect_multi_state <- function(...) {
UseMethod("evaluate_effect_multi_state")
}
#' Evaluate a component effect
#'
#' Calculate latent component effects given some data and the state of the
#' component's internal random variables.
#'
#' @export
#' @keywords internal
#' @param component A `bru_mapper`, `bru_component`, `comp_simple`, or
#' `comp_simple_list`.
#' @param input Pre-evaluated component input
#' @param state Specification of one (for `evaluate_effect_single_state`) or
#' several (for `evaluate_effect_multi_State`) latent variable states:
#' * `evaluate_effect_single_state.bru_mapper`:
#' A vector of the latent component state.
#' * `evaluate_effect_single_state.*_list`: list of named state vectors.
#' * `evaluate_effect_multi_state.*_list`: list of lists of named state vectors.
#' @param ... Optional additional parameters, e.g. `inla_f`. Normally unused.
#' @param label Option label used for any warning messages, specifying the
#' affected component.
#' @author Fabian E. Bachl \email{bachlfab@@gmail.com} and
#' Finn Lindgren \email{finn.lindgren@@gmail.com}
#' @rdname evaluate_effect
#' @keywords internal
evaluate_effect_single_state.bru_mapper <- function(component, input, state,
...,
label = NULL) {
values <- ibm_eval(component, input = input, state = state, ...)
not_ok <- ibm_invalid_output(
component,
input = input,
state = state
)
if (any(not_ok)) {
if (is.null(label)) {
warning("Inputs for a mapper give some invalid outputs.",
immediate. = TRUE
)
} else {
warning("Inputs for '", label, "' give some invalid outputs.",
immediate. = TRUE
)
}
}
as.vector(as.matrix(values))
}
#' @return * `evaluate_effect_single_state.component_list`: A list of evaluated
#' component effect values
#' @export
#' @rdname evaluate_effect
#' @keywords internal
evaluate_effect_single_state.comp_simple_list <- function(components,
input,
state,
...) {
result <- list()
for (label in names(components)) {
result[[label]] <- evaluate_effect_single_state(
components[[label]],
input = input[[label]],
state = state[[label]],
...,
label = label
)
}
result
}
#' @return * `evaluate_effect_multi.comp_simple_list`: A list of lists of
#' evaluated component effects, one list for each state
#' @export
#' @rdname evaluate_effect
#' @keywords internal
evaluate_effect_multi_state.comp_simple_list <- function(components,
input,
state,
...) {
lapply(
state,
function(x) {
evaluate_effect_single_state(
components,
input = input,
state = x,
...
)
}
)
}
#' @export
#' @rdname evaluate_effect
#' @keywords internal
evaluate_effect_single_state.component_list <- function(components,
input,
state,
...) {
comp_simple <- evaluate_comp_simple(components, input = input, ...)
evaluate_effect_single_state(comp_simple, input = input, state = state, ...)
}
#' @export
#' @rdname evaluate_effect
#' @keywords internal
evaluate_effect_multi_state.component_list <- function(components,
input,
state,
...) {
comp_simple <- evaluate_comp_simple(components, input = input, ...)
evaluate_effect_multi_state(comp_simple, input = input, state = state, ...)
}
#' Evaluate component effects or expressions
#'
#' Evaluate component effects or expressions, based on a bru model and one or
#' several states of the latent variables and hyperparameters.
#'
#' @param data A `list`, `data.frame`, or `Spatial*DataFrame`, with coordinates
#' and covariates needed to evaluate the model.
#' @param state A list where each element is a list of named latent state
#' information, as produced by [evaluate_state()]
#' @param effects A list where each element is list of named evaluated effects,
#' as computed by [evaluate_effect_multi_state.component_list()]
#' @param predictor Either a formula or expression
#' @param used A [bru_used()] object, or NULL (default)
#' @param format character; determines the storage format of the output.
#' Available options:
#' * `"auto"` If the first evaluated result is a vector or single-column matrix,
#' the "matrix" format is used, otherwise "list".
#' * `"matrix"` A matrix where each column contains the evaluated predictor
#' expression for a state.
#' * `"list"` A list where each column contains the evaluated predictor
#' expression for a state.
#'
#' Default: "auto"
#' @details For each component, e.g. "name", the state values are available as
#' `name_latent`, and arbitrary evaluation can be done with `name_eval(...)`,
#' see [component_eval()].
#' @return A list or matrix is returned, as specified by `format`
#' @keywords internal
#' @rdname evaluate_predictor
evaluate_predictor <- function(model,
state,
data,
effects,
predictor,
used = NULL,
format = "auto") {
stopifnot(inherits(model, "bru_model"))
format <- match.arg(format, c("auto", "matrix", "list"))
pred.envir <- environment(predictor)
if (inherits(predictor, "formula")) {
predictor <- parse(
text = as.character(predictor)[length(as.character(predictor))]
)
}
formula.envir <- environment(model$formula)
enclos <-
if (!is.null(pred.envir)) {
pred.envir
} else if (!is.null(formula.envir)) {
formula.envir
} else {
parent.frame()
}
used <- bru_used(used, labels = names(model$effects))
envir <- new.env(parent = enclos)
# Find .data. first,
# then data variables,
# then pred.envir variables (via enclos),
# then formula.envir (via enclos if pred.envir is NULL):
# for (nm in names(pred.envir)) {
# assign(nm, pred.envir[[nm]], envir = envir)
# }
#
# Note: Since 2.7.0.9019, no longer converts Spatial*DataFrame to data frame
# here; coordinates must be accessed via sp::coordinates() if needed.
for (nm in names(data)) {
assign(nm, data[[nm]], envir = envir)
}
assign(".data.", data, envir = envir)
# Rename component states from label to label_latent
state_names <- as.list(expand_labels(
names(state[[1]]),
names(model$effects),
suffix = "_latent"
))
names(state_names) <- names(state[[1]])
# Construct _eval function names
eval_names <- as.list(expand_labels(
intersect(names(state[[1]]), names(model$effects)),
intersect(names(state[[1]]), names(model$effects)),
suffix = "_eval"
))
names(eval_names) <- intersect(names(state[[1]]), names(model$effects))
eval_fun_factory <-
function(.comp, .envir, .enclos) {
.is_offset <- .comp$main$type %in% c("offset", "const")
.is_iid <- .comp$main$type %in% c("iid")
.mapper <- .comp$mapper
.label <- paste0(.comp$label, "_latent")
.iid_precision <- paste0("Precision_for_", .comp$label)
.iid_cache <- list()
.iid_cache_index <- NULL
eval_fun <- function(main, group = NULL,
replicate = NULL,
weights = NULL,
.state = NULL) {
if (is.null(group)) {
group <- rep(1, NROW(main))
}
if (is.null(replicate)) {
replicate <- rep(1, NROW(main))
}
if (!.is_offset && is.null(.state)) {
.state <- eval(
parse(text = .label),
envir = .envir,
enclos = .enclos
)
}
.input <- list(
mapper = list(
main = main,
group = group,
replicate = replicate
),
scale = weights
)
.values <- ibm_eval(
.mapper,
input = .input,
state = .state
)
if (!.is_iid) {
not_ok <- ibm_invalid_output(
.mapper[["mappers"]][[1]],
input = .input[[1]],
state = .state
)
if (any(not_ok)) {
warning("Inputs for `ibm_eval()` for '",
.comp[["label"]],
'" give some invalid outputs.')
}
} else { # .is_iid, invalid indices give new samples
# Check for known invalid output elements, based on the
# initial mapper (subsequent mappers in the component pipe
# are assumed to keep the same length and validity)
not_ok <- ibm_invalid_output(
.mapper[["mappers"]][[1]],
input = .input[[1]],
state = .state
)
if (any(not_ok)) {
.cache_state_index <- eval(
parse(text = ".cache_state_index"),
envir = .envir,
enclos = .enclos
)
if (!identical(.cache_state_index, .iid_cache_index)) {
.iid_cache_index <<- .cache_state_index
.iid_cache <<- list()
}
key <- as.character(main[not_ok])
not_cached <- !(key %in% names(.iid_cache))
if (any(not_cached)) {
.prec <- eval(
parse(text = .iid_precision),
envir = .envir,
enclos = .enclos
)
for (k in unique(key[not_cached])) {
.iid_cache[k] <<- rnorm(1, mean = 0, sd = .prec^-0.5)
}
}
.values[not_ok] <- vapply(
key,
function(k) .iid_cache[[k]],
0.0
)
}
}
as.matrix(.values)
}
eval_fun
}
for (nm in names(eval_names)) {
assign(eval_names[[nm]],
eval_fun_factory(model$effects[[nm]], .envir = envir, .enclos = enclos),
envir = envir
)
}
# Remove problematic objects:
problems <- c(".Random.seed")
remove(list = intersect(names(envir), problems), envir = envir)
n <- length(state)
for (k in seq_len(n)) {
# Keep track of the iteration index so the iid cache can be invalidated
assign(".cache_state_index", k, envir = envir)
for (nm in names(state[[k]])) {
assign(state_names[[nm]], state[[k]][[nm]], envir = envir)
}
for (nm in names(effects[[k]])) {
assign(nm, effects[[k]][[nm]], envir = envir)
}
result_ <- eval(predictor, envir = envir, enclos = enclos)
if (k == 1) {
if (identical(format, "auto")) {
if ((is.vector(result_) && !is.list(result_)) ||
(is.matrix(result_) && (NCOL(result_) == 1))) {
format <- "matrix"
} else {
format <- "list"
}
}
if (identical(format, "matrix")) {
result <- matrix(0.0, NROW(result_), n)
rownames(result) <- row.names(as.matrix(result_))
} else if (identical(format, "list")) {
result <- vector("list", n)
}
}
if (identical(format, "list")) {
result[[k]] <- result_
} else {
result[, k] <- result_
}
}
result
}
#' Evaluate component values in predictor expressions
#'
#' In predictor expressions, `name_eval(...)` can be used to evaluate
#' the effect of a component called "name".
#'
#' @param main,group,replicate,weights Specification of where to evaluate a
#' component. The four inputs are passed on to the joint `bru_mapper` for the
#' component, as
#' ```
#' list(mapper = list(
#' main = main,
#' group = group,
#' replicate = replicate),
#' scale = weights)
#' ```
#' @param .state The internal component state. Normally supplied automatically
#' by the internal methods for evaluating inlabru predictor expressions.
#' @return A vector of values for a component
#' @examples
#' \dontrun{
#' if (bru_safe_inla()) {
#' mesh <- fmesher::fm_mesh_2d_inla(
#' cbind(0, 0),
#' offset = 2, max.edge = 0.25
#' )
#' spde <- INLA::inla.spde2.pcmatern(mesh,
#' prior.range = c(0.1, 0.01),
#' prior.sigma = c(2, 0.01)
#' )
#' data <- sp::SpatialPointsDataFrame(
#' matrix(runif(10), 5, 2),
#' data = data.frame(z = rnorm(5))
#' )
#' fit <- bru(z ~ -1 + field(coordinates, model = spde),
#' family = "gaussian", data = data
#' )
#' pred <- predict(
#' fit,
#' data = data.frame(x = 0.5, y = 0.5),
#' formula = ~ field_eval(cbind(x, y))
#' )
#' }
#' }
component_eval <- function(main,
group = NULL,
replicate = NULL,
weights = NULL,
.state = NULL) {
stop(paste0(
"In your predictor expression, use 'mylabel_eval(...)' instead of\n",
"'component_eval(...)'. See ?component_eval for more information."
))
}
#' Compute all component linearisations
#'
#' Computes individual `bru_mapper_taylor` objects for included components
#' for each model likelihood
#'
#' @param model A `bru_model` object
#' @param input A list of named lists of component inputs
#' @param state A named list of component states
#' @param inla_f Controls the input data interpretations
#' @return A list (class 'comp_simple') of named lists (class
#' 'comp_simple_list') of `bru_mapper_taylor` objects, one for each included
#' component
#' @rdname evaluate_comp_lin
#' @keywords internal
evaluate_comp_lin <- function(model, input, state, inla_f = FALSE) {
stopifnot(inherits(model, "bru_model"))
mappers <-
lapply(
input,
function(inp) {
included <- parse_inclusion(
names(model[["effects"]]),
names(inp),
NULL
)
mappers <- comp_lin_eval(
model[["effects"]][included],
input = inp[included],
state = state[included],
inla_f = inla_f
)
class(mappers) <- c("comp_simple_list", class(mappers))
mappers
}
)
class(mappers) <- c("comp_simple_list_list", class(mappers))
mappers
}
#' Compute simplified component mappings
#'
#' Computes individual `bru_mapper_taylor` objects for included linear
#' components for each model likelihood, and keeps non-linear mappers intact.
#'
#' @param model A `bru_model` object
#' @param input A list of named lists of component inputs
#' @param inla_f Controls the input data interpretations
#' @return A list (class 'comp_simple_list_list') of named lists (class
#' 'comp_simple_list') of `bru_mapper` objects, one for each included
#' component
#' @export
#' @keywords internal
#' @rdname evaluate_comp_simple
evaluate_comp_simple <- function(...) {
UseMethod("evaluate_comp_simple")
}
#' @export
#' @rdname evaluate_comp_simple
evaluate_comp_simple.component_list <- function(components, input,
inla_f = FALSE, ...) {
mappers <- lapply(
components,
function(x) {
label <- x[["label"]]
ibm_simplify(
x[["mapper"]],
input[[label]],
state = NULL,
inla_f = inla_f
)
}
)
class(mappers) <- c("comp_simple_list", class(mappers))
mappers
}
#' @export
#' @rdname evaluate_comp_simple
evaluate_comp_simple.bru_model <- function(model, input, ...) {
mappers <-
lapply(
input,
function(inp) {
included <- parse_inclusion(
names(model[["effects"]]),
names(inp),
NULL
)
evaluate_comp_simple(model[["effects"]][included], input = inp, ...)
}
)
class(mappers) <- c("comp_simple_list_list", class(mappers))
mappers
}
#' Subsetting of comp_simple_list objects, retaining class
#' @export
#' @param x `comp_simple_list` object from which to extract element(s)
#' @param i indices specifying elements to extract
#' @keywords internal
#' @rdname evaluate_comp_simple_list_subsetting
`[.comp_simple_list` <- function(x, i) {
env <- environment(x)
object <- NextMethod()
class(object) <- c("comp_simple_list", "list")
environment(object) <- env
object
}
#' Compute all component inputs
#'
#' Computes the component inputs for included components
#' for each model likelihood
#'
#' @param model A `bru_model` object
#' @param lhoods A `bru_like_list` object
#' @param inla_f logical
#' @rdname evaluate_inputs
#' @keywords internal
evaluate_inputs <- function(model, lhoods, inla_f) {
stopifnot(inherits(model, "bru_model"))
lapply(
lhoods,
function(lh) {
included <- bru_used(lh)[["effect"]]
input_eval(
model$effects[included],
data = lh[["data"]],
inla_f = inla_f
)
}
)
}
#' Compute all index values
#'
#' Computes the index values matrices for included components
#'
#' @param model A `bru_model` object
#' @param lhoods A `bru_like_list` object. Deprecated and ignored
#' @return A named list of `idx_full` and `idx_inla`,
#' named list of indices, and `inla_subset`, and `inla_subset`,
#' a named list of logical subset specifications for extracting the `INLA::f()`
#' compatible index subsets.
#' @rdname evaluate_index
#' @keywords internal
evaluate_index <- function(model, lhoods) {
stopifnot(inherits(model, "bru_model"))
included <- bru_used(lhoods)
included <- union(included[["effect"]], included[["latent"]])
list(
idx_full = index_eval(model[["effects"]][included], inla_f = FALSE),
idx_inla = index_eval(model[["effects"]][included], inla_f = TRUE),
inla_subset = inla_subset_eval(model[["effects"]][included])
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.