R/legacy_methods.R

In inlabru: Bayesian Latent Gaussian Modelling using INLA and Extensions

#' 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 data_extra Additional data for the predictor evaluation
#' @param input Precomputed inputs list for the components
#' @param comp_simple Precomputed [bm_list] of simplified mappers for the
#' components
#' @param predictor A formula or a [bru_pred_expr] 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:1"
#' @param used A [bru_used()] object, or NULL (default)
#' @param n_pred integer. If provided, scalar predictor results are expanded to
#' vectors of length `n_pred`.
#' @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,
  data_extra = NULL,
  input = NULL,
  comp_simple = NULL,
  predictor = NULL,
  format = NULL,
  used = NULL,
  n_pred = NULL,
  ...
) {
  stopifnot(identical(bru_options_get("bru_method")$autodiff, "pandemic"))
  comp_lst <- as_bru_comp_list(model)
  if (inherits(predictor, "bru_pred_expr")) {
    if (!is.null(used)) {
      warning(
        paste0(
          "Overriding `used` argument to `evaluate_model()` ",
          "in favour of `bru_used(predictor)`."
        )
      )
    }
    used <- bru_used(predictor, labels = names(comp_lst))
    predictor <- bru_pred_expr(predictor, format = "formula")
  } else {
    used <- bru_used(used, labels = names(comp_lst))
  }

  if (is.null(state)) {
    stop("Not enough information to evaluate model states.")
  }
  if (is.null(input)) {
    input <- bru_input(
      comp_lst[used$effect],
      data = data
    )
  }
  if (is.null(comp_simple) && !is.null(input)) {
    comp_simple <- ibm_simplify(
      comp_lst[used$effect],
      input = input,
      inla_f = TRUE
    )
  }
  if (is.null(comp_simple)) {
    effects <- NULL
  } else {
    effects <-
      lapply(
        state,
        function(x) {
          evaluate_effect_single_state(
            comp_simple,
            state = x,
            input = input
          )
        }
      )
  }

  if (is.null(predictor)) {
    return(effects)
  }

  values <- evaluate_predictor(
    model,
    state = state,
    data = data,
    data_extra = data_extra,
    effects = effects,
    predictor = predictor,
    format = format,
    used = used,
    n_pred = n_pred
  )

  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,
  ...
) {
  stopifnot(identical(bru_options_get("bru_method")$autodiff, "pandemic"))
  # 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(
      as_bru_comp_list(model),
      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(...) {
  stopifnot(identical(bru_options_get("bru_method")$autodiff, "pandemic"))
  UseMethod("evaluate_effect_single_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_comp], or
#'   [bm_list].
#' @param input Pre-evaluated component input
#' @param state Specification of one latent variable state:
#' * `evaluate_effect_single_state.bru_mapper`:
#'   A vector of the latent component state.
#' * `evaluate_effect_single_state.*_list`: list of named state vectors.
#' @param \dots 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.bm_list`: A list of
#'   evaluated component effect values
#' @export
#' @rdname evaluate_effect
#' @keywords internal
evaluate_effect_single_state.bm_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
}

#' @export
#' @rdname evaluate_effect
#' @keywords internal
evaluate_effect_single_state.bru_comp_list <- function(
  components,
  input,
  state,
  ...
) {
  comp_simple <- ibm_simplify(components, input = input, state = state, ...)
  evaluate_effect_single_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 data_extra Additional data for the predictor evaluation. Variables
#' with the same name as in `data` will be ignored, unless accessed via
#' `.data_extra.[["name"]]` or `.data_extra.$name`, or via pronouns;
#'  see Details.
#' @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,
#' each computed by [evaluate_effect_single_state.bru_comp_list()]
#' @param predictor Either a formula or [bru_pred_expr] 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"
#' @param n_pred integer. If provided, scalar predictor results are expanded to
#' vectors of length `n_pred`.
#' @details For each component, e.g. "name", the latent state values are
#'   available as `name_latent`, and arbitrary evaluation can be done with
#'   `name_eval(...)`, see [bru_comp_eval()].
#'
#'   The evaluation supports several [rlang::as_data_pronoun()] data masking
#'   pronouns, to access variables from different data sources, and some of
#'   these also have corresponding full objects, with an appended `.` in the
#'   name. The full objects can be passed as arguments to functions.
#'   \describe{
#'   \item{.effect/.effect.}{refers to the `effects` vectors}
#'   \item{.latent/.latent.}{refers to the latent state vectors}
#'   \item{.data/.data.}{refers to the main `data` argument}
#'   \item{.data_extra/.data_extra.}{refers to the `data_extra` argument}
#'   \item{.env}{refers to the evaluation environment of the predictor}
#'   }
#' @return A list or matrix is returned, as specified by `format`
#' @keywords internal
#' @rdname evaluate_predictor
evaluate_predictor <- function(
  model,
  state,
  data,
  data_extra,
  effects,
  predictor,
  used = NULL,
  format = "auto",
  n_pred = NULL
) {
  stopifnot(
    identical(bru_options_get("bru_method")$autodiff, "pandemic"),
    inherits(model, "bru_model")
  )
  format <- match.arg(format, c("auto", "matrix", "list"))
  comp_lst <- as_bru_comp_list(model)
  if (inherits(predictor, "bru_pred_expr")) {
    if (!is.null(used)) {
      warning(
        paste0(
          "Overriding `used` argument to `evaluate_predictor()` ",
          "in favour of `bru_used(predictor)`."
        )
      )
    }
    used <- bru_used(predictor, labels = names(comp_lst))
    predictor <- bru_pred_expr(predictor, format = "formula")
  }
  pred.envir <- environment(predictor)
  if (inherits(predictor, "formula")) {
    pred_text <- as.character(predictor)
    pred_text <- pred_text[length(pred_text)]
    predictor <- rlang::parse_expr(pred_text)
  }
  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(comp_lst))

  # General evaluation environment
  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.

  # Rename component states from label to label_latent
  state_names <- as.list(expand_labels(
    names(state[[1]]),
    names(comp_lst),
    suffix = "_latent"
  ))
  names(state_names) <- names(state[[1]])

  # Construct _eval function names
  eval_names <- as.list(expand_labels(
    intersect(names(state[[1]]), names(comp_lst)),
    intersect(names(state[[1]]), names(comp_lst)),
    suffix = "_eval"
  ))
  names(eval_names) <- intersect(names(state[[1]]), names(comp_lst))

  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
      ) {
        n_input <- ibm_n_output(
          .mapper[["mappers"]][["core"]][["mappers"]][["main"]],
          input = main
        )
        if (is.null(group)) {
          group <- rep(1, n_input)
        }
        if (is.null(replicate)) {
          replicate <- rep(1, n_input)
        }
        if (!.is_offset && is.null(.state)) {
          .state <- rlang::eval_tidy(
            rlang::parse_expr(.label),
            data = data_mask,
            env = .envir
          )
        }
        .input <- list(
          core = list(
            main = main,
            group = group,
            replicate = replicate
          ),
          scale = weights
        )

        .values <- ibm_eval(
          .mapper,
          input = .input,
          state = .state
        )
        if (.is_iid) {
          # .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 <- rlang::eval_tidy(
              rlang::parse_expr(".cache_state_index"),
              data = data_mask,
              env = .envir
            )
            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 <- rlang::eval_tidy(
                rlang::parse_expr(.iid_precision),
                data = data_mask,
                env = .envir
              )
              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
            )
          }
        } else {
          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.'
            )
          }
        }

        as.matrix(.values)
      }
      eval_fun
    }
  eval_list <- list()
  for (nm in names(eval_names)) {
    eval_list[[eval_names[[nm]]]] <-
      eval_fun_factory(
        comp_lst[[nm]],
        .envir = envir,
        .enclos = enclos
      )
  }

  # Remove problematic objects:
  problems <- c(".Random.seed")
  remove(list = intersect(names(envir), problems), envir = envir)

  n <- length(state)
  for (k in seq_len(n)) {
    state_df <- stats::setNames(state[[k]], state_names[names(state[[k]])])
    data_mask <- bru_data_mask(
      list(
        effect = effects[[k]],
        state_df,
        latent = state[[k]],
        data = data,
        data_extra = data_extra,
        eval_list,
        # Keep track of the iteration index so the iid cache can be
        # invalidated
        list(.cache_state_index = k)
      )
    )

    result_ <- rlang::eval_tidy(predictor, data = data_mask, env = envir)
    if (!is.null(n_pred) && is.numeric(result_) && length(result_) == 1) {
      result_ <- rep(result_, n_pred)
    }
    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
}