R/fit_inla_model.R

In mbg: Model-Based Geostatistics

#' Fit INLA model
#'
#' @description Fit an INLA model based on a constructed data stack and formula
#'
#' @details Using [INLA::inla()] with reasonable defaults and settings tuned to predict
#'   across a grid.
#'
#' @param formula (character) INLA formula to fit. Generated in
#'   [prepare_inla_data_stack()], then interpreted using [stats::as.formula()] within the
#'   call to [INLA::inla()].
#' @param data_stack Stacked data, covariates, and spatial index. Generated in
#'   [prepare_inla_data_stack()].
#' @param spde SPDE object generated by [prepare_inla_data_stack()].
#' @param samplesize_vec (`integer(N)`, default 1) Sample sizes for each data observation.
#'   Only used for binomial data models.
#' @param precision_vec (`numeric(N)`, default 1) Precision for each data observation.
#'   Only used for gaussian data models.
#' @param family (character, default 'binomial') GLM family to use. For more information,
#'   see [stats::family()].
#' @param link (character, default 'logit') Link function to use, typically related to the
#'   GLM `family`.
#' @param fixed_effects_pc_prior A named list specifying the penalized complexity prior
#'   for all fixed effects except for the intercept. The two named items are "threshold",
#'   the test threshold for the size of each fixed effect, and "prob_above", the prior
#'   probability that the beta for each covariate will EXCEED that threshold.
#' @param verbose (`logical(1)`, default TRUE) Log progress for INLA model fitting?
#'
#' @return A fitted INLA model object created by [INLA::inla()]
#'
#' @concept model_fit
#'
#' @seealso [mbg::MbgModelRunner]
#'
#' @importFrom stats as.formula
#' @export
fit_inla_model <- function(
  formula, data_stack, spde, samplesize_vec = 1, precision_vec = 1,
  family = 'binomial', link = 'logit',
  fixed_effects_pc_prior = list(threshold = 3, prob_above = 0.05),
  verbose = TRUE
){
  spde <- spde
  if(verbose) logging_start_timer("MBG model fitting")
  inla_model <- INLA::inla(
    formula = stats::as.formula(formula),
    family = family,
    control.family = list(link = link),
    Ntrials = if(family == 'binomial') samplesize_vec else NULL,
    scale = if(family == 'gaussian') precision_vec else NULL,
    data = INLA::inla.stack.data(data_stack),
    control.compute = list(
      config = TRUE,
      waic = TRUE
    ),
    control.fixed = list(
      prec = list(
        prior = 'pc.prec',
        param = c(fixed_effects_pc_prior$threshold, fixed_effects_pc_prior$prob_above)
      )
    ),
    control.predictor = list(
      compute = FALSE,
      link = 1,
      A = INLA::inla.stack.A(data_stack)
    )
  )
  if(verbose) logging_stop_timer()
  return(inla_model)
}