4 Example models

This chapter walks through two model files in the package. The IMM is built on a native brms family, brms::mixture(), and needs no Stan code of its own; it shows shared classes, check_data() dispatch and non-linear formulas. The DDM is a custom family with its own Stan code; it shows a flat specification, a model without a location parameter, and a custom first formula line. Between them they cover what use_model_template() asks you to fill in.

4.1 The Interference Measurement Model (IMM)

The IMM lives in R/model_imm.R.

4.1.1 Model definition

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).",
      citation = glue(
        "Oberauer, K., & Lin, H.-Y. (2017). An interference model \\
          of visual working memory. Psychological Review, 124(1), 21-59. \\
          https://doi.org/10.1037/rev0000044"
      ),
      version = version,
      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 = "normal(0, 0.5)", effects = "normal(0, 0.25)", sd = "exponential(4)"),
        kappa = list(main = "normal(2, 1)", effects = "normal(0, 1)", sd = "exponential(1)"),
        a = list(main = "normal(0, 1)", effects = "normal(0, 1)", sd = "exponential(1)"),
        c = list(main = "normal(0, 1)", effects = "normal(0, 1)", sd = "exponential(1)"),
        s = list(main = "normal(0, 1)", effects = "normal(0, 1)", sd = "exponential(1)")
      ),
      # central 50% of the main default prior on the native scale
      init_ranges = list(
        mu1 = c(-0.65, 0.65),
        kappa = c(3.8, 15),
        a = c(0.51, 2),
        c = c(0.51, 2),
        s = c(0.51, 2)
      )
    ),
    # 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
    out$init_ranges$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$init_ranges$a <- NULL
  }

  out <- set_links(out, links)
  out
}

The fields are described in Chapter 3. What the IMM adds:

Versions without a version table. The IMM predates version tables. It builds the full specification and then removes parameters for the bsc and abc versions. New models use a version table instead (Section A.6).
More response information than one column. other_vars names the non-target feature and distance columns and the set size; with regex = TRUE the user can give a pattern instead of the column names.
Fixed parameters that identify the model. mu2 and kappa2 are fixed internally. mu1 is in both fixed_parameters and parameters, so it is fixed by default and users may estimate it by giving it a formula.
Links that cannot be changed. The IMM writes its links into the mixture family and the non-linear formulas, so settable_links.imm() returns character(0) and a link a user passes is refused (Section A.1).

4.1.2 `check_data()` and the dispatch chain

The bsc and full versions need an extra check of the distance columns:

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
}

Everything else is inherited. For imm_full, with class c("bmmodel", "circular", "non_targets", "imm", "imm_full"), the methods run in this order:

check_data.bmmodel(): the data can be coerced to a data frame and has rows; sets the data_name and checked attributes.
check_data.circular(): the response is in radians, wrapped to $(-\pi, \pi]$.
check_data.non_targets(): checks nt_features and set_size, creates index variables, and stores max_set_size and lure_idx_vars as data attributes.
No method exists for imm, so dispatch passes on.
check_data.imm_full(): checks nt_distances.

General methods run first, so the model’s own method can rely on their work: by the time check_data.imm_full() runs, max_set_size is already an attribute. A new circular model with non-targets inherits all of these checks and only writes its own.

4.1.3 `configure_model()`

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 method reads what check_data() stored (Section 3.6), builds the formula with bmf2bf(), adds the non-linear terms and sets the family. It returns no stanvars, because brms::mixture() needs no Stan code from bmm.

The three IMM versions differ only in the theta expression, so what they share lives in private helpers:

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"
  )
}

Each version builds theta_exprs, one expression per non-target position, and passes it to .add_imm_mixture_terms(). This is the pattern to copy when your versions share a formula skeleton: put the skeleton in a helper prefixed with a dot, and let each version supply what differs. Do not start a helper’s name with .model_, because model_names() lists every object with that prefix as a model.

4.1.4 From `bmmformula` to `brmsformula`

bmf2bf() builds the first line of the brms formula from the response and adds one component per parameter. Most models use the default:

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
}

bmf2bf.bmmodel() runs first, because "bmmodel" is first in the class vector, and asks the rest of the chain for a first line through NextMethod(). If no model method answers, bmf2bf.default() returns NULL and bmf2bf.bmmodel() builds resp ~ 1 itself. A model that needs a different first line defines its own method, as the DDM below does. Each parameter formula then becomes brms::nlf() if one of its predictors is also a parameter and brms::lf() otherwise.

4.2 The Diffusion Decision Model (DDM)

The DDM lives in R/model_ddm.R.

4.2.1 Model definition

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 as a proportion of the boundary separation, measured from the lower boundary (0.5 = unbiased, > 0.5 = biased towards the upper response)."
  ),
  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)", sd = "exponential(1)"),
    bound = list(main = "normal(0,0.5)", effects = "normal(0,0.5)", sd = "exponential(2)"),
    ndt   = list(main = "normal(-1.5,0.5)", effects = "normal(0,0.3)", sd = "exponential(2)"),
    zr    = list(main = "normal(0,0.5)", effects = "normal(0,0.3)", sd = "exponential(2)")
  ),
  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",
      citation = glue(
        "Ratcliff, R. (1978). A theory of memory retrieval. Psychological \\
        Review, 85(2), 59-108. https://doi.org/10.1037/0033-295X.85.2.59"
      ),
      version = "NA",
      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 <- set_links(out, links)
  out
}

This is the layout use_model_template() generates: the specification in .ddm_defaults, read field by field by the constructor. The DDM differs from the IMM in three ways:

No natural mu. The family gets a dummy mu, fixed to 0 through fixed_parameters and left out of parameters (Section A.2).
Two response variables, rt and response, which need their own first formula line.
init_ranges keep the non-decision time below the fastest response at the start of sampling (Section A.3).

Its class chain is short, c("bmmodel", "ddm"), because the response time models share few data checks.

4.2.2 `bmf2bf()`

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 produces rt | dec(response) ~ 1, where dec() tells brms which column codes the decision (0 = lower, 1 = upper boundary).

4.2.3 `configure_model()`

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,
      posterior_epred = posterior_epred_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)
}

The pieces of a custom family:

dpars lists the dummy mu first, with an identity link, then the model parameters with the links from model$links, so a link the user sets reaches the family.
vars = "dec[n]" passes the decision column to the Stan likelihood. brms connects the family to Stan by name: the family "ddm" calls ddm_lpdf(), defined in inst/stan_chunks/ddm_functions.stan, which configure_model() reads and passes as stanvars.
log_lik, posterior_predict and posterior_epred are plain R functions (Chapter 6). The first two are built on dddm() and rddm() from R/distributions.R; posterior_epred_ddm() returns the expected response time.

4.1 The Interference Measurement Model (IMM)

4.1.1 Model definition

4.1.2 check_data() and the dispatch chain

4.1.3 configure_model()

4.1.4 From bmmformula to brmsformula