4  Example model file

All models in the package are defined as S3 classes and follow a strict template. This allows us to implement general methods for handling model fitting, data checking, and post-processing. Each model has an internal function that defines the model and its parameters, and a user-facing alias. Let’s look at how two models are implemented - the IMM model, which uses both general class and specific model methods, but no custom stan code, and the SDM model, which depends heavily on custom stan code. If you use the use_model_template() function, templates for all sections below will be automatically generated for your model.

4.1 The Interference Measurement Model (IMM)

The model is defined in the file R/model_imm.R. Let’s go through the different parts.

4.1.1 Model definition

The full IMM model is defined in the following internal model class:

R/model_imm.R

.model_imm <-
  function(resp_error = NULL, nt_features = NULL, nt_distances = NULL,
           set_size = NULL, regex = FALSE, version = "full", links = NULL,
           call = NULL, ...) {
    out <- structure(
      list(
        resp_vars = nlist(resp_error),
        other_vars = nlist(nt_features, nt_distances, set_size),
        domain = "Visual working memory",
        task = "Continuous reproduction",
        name = "Interference measurement model by Oberauer and Lin (2017).",
        version = version,
        citation = glue(
          "Oberauer, K., & Lin, H.Y. (2017). An interference model \\
          of visual working memory. Psychological Review, 124(1), 21-59"
        ),
        requirements = glue(
          '- The response vairable should be in radians and \\
          represent the angular error relative to the target
          - The non-target features should be in radians and be \\
          centered relative to the target'
        ),
        parameters = list(
          mu1 = glue(
            "Location parameter of the von Mises distribution for memory \\
            responses (in radians). Fixed internally to 0 by default."
          ),
          kappa = "Concentration parameter of the von Mises distribution",
          a = "General activation of memory items",
          c = "Context activation",
          s = "Spatial similarity gradient"
        ),
        links = list(
          mu1 = "tan_half",
          kappa = "log",
          a = "log",
          c = "log",
          s = "log"
        ),
        fixed_parameters = list(mu1 = 0, mu2 = 0, kappa2 = -100),
        default_priors = list(
          mu1 = list(main = "student_t(1, 0, 1)"),
          kappa = list(main = "normal(2, 1)", effects = "normal(0, 1)"),
          a = list(main = "normal(0, 1)", effects = "normal(0, 1)"),
          c = list(main = "normal(0, 1)", effects = "normal(0, 1)"),
          s = list(main = "normal(0, 1)", effects = "normal(0, 1)")
        )
      ),
      # attributes
      regex = regex,
      regex_vars = c('nt_features', 'nt_distances'),
      class = c("bmmodel", "circular", "non_targets", "imm", paste0('imm_',version)),
      call = call
    )

    # add version specific information
    if (version == "abc") {
      out$parameters$s <- NULL
      out$links$s <- NULL
      out$default_priors$s <- NULL
      attributes(out)$regex_vars <- c('nt_features')
    } else if (version == "bsc") {
      out$parameters$a <- NULL
      out$links$a <- NULL
      out$default_priors$a <- NULL
    }

    out$links[names(links)] <- links
    out
  }

Here is a brief explanation of the different components of the model definition:

resp_vars: a list of response variables that the model will be fitted to. These variables will be used to construct the brmsformula passed to brms together with the bmmformula and the parameters of the model. The user has to provide these variables in the data frame that is passed to the bmm() function

other_vars: a list of additional variables that are required for the model. This is used to check if the data contains all necessary information for fitting the model. In the example above, the IMM model requires the names of the variables specifying the non-target features relative to the target, the variables specifying the distance of the non-targets to the target, and the set_size. The user has to provide these variables in the data frame that is passed to the bmm() function

domain, task, name, citation, requirements: contains information about the model, such as the domain, task, name, citation, requirements. This information is used for generating help pages

version: if the model has multiple versions, this argument is specified by the user. Then it is used to dynamically adjust some information in the model object. In the case of the imm model, we have three versions - full, bsc and abc. As you can see at the end of the script, some parameters are deleted depending on the model version.

parameters: contains a named list of all parameters in the model that can be estimated by the user and their description. This information is used internally to check if the bmmformula contains linear model formulas for all model parameters, and to decide what information to include in the summary of bmmfit objects.

links: a named list providing the link function for each parameter. For example, kappa in the imm models has to be positive, so it is sampled on the log scale. This information is used in defining the model family and for the summary methods. If you want the user to be able to specify custom link functions, the next to last line of the script replaces the links with those provided by the user

fixed_parameters in the imm several parameters are fixed to constant values internally to identify the model. Only one of them, mu1 is also part of the parameters block - this is the only fixed parameters that users can choose to estimate instead of leaving it fixed. mu2 and kappa2 cannot be freely estimated.

default_priors a list of lists for each parameter in the model. Each prior has two components: main, the prior that will be put on the Intercept or on each level of a factor if the intercept is suppressed; effects, the prior to put on the regression coefficients relative to the intercept. The priors are described as in the set_prior function from brms. This information is used by the configure_prior() S3 method to automatically set the default priors for the model. The priors that you put here will be used by bmm() unless the users chooses to overwrite them.

regex: For the imm models, the nt_features and nt_distances variables can be specified with regular expressions, if the user sets regex = TRUE

call: this automatically records how the model was called so that the call can be printed in the summary after fitting. Leave it as is.

class: is the most important part. It contains the class of the model. This is used by generic S3 methods to perform data checks and model configuration. The classes should be ordered from most general to most specific. A general class exists when the same operations can be performed on multiple models. For example, the ‘3p’, ‘imm_abc’, ‘imm_bsc’ and ‘imm_full’ models all have non-targets and set_size arguments, so the same data checks can be performed on all of them, represented by the class non_targets. The first class should always be bmmodel, which is the main class for all models. The last class should be the specific model name, in this case imm_full, imm_abc or imm_bsc, which is automatically constructed if a version argument is provided. Otherwise the last class will be just the name of the model.

4.1.2 Model alias

The model alias is a user-facing function that calls the internal model function. It is defined as follows:

R/model_imm.R

#' @title `r .model_imm()$name`
#' @description Three versions of the `r .model_imm()$name` - the full, bsc, and abc.
#' `IMMfull()`, `IMMbsc()`, and `IMMabc()` are deprecated and will be removed in the future.
#' Please use `imm(version = 'full')`, `imm(version = 'bsc')`, or `imm(version = 'abc')` instead.
#'
#' @name imm
#' @details `r model_info(.model_imm(), components =c('domain', 'task', 'name', 'citation'))`
#' #### Version: `full`
#' `r model_info(.model_imm(version = "full"), components = c('requirements', 'parameters', 'fixed_parameters', 'links', 'prior'))`
#' #### Version: `bsc`
#' `r model_info(.model_imm(version = "bsc"), components = c('requirements', 'parameters', 'fixed_parameters', 'links', 'prior'))`
#' #### Version: `abc`
#' `r model_info(.model_imm(version = "abc"), components =c('requirements', 'parameters', 'fixed_parameters', 'links', 'prior'))`
#'
#' Additionally, all imm models have an internal parameter that is fixed to 0 to
#' allow the model to be identifiable. This parameter is not estimated and is not
#' included in the model formula. The parameter is:
#'
#'   - b = "Background activation (internally fixed to 0)"
#'
#' @param resp_error The name of the variable in the provided dataset containing
#'   the response error. The response Error should code the response relative to
#'   the to-be-recalled target in radians. You can transform the response error
#'   in degrees to radian using the `deg2rad` function.
#' @param nt_features A character vector with the names of the non-target
#'   variables. The non_target variables should be in radians and be centered
#'   relative to the target. Alternatively, if regex=TRUE, a regular
#'   expression can be used to match the non-target feature columns in the
#'   dataset.
#' @param nt_distances A vector of names of the columns containing the distances
#'   of non-target items to the target item. Alternatively, if regex=TRUE, a regular
#'   expression can be used to match the non-target distances columns in the
#'   dataset. Only necessary for the `bsc` and `full` versions.
#' @param set_size Name of the column containing the set size variable (if
#'   set_size varies) or a numeric value for the set_size, if the set_size is
#'   fixed.
#' @param regex Logical. If TRUE, the `nt_features` and `nt_distances` arguments
#'   are interpreted as a regular expression to match the non-target feature
#'   columns in the dataset.
#' @param version Character. The version of the IMM model to use. Can be one of
#'  `full`, `bsc`, or `abc`. The default is `full`.
#' @param ... used internally for testing, ignore it
#' @return An object of class `bmmodel`
#' @keywords bmmodel
#' @examplesIf isTRUE(Sys.getenv("BMM_EXAMPLES"))
#' # load data
#' data <- oberauer_lin_2017
#'
#' # define formula
#' ff <- bmmformula(
#'   kappa ~ 0 + set_size,
#'   c ~ 0 + set_size,
#'   a ~ 0 + set_size,
#'   s ~ 0 + set_size
#' )
#'
#' # specify the full IMM model with explicit column names for non-target features and distances
#' # by default this fits the full version of the model
#' model1 <- imm(resp_error = "dev_rad",
#'               nt_features = paste0('col_nt', 1:7),
#'               nt_distances = paste0('dist_nt', 1:7),
#'               set_size = 'set_size')
#'
#' # fit the model
#' fit <- bmm(formula = ff,
#'            data = data,
#'            model = model1,
#'            cores = 4,
#'            backend = 'cmdstanr')
#'
#' # alternatively specify the IMM model with a regular expression to match non-target features
#' # this is equivalent to the previous call, but more concise
#' model2 <- imm(resp_error = "dev_rad",
#'               nt_features = 'col_nt',
#'               nt_distances = 'dist_nt',
#'               set_size = 'set_size',
#'               regex = TRUE)
#'
#' # fit the model
#' fit <- bmm(formula = ff,
#'            data = data,
#'            model = model2,
#'            cores = 4,
#'            backend = 'cmdstanr')
#'
#' # you can also specify the `bsc` or `abc` versions of the model to fit a reduced version
#' model3 <- imm(resp_error = "dev_rad",
#'               nt_features = 'col_nt',
#'               set_size = 'set_size',
#'               regex = TRUE,
#'               version = 'abc')
#' fit <- bmm(formula = ff,
#'            data = data,
#'            model = model3,
#'            cores = 4,
#'            backend = 'cmdstanr')
#' @export
imm <- function(resp_error, nt_features, nt_distances, set_size, regex = FALSE, version = "full", ...) {
  call <- match.call()
  dots <- list(...)
  if ("setsize" %in% names(dots)) {
    set_size <- dots$setsize
    warning2("The argument 'setsize' is deprecated. Please use 'set_size' instead.")
  }
  if (version == "abc") nt_distances <- NULL
  stop_missing_args()
  .model_imm(resp_error = resp_error, nt_features = nt_features,
             nt_distances = nt_distances, set_size = set_size, regex = regex,
             version = version, call = call, ...)
}

The details will be filled out automatically from the model definition. This does some fancy formatting to include documentation about all versions of the model in the same help file.

4.1.3 check_data() methods

Each model should have a check_data.modelname() method that checks if the data contains all necessary information for fitting the model. For the IMM, the bsc and full versions require a special check for the nt_distances variables:

R/model_imm.R

#' @export
check_data.imm_bsc <- function(model, data, formula) {
  data <- .check_data_imm_dist(model, data, formula)
  NextMethod("check_data")
}

#' @export
check_data.imm_full <- function(model, data, formula) {
  data <- .check_data_imm_dist(model, data, formula)
  NextMethod("check_data")
}

.check_data_imm_dist <- function(model, data, formula) {
  nt_distances <- model$other_vars$nt_distances
  max_set_size <- attr(data, "max_set_size")

  stopif(
    !isTRUE(all.equal(length(nt_distances), max_set_size - 1)),
    "The number of columns for non-target distances in the argument \\
    'nt_distances' should equal max(set_size)-1})"
  )

  # replace NA values with 999 so they have 0 effect through the distance formula
  data[, nt_distances][is.na(data[, nt_distances])] <- 999

  stopif(
    any(data[, nt_distances] < 0),
    "All non-target distances to the target need to be postive."
  )
  data
}

The IMM models share methods with the mixture3p model, all of which are of class non_targets so the check_data.non_targets method is defined in the general file R/helpers-data.R. If you are adding a new model, you should check if the data requirements are similar to any existing model and define the check_data method only for the methods that are unique to your model.

The check_data.mymodel() function should always take the arguments model, data, and formula and return the data with the necessary transformations. It should also call NextMethod("check_data") to call the check_data method of the more general class.

To understand what validation your model inherits, trace the NextMethod dispatch chain for your class. For example, the imm_full model with class c("bmmodel", "circular", "non_targets", "imm", "imm_full") dispatches as:

1. check_data.bmmodel — generic checks (data is coercible to a data frame and has rows), then sets the data_name and checked attributes
2. → NextMethod → check_data.circular — validates response is in radians, wraps to \((-\pi, \pi]\)
3. → NextMethod → check_data.non_targets — validates nt_features, set_size, creates index variables, stores data attributes
4. → NextMethod → no check_data method is defined for the imm class, so the chain passes straight through it
5. → NextMethod → check_data.imm_full — validates nt_distances (model-specific)

The order is general first, specific last — the reverse of what most people expect. S3 dispatch walks the class vector from left to right, and "bmmodel" is the leftmost class, so check_data.bmmodel runs before check_data.imm_full, not after. A class with no method for the generic is simply skipped.

This matters for two reasons. Your model-specific method runs last, so it can rely on everything the general methods have already done: by the time check_data.imm_full runs, the response has been wrapped to radians and max_set_size is already an attribute on the data. And it means a fall-through fallback must go on .default, not on .bmmodel — see Chapter 6.

This also means a new circular model with non-targets automatically gets radian validation, non-target checks, and column existence checks — you only need to write methods for your model-specific validations.

4.1.4 configure_model() methods

The configure_model.mymodel() method is where you specify the model formula, the family, any custom code. The method is defined as follows for the IMM model:

(we show only the IMMfull version)

R/model_imm.R

configure_model.imm_full <- function(model, data, formula) {
  # retrieve arguments from the data check
  max_set_size <- attr(data, "max_set_size")
  lure_idx <- attr(data, "lure_idx_vars")
  nt_features <- model$other_vars$nt_features
  nt_distances <- model$other_vars$nt_distances

  formula <- .base_imm_formula(model, formula) +
    brms::nlf(theta1 ~ log(exp(c) + exp(a))) +
    brms::nlf(expS ~ exp(s))

  theta_exprs <- glue("{lure_idx} * log(exp(c-expS*{nt_distances}) + exp(a)) + (1 - {lure_idx}) * (-100)")
  formula <- .add_imm_mixture_terms(formula, max_set_size, lure_idx, nt_features, theta_exprs)
  formula$family <- .imm_mixture_family(max_set_size)

  nlist(formula, data)
}

The three IMM versions differ only in the theta expression and in which extra non-linear terms they add, so everything they share lives in three private helpers at the top of R/model_imm.R:

R/model_imm.R

.base_imm_formula <- function(model, formula) {
  bmf2bf(model, formula) +
    brms::lf(kappa2 ~ 1) +
    brms::lf(mu2 ~ 1) +
    brms::nlf(kappa1 ~ kappa)
}

.add_imm_mixture_terms <- function(formula, max_set_size, lure_idx, nt_features, theta_exprs) {
  kappa_nts <- paste0("kappa", 3:(max_set_size + 1))
  theta_nts <- paste0("theta", 3:(max_set_size + 1))
  mu_nts <- paste0("mu", 3:(max_set_size + 1))

  for (i in 1:(max_set_size - 1)) {
    formula <- formula +
      glue_nlf("{kappa_nts[i]} ~ kappa") +
      glue_nlf("{theta_nts[i]} ~ {theta_exprs[i]}") +
      glue_nlf("{mu_nts[i]} ~ {nt_features[i]}")
  }

  formula
}

.imm_mixture_family <- function(max_set_size) {
  brms::mixture(
    brms::von_mises("tan_half"), brms::von_mises("identity"),
    nmix = c(1, max_set_size),
    order = "none"
  )
}

The version-specific method then builds theta_exprs as a glue vector — one expression per non-target position — and hands it to .add_imm_mixture_terms(), which loops. imm_abc passes "{lure_idx} * a + (1 - {lure_idx}) * (-100)", imm_full the distance-weighted expression above. This is the pattern to copy if your model has versions that share a formula skeleton: put the skeleton in a helper prefixed with a dot, and let each version supply only what differs. Note the dot prefix, but not .model_ — supported_models() scans for that prefix and would list your helper as a model.

The configure_model method should always take the arguments model, data, and formula (as a bmmformula) and return a named list with the formula (as a brmsformula) and the data. The brmsfamily should be stored within the formula.

Inside the configure_model method the brmsformula is generated using the bmf2bf function. This function converts the bmmformula passed to bmm() function into a brmsformula based on the information for the response variables provided in the bmmmodel object. There is a general method in R/bmmformula.R to construct the formula for all models with a single response variable.

R/bmmformula.R

# If no model specific methods exist, it will construct a base brms formula from the first resp_vars
# We do it this way because most models require just one response variable
#' @export
bmf2bf.bmmodel <- function(model, formula = bmmformula()) {
  brms_formula <- NextMethod("bmf2bf") %||% brms::bf(glue("{model$resp_vars[[1]]} ~ 1"))
  components <- lapply(formula, function(x) if (is_nl(x)) brms::nlf(x) else brms::lf(x))
  Reduce(`+`, components, init = brms_formula)
}

#' @export
bmf2bf.default <- function(model, formula) {
  NULL
}

Read the dispatch order carefully here, because it is the reverse of what the names suggest. bmf2bf() is called with the model object, whose class vector starts with "bmmodel", so bmf2bf.bmmodel runs first. Its NextMethod("bmf2bf") walks on to the rest of the chain, and only reaches bmf2bf.default if no model-specific method intercepted.

That is why bmf2bf.default returns NULL rather than a formula. It means “nobody claimed the first line”. bmf2bf.bmmodel sees the NULL, falls through the %||%, and builds the single-response first line resp ~ 1 itself. A model with several response variables, or one needing vreal(), vint() or trials(), defines bmf2bf.mymodel that returns the first line it wants; NextMethod() then returns that instead of NULL and the %||% fallback never fires.

Once the first line exists, the rest is mechanical: each parameter formula becomes brms::nlf() if it is marked non-linear and brms::lf() otherwise, and Reduce() adds them all to the base formula. Whether a component is marked non-linear is decided in assign_nl_attr() (R/bmmformula.R:375-382) purely by whether one of its predictors is also a predicted parameter — which is the mechanism behind the name-clash warning described in the next chapter.

This conversion from a bmmformula object into a brmsformula object is done to avoid users having to specify complicated and long formulas or specifying all additional response information in the brmsformula themselves. For more detailed information on the use of additional response information in a brmsformula please see the brmsformula documentation.

4.2 The Signal Discrimination Model (SDM)

The SDM model is defined in the file R/model_sdm.R. The SDM model differs in the configuration compared to the IMM model, as it requires custom STAN code. Let’s go through the different parts. As before, we start with the model definition.

4.2.1 Model definition

R/model_sdm.R

.model_sdm <- function(resp_error = NULL, links = NULL, version = "simple", call = NULL, ...) {
  out <- structure(
    list(
      resp_vars = nlist(resp_error),
      other_vars = nlist(),
      domain = "Visual working memory",
      task = "Continuous reproduction",
      name = "Signal Discrimination Model (SDM) by Oberauer (2023)",
      citation = glue(
        "Oberauer, K. (2023). Measurement models for visual working memory - \\
        A factorial model comparison. Psychological Review, 130(3), 841-852"
      ),
      version = version,
      requirements = glue(
        "- The response variable should be in radians and represent the angular \\
        error relative to the target"
      ),
      parameters = list(
        mu = glue("Location parameter of the SDM distribution (in radians; \\
                  by default fixed internally to 0)"),
        c = "Memory strength parameter of the SDM distribution",
        kappa = "Precision parameter of the SDM distribution"
      ),
      links = list(
        mu = "tan_half",
        c = "log",
        kappa = "log"
      ),
      fixed_parameters = list(mu = 0),
      default_priors = list(
        mu = list(main = "student_t(1, 0, 1)"),
        kappa = list(main = "student_t(5, 1.75, 0.75)", effects = "normal(0, 1)"),
        c = list(main = "student_t(5, 2, 0.75)", effects = "normal(0, 1)")
      ),
      init_ranges = list(
        mu = c(-0.5, 0.5),
        kappa = c(2.5, 3.5),
        c = c(4, 6)
      )
    ),
    class = c("bmmodel", "circular", "sdm", paste0("sdm_", version)),
    call = call
  )
  out$links[names(links)] <- links
  out
}

Note that compared to the earlier IMM example, the SDM model now includes init_ranges. This provides initial value ranges for MCMC sampling, which the create_initfun() step uses to generate appropriate starting values. The SDM model benefits from custom initialization because its custom Stan likelihood can be sensitive to extreme starting values.

The model definition is similar to the IMM model, but the SDM model only requires the user to specify the response error, but not additional variables such as non-target variables. The class is also different, as the SDM model is not a subclass of the IMM model. We’ll skip the alias for the SDM model, as it is similar for every model.

4.2.2 check_data() methods

The SDM shares a class with other circular models, so most of the data checks are performed by check_data.circular method, defined in the general file R/helpers-data.R. The sdm however, samples much more quickly in Stan, if the data is sorted by the predictor variables, so we have the following custom data check method for the sdm:

R/model_sdm.R

#' @export
check_data.sdm <- function(model, data, formula) {
  # data sorted by predictors is necessary for speedy computation of normalizing constant
  data <- order_data_query(model, data, formula)
  attr(data, "sdm_run_metadata") <- sdm_run_metadata(data, formula, model)
  NextMethod("check_data")
}

This is a compact example of the attribute bridge. Having sorted the data, the method computes where each run of identical predictor values starts and how long it is — the information the Stan likelihood needs to compute its normalizing constant once per run rather than once per row — and hangs the result on the data frame as sdm_run_metadata. configure_model.sdm reads it back out a few pipeline steps later. The computation itself is a separate function, sdm_run_metadata(), and that matters: it is the reason configure_model.sdm can recompute the metadata when it is called directly and the attribute is absent.

One detail in that function is worth copying. brms drops rows with missing values in any model variable before fitting, so run boundaries computed on the full data frame would not line up with the rows Stan actually sees. sdm_run_metadata() therefore applies stats::complete.cases() over the response and predictor columns before computing the runs (R/model_sdm.R:202-211). Any index you precompute in check_data and pass through to Stan has the same problem.

4.2.3 configure_model() methods

The configure_model method for the SDM model is different compared to the IMM model, as it requires custom STAN code. The method is defined as follows:

R/model_sdm.R

#' @export
configure_model.sdm <- function(model, data, formula) {
  # note - c has a log link, but I've coded it manually for computational efficiency
  sdm_simple <- brms::custom_family(
    name = "sdm_simple",
    dpars = c("mu", "c", "kappa"),
    links = c("tan_half", "identity", "log"),
    lb = c(NA, NA, NA),
    ub = c(NA, NA, NA),
    type = "real", loop = FALSE,
    log_lik = log_lik_sdm_simple,
    posterior_predict = posterior_predict_sdm_simple
  )

  # prepare initial stanvars to pass to brms, model formula and priors
  sc_path <- system.file("stan_chunks", package = "bmm")
  stan_funs <- read_lines2(paste0(sc_path, "/sdm_simple_funs.stan"))
  stan_tdata <- read_lines2(paste0(sc_path, "/sdm_simple_tdata.stan"))
  likelihood_file <- if (sdm_use_threaded_likelihood()) {
    "sdm_simple_likelihood_threaded.stan"
  } else {
    "sdm_simple_likelihood.stan"
  }
  stan_likelihood <- read_lines2(paste0(sc_path, "/", likelihood_file))
  stan_tdata_pll_args <- if (sdm_use_threaded_likelihood()) {
    "data matrix COSN"
  }
  run_metadata <- attr(data, "sdm_run_metadata")
  if (is.null(run_metadata)) {
    # Guard direct configure_model.sdm() calls that bypass check_data.sdm().
    run_metadata <- sdm_run_metadata(data, formula, model)
  }
  stanvars <- brms::stanvar(scode = stan_funs, block = "functions") +
    brms::stanvar(scode = stan_tdata, block = "tdata", pll_args = stan_tdata_pll_args) +
    brms::stanvar(x = run_metadata$G_sdm_runs, name = "G_sdm_runs") +
    sdm_stanvar_int_array(run_metadata$sdm_run_start, "sdm_run_start", "G_sdm_runs") +
    sdm_stanvar_int_array(run_metadata$sdm_run_count, "sdm_run_count", "G_sdm_runs") +
    brms::stanvar(scode = stan_likelihood, block = "likelihood", position = "end")

  # construct main brms formula from the bmm formula
  formula <- bmf2bf(model, formula)
  formula$family <- sdm_simple

  # return the list
  nlist(formula, data, stanvars)
}

Note that configure_model no longer returns init — initial values are now handled separately by create_initfun(), which uses the init_ranges defined in the model object.

This is the longest configure_model method in the package, so it is worth separating the parts that every custom-family model needs from the parts that are specific to SDM.

Every custom-family model needs the first two. The brms::custom_family() call declares the parameters (dpars), their link functions, bounds, whether the likelihood is vectorized (loop = FALSE here), and the R functions that implement log_lik and posterior_predict. For more information, see the brms custom families vignette. And the Stan code is read from inst/stan_chunks/ through system.file() — not a relative path — so that it is found once the package is installed rather than only when it is loaded from source. Each chunk becomes its own brms::stanvar(), tagged with the Stan block it belongs in.

The rest is SDM-specific, and all of it exists to make the normalizing constant cheap. Three things are going on:

• Data passed to Stan, not just code. brms::stanvar(x = ..., name = ...) passes an R value into the Stan data block. SDM passes the run count and, through the helper sdm_stanvar_int_array(), the run start and length arrays computed back in check_data.sdm. The helper exists because stanvar() would declare an integer array as vector, so it overwrites the generated scode with array[G_sdm_runs] int <name>;.
• Two likelihood chunks. Which one is read depends on sdm_use_threaded_likelihood(). See the threading section of the architecture chapter for why there are two and what threading(force = TRUE) does.
• pll_args. Data that the threaded likelihood needs must also be declared in the partial-log-likelihood signature. That is what pll_args = "data matrix COSN" on the tdata stanvar does, and what sdm_stanvar_int_array() sets for the index arrays. Without it the variable is out of scope inside partial_log_lik and the model does not compile.

Converting the bmmformula to a brmsformula and collecting all arguments is done entirely using the bmf2bf method.

4.2.4 Post-processing methods

Unlike the imm model, the sdm model requires some special post-processing because of the way the link functions are coded. These methods are applied after the brmsfit object is returned, at the very end of the bmm() pipeline:

R/model_sdm.R

#' @export
postprocess_brm.sdm <- function(model, fit, ...) {
  # manually set link_c to "log" since I coded it manually
  fit$family$link_c <- "log"
  fit$formula$family$link_c <- "log"
  fit
}

#' @export
revert_postprocess_brm.sdm <- function(model, fit, ...) {
  fit$family$link_c <- "identity"
  fit$formula$family$link_c <- "identity"
  fit
}

The revert_postprocess_brm method is the inverse of postprocess_brm. It is called internally when functions like conditional_effects() or emmeans() need to temporarily undo post-processing to work with the model in its original brms representation. If your postprocess_brm method modifies the fitted object (e.g., changing link functions, renaming parameters), you should also define a revert_postprocess_brm method that reverses those changes.

We also have a couple of special functions for custom families in brms (see the log_likand posterior_predict argument in the call to brms::custom_familiy()), which allow other typical tools from brms such posterior_predict of bridgesampling to work:

R/model_sdm.R

log_lik_sdm_simple <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  c <- brms::get_dpar(prep, "c", i = i)
  kappa <- brms::get_dpar(prep, "kappa", i = i)
  y <- prep$data$Y[i]
  dsdm(y, mu, c, kappa, log = T)
}

posterior_predict_sdm_simple <- function(i, prep, ...) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  c <- brms::get_dpar(prep, "c", i = i)
  kappa <- brms::get_dpar(prep, "kappa", i = i)
  rsdm(length(mu), mu, c, kappa)
}

4.3 The Diffusion Decision Model (DDM)

The DDM model is defined in R/model_ddm.R.

4.3.1 Model definition

It demonstrates several patterns that differ from the circular models above: suppressing mu for a model with no natural location parameter, init_ranges for MCMC initialization, and a custom bmf2bf method for response variable wrapping.

The DDM uses a defaults table to store parameters, links, priors, and init_ranges in one place:

R/model_ddm.R

.ddm_defaults <- list(
  parameters = list(
    drift = "Drift rate = Average rate of evidence accumulation of the decision processes",
    bound = "Boundary separation = Distance between the decision boundaries that need to be reached",
    ndt   = "Non-decision time = Additional time required beyond the evidence accumulation process",
    zr    = "Relative starting point = Starting point between the decision thresholds relative to the upper bound."
  ),
  links = list(
    drift = "identity", bound = "log", ndt = "log", zr = "logit"
  ),
  fixed_parameters = list(
    zr = 0,
    mu = 0
  ),
  priors = list(
    drift = list(main = "cauchy(0,1)", effects = "normal(0,0.5)"),
    bound = list(main = "normal(0,0.5)", effects = "normal(0,0.5)"),
    ndt   = list(main = "normal(-1.5,0.5)", effects = "normal(0,0.3)"),
    zr    = list(main = "normal(0,0.5)", effects = "normal(0,0.3)")
  ),
  init_ranges = list(
    mu = c(-0.1, 0.1),
    drift  = c(-1, 1),
    bound  = c(1, 2),
    ndt    = c(0.01, 0.05),
    zr     = c(0.45, 0.55)
  )
)

.model_ddm <- function(rt = NULL, response = NULL, links = NULL, call = NULL, ...) {
  out <- structure(
    list(
      resp_vars = nlist(rt, response),
      other_vars = nlist(),
      domain = "Decision Making / Response times",
      task = "Two-Alternative Force Choice RT",
      name = "Diffusion Decision Model",
      version = "NA",
      citation = glue(
        "Ratcliff, R. (1978). A theory of memory retrieval. Psychological Review, 85(2), 59-108. https://doi.org/10/fjwm2f;"
      ),
      requirements = glue(
        "- The response time should be in seconds and \\
          represent the time between onset of the target stimulus until the response execution
          - The response should be coded numerically: \\
          0 = lower response, 1 = upper response"
      ),
      parameters = .ddm_defaults[["parameters"]],
      links = .ddm_defaults[["links"]],
      fixed_parameters = .ddm_defaults[["fixed_parameters"]],
      default_priors = .ddm_defaults[["priors"]],
      init_ranges = .ddm_defaults[["init_ranges"]]
    ),
    class = c("bmmodel", "ddm"),
    call = call
  )
  out$links[names(links)] <- links
  out
}

Key differences from the circular models:

• No natural mu: The DDM has no location parameter. A dummy mu is included in the family’s dpars and fixed to 0 via fixed_parameters, which satisfies brms’ requirement that every custom family has a mu parameter.
• Two response variables: resp_vars = nlist(rt, response) — reaction time and binary choice. This requires a custom bmf2bf method.
• init_ranges: Specifies plausible starting ranges for each parameter to ensure stable MCMC sampling. The create_initfun() step uses these.
• Short class chain: c("bmmodel", "ddm") — no intermediate domain class. RT models don’t share as many check_data patterns as the circular models do.

4.3.2 User-facing alias

The DDM alias follows the standard pattern — calling stop_missing_args() and then the internal constructor:

R/model_ddm.R

ddm <- function(rt, response, links = NULL, ...) {
  call <- match.call()
  stop_missing_args()
  .model_ddm(rt = rt, response = response, links = links, call = call, ...)
}

4.3.3 bmf2bf method

Because the DDM has two response variables (rt and response), it needs a custom bmf2bf method to construct the first line of the brms formula with the proper response wrapping:

R/model_ddm.R

#' @export
bmf2bf.ddm <- function(model, formula) {
  rt <- model$resp_vars$rt
  response <- model$resp_vars$response
  brms::bf(paste0(rt, " | dec(", response, ") ~ 1"))
}

This creates a formula like rt | dec(response) ~ 1, where dec() tells brms which variable codes the binary decision (0 = lower, 1 = upper boundary).

4.3.4 configure_model method

R/model_ddm.R

#' @export
configure_model.ddm <- function(model, data, formula) {
  formula <- bmf2bf(model, formula)

  ddm_family <- function(link_drift, link_bound, link_ndt, link_zr) {
    brms::custom_family(
      "ddm",
      dpars = c("mu", "drift", "bound", "ndt", "zr"),
      links = c("identity", link_drift, link_bound, link_ndt, link_zr),
      lb = c(NA, NA, 0.1, 0, 0),
      ub = c(NA, NA, NA, NA, 1),
      type = "real",
      vars = "dec[n]",
      loop = TRUE,
      log_lik = log_lik_ddm,
      posterior_predict = posterior_predict_ddm
    )
  }

  formula$family <- ddm_family(
    link_drift = model$links$drift,
    link_bound = model$links$bound,
    link_ndt = model$links$ndt,
    link_zr = model$links$zr
  )

  sc_path <- system.file("stan_chunks", package = "bmm")
  stan_functions <- read_lines2(paste0(sc_path, "/ddm_functions.stan"))
  stanvars <- brms::stanvar(scode = stan_functions, block = "functions")

  nlist(formula, data, stanvars)
}

Note that:

• mu is included in dpars with an identity link and fixed to 0 (the mu-suppression pattern)
• vars = "dec[n]" passes the decision variable to the Stan likelihood
• The log_lik and posterior_predict functions use the dddm() and rddm() distribution functions from R/distributions.R
• configure_model returns nlist(formula, data, stanvars) without init — initialization is handled by create_initfun()

4.3.5 check_data method

The DDM has thorough data validation for RT and response variables:

R/model_ddm.R

#' @export
check_data.ddm <- function(model, data, formula) {
  rt_var <- model$resp_vars$rt
  response_var <- model$resp_vars$response

  # check variables exist
  stopif(not_in(rt_var, colnames(data)),
         "The RT variable '{rt_var}' is not present in the data.")
  stopif(not_in(response_var, colnames(data)),
         "The response variable '{response_var}' is not present in the data.")

  # remove NAs with warning
  if (any(is.na(data[, rt_var]))) {
    data <- data[!is.na(data[, rt_var]), ]
    warning2("Some values in {rt_var} were NA. These were removed.")
  }

  # validate RT: must be positive, in seconds (warn if > 10s or < 0.1s)
  stopif(any(data[, rt_var] < 0, na.rm = TRUE),
         "Some reaction times are lower than zero.")
  warnif(any(data[, rt_var] > 10, na.rm = TRUE),
         "Your data contains reaction times larger than 10 seconds...")

  # validate response: accepts integer (0/1), logical, or character ("upper"/"lower")
  # coerces all formats to 0/1 integer
  # ... (full validation with type coercion)

  NextMethod("check_data")
}

The DDM’s check_data is a good example of thorough validation: it checks variable existence, removes NAs with warnings, validates RT ranges, and flexibly accepts multiple response coding formats (integer, logical, character) while coercing them all to 0/1.

4.4 The EZ-Diffusion Model (EZDM)

The EZDM model (R/model_ezdm.R) demonstrates two additional patterns: the version table pattern for managing multiple model versions, and working with aggregate-level data (summary statistics rather than trial-level observations).

4.4.1 Version table

Instead of using if/else logic to handle version-specific parameters (as the IMM does), the EZDM uses a lookup table:

R/model_ezdm.R

.ezdm_version_table <- list(
  "3par" = list(
    parameters = list(
      drift = "Drift rate = ...",
      bound = "Boundary separation = ...",
      ndt = "Non-decision time = ...",
      s = "The diffusion constant = ..."
    ),
    links = list(drift = "identity", bound = "log", ndt = "log", s = "log"),
    fixed_parameters = list(s = 0, mu = 0),
    priors = list(
      drift = list(main = "cauchy(0,1)", effects = "normal(0,0.5)"),
      bound = list(main = "normal(0,0.5)", effects = "normal(0,0.5)"),
      ndt = list(main = "normal(-1.5,0.5)", effects = "normal(0,0.3)"),
      s = list(main = "normal(0,1)", effects = "normal(0,0.3)")
    ),
    init_ranges = list(
      mu = c(0, 1), drift = c(-1, 1), bound = c(1, 2),
      ndt = c(0.25, 0.5), s = c(0.99, 1.01)
    )
  ),
  "4par" = list(
    parameters = list(...),  # adds zr parameter
    links = list(...),       # adds zr = "logit"
    ...
  )
)

The model constructor then indexes into this table:

.model_ezdm <- function(mean_rt = NULL, var_rt = NULL, n_upper = NULL,
                         n_trials = NULL, version = "3par", links = NULL,
                         call = NULL, ...) {
  out <- structure(
    list(
      resp_vars = nlist(mean_rt, var_rt, n_upper),
      other_vars = nlist(n_trials),
      ...
      parameters = .ezdm_version_table[[version]][["parameters"]],
      links = .ezdm_version_table[[version]][["links"]],
      fixed_parameters = .ezdm_version_table[[version]][["fixed_parameters"]],
      default_priors = .ezdm_version_table[[version]][["priors"]],
      init_ranges = .ezdm_version_table[[version]][["init_ranges"]]
    ),
    class = c("bmmodel", "ezdm"),
    call = call
  )
  if (!is.null(version)) class(out) <- c(class(out), paste0("ezdm_", version))
  ...
}

This pattern is cleaner and easier to extend than conditional logic. New models with multiple versions should prefer this approach.

4.4.2 bmf2bf for aggregate data

The EZDM operates on aggregate statistics (mean RT, variance of RT, number of upper responses, total trials) rather than trial-level data. This requires version-specific bmf2bf methods that pack multiple response variables into brms’ special response terms:

R/model_ezdm.R

#' @export
bmf2bf.ezdm_3par <- function(model, formula) {
  mean_rt <- model$resp_vars$mean_rt
  var_rt <- model$resp_vars$var_rt
  n_upper <- model$resp_vars$n_upper
  n_trials <- model$other_vars$n_trials

  brms::bf(paste0(mean_rt, " | vreal(", var_rt, ") + vint(", n_upper,
                  ") + trials(", n_trials, ") ~ 1"))
}

This creates a formula like mean_rt | vreal(var_rt) + vint(n_upper) + trials(n_trials) ~ 1, where:

• vreal() passes real-valued additional response information to Stan
• vint() passes integer-valued additional response information
• trials() passes the trial count

These special brms terms make the additional variables available in the Stan likelihood via vreal1[n], vint1[n], and trials[n].

4.4.3 Version-specific configure_model

Each EZDM version has its own configure_model method because the Stan likelihood functions differ:

R/model_ezdm.R

#' @export
configure_model.ezdm_3par <- function(model, data, formula) {
  formula <- bmf2bf(model, formula)
  links <- model$links

  formula$family <- brms::custom_family(
    "ezdm_3par",
    dpars = c("mu", "drift", "bound", "ndt", "s"),
    links = c("identity", links$drift, links$bound, links$ndt, links$s),
    lb = c(NA, NA, 0, 0, 0),
    ub = c(NA, NA, NA, NA, NA),
    type = "real",
    log_lik = log_lik_ezdm_3par,
    posterior_predict = posterior_predict_ezdm_3par,
    loop = TRUE,
    vars = c("vreal1[n]", "vint1[n]", "trials[n]")
  )

  sc_path <- system.file("stan_chunks", package = "bmm")
  stan_functions <- read_lines2(paste0(sc_path, "/ezdm_3par_functions.stan"))
  stanvars <- brms::stanvar(scode = stan_functions, block = "functions")

  nlist(formula, data, stanvars)
}

Note the vars argument in custom_family — it lists the additional response terms that the Stan likelihood needs access to.

---

We will now look at how to construct all these parts for a new model. Hint: you don’t have to do it manually, you can use the use_model_template() function to generate templates for your model.