3  BMM code structure

You do not have to edit any of the files described here to add a model: use_model_template() (Chapter 5) writes a model file, and the package finds it. This chapter explains what happens to that file when a user fits a model, so that you know which part of it runs when.

3.1 What happens when a user calls bmm()

bmm() in R/bmm.R knows nothing about specific models. It passes the user’s formula, data and model to a pipeline of S3 generics, and each model supplies methods for the steps where it differs from the default:

bmm(formula, data, model)
|- configure_fit()              R/helpers-model.R
|   |- check_model()            is the model valid, are the user's links allowed
|   |- check_data()             validate and transform the data     <- model method
|   |- check_formula()          one formula per parameter
|   |- configure_model()        brms formula, family, Stan code     <- model method
|   |- configure_prior()        default priors + the user's prior
|   `- create_initfun()         initial values for the sampler
|- brms::brm()                  compile and sample
`- postprocess_brm()            rename, store metadata              <- optional model method

configure_fit() holds the whole configuration sequence:

R/helpers-model.R

configure_fit.bmmformula <- function(formula, data = NULL, model = NULL, prior = NULL,
                                     until = "init", frame_args = list()) {
  user_formula <- formula
  # the user's expression for the data lives in the caller of the generic
  # (bmm() or an extractor). It is set before check_data(), whose methods print
  # it (order_data_query()), and again after, because some methods rebuild the
  # data frame without it
  data_name <- substitute_name(data, envir = parent.frame())
  if (!missing(data) && !is.null(data) && (is.list(data) || is.atomic(data))) {
    attr(data, "data_name") <- data_name
  }
  model <- check_model(model, data, formula)
  data <- check_data(model, data, formula)
  attr(data, "data_name") <- data_name
  formula <- check_formula(model, data, formula)
  config_args <- configure_model(model, data, formula)
  if (until == "model") {
    return(nlist(config_args, model, user_formula))
  }
  prior <- brms::do_call(
    configure_prior, c(list(model, data, config_args$formula, prior), frame_args)
  )
  # the prior decides which parameters exist, so the inits are built from it
  if (until == "init") {
    config_args$init <- brms::do_call(
      create_initfun, c(list(model, data, config_args$formula, prior), frame_args)
    )
  }
  nlist(config_args, prior, model, user_formula)
}

The same function serves stancode(), standata() and default_prior(), which stop early through until. So the methods you write for check_data() and configure_model() are also what a user sees when they inspect the Stan code or the default priors of your model without fitting it.

Three details of bmm() itself concern model authors:

• Caching. With file_refit = "on_change", a cached fit is reused only if the Stan code and data are unchanged. configure_model() must therefore generate the same Stan code for the same inputs.
• Model-frame arguments. brms_frame_args() picks out the brm() arguments that change the model frame (data2, knots, drop_unused_levels) and passes them to configure_prior() and create_initfun(). This is why those methods need ... in their signature.
• Threading. local_brms_threads() records whether the user asked for within-chain threading, because a model may need to generate different Stan code for it (Section A.5).

3.2 The bmmformula class

Users specify what to estimate with a bmmformula (alias bmf()): a list of formulas, one per model parameter. Unlike a brmsformula, it has no response variable on the left-hand side, because the response is named when the model object is constructed.

f <- bmf(kappa ~ 0 + set_size + (1 | id), c ~ 1 + (1 | id), a ~ 1)
f <- bmf(kappa ~ condition, a = 0.5)  # '=' fixes 'a' to 0.5 instead of estimating it

A parameter the user does not mention gets an intercept-only formula, with a message. During configuration, the bmf2bf() method converts the bmmformula to a brmsformula. It adds the response as the first formula line and turns each parameter formula into a linear (brms::lf()) or non-linear (brms::nlf()) component.

3.3 Models

Each model is an S3 class with an internal constructor, .model_<name>(), and a user-facing alias, <name>(). The parameters, links, fixed parameters, priors and initialization ranges live in a separate list, .<name>_defaults, or, for a model with several versions, .<name>_version_table (Section A.6). The constructor reads its fields from that list:

.my_model_defaults <- list(
  parameters = list(
    par1 = "Parameter 1 = description of parameter 1",
    par2 = "Parameter 2 = description of parameter 2"
  ),
  links = list(par1 = "identity", par2 = "log"),
  fixed_parameters = list(mu = 0),
  priors = list(
    par1 = list(main = "normal(0, 1)", effects = "normal(0, 0.5)", sd = "exponential(1)"),
    par2 = list(main = "normal(0, 0.5)", effects = "normal(0, 0.5)", sd = "exponential(2)")
  ),
  init_ranges = list(par1 = c(-1, 1), par2 = c(0.5, 1.5))
)

.model_my_model <- function(resp_var1 = NULL, required_arg1 = NULL,
                            links = NULL, call = NULL, ...) {
  out <- structure(
    list(
      resp_vars = nlist(resp_var1),
      other_vars = nlist(required_arg1),
      domain = "",
      task = "",
      name = "",
      citation = character(),
      version = "NA",
      requirements = "",
      parameters = .my_model_defaults[["parameters"]],
      links = .my_model_defaults[["links"]],
      fixed_parameters = .my_model_defaults[["fixed_parameters"]],
      default_priors = .my_model_defaults[["priors"]],
      init_ranges = .my_model_defaults[["init_ranges"]]
    ),
    class = c("bmmodel", "my_model"),
    call = call
  )
  set_links(out, links)
}

The fields of the specification:

• parameters: the parameters a user can estimate or fix, each with a description. A dummy parameter that exists only for brms stays out (Section A.2).
• links: the link function of each parameter, which maps the scale Stan samples on to the scale the parameter is interpreted on. Users can change links through the alias’s links argument; set_links() validates their choice (Section A.1).
• fixed_parameters: parameters held at a constant, such as list(mu = 0).
• priors: per parameter, a prior for the intercept (main), for regression coefficients (effects) and for random-effect standard deviations (sd); see Section 3.8.
• init_ranges: optional c(lower, upper) windows for initial values, for parameters whose likelihood is undefined in part of their range (Section A.3).

The metadata fields (domain, task, name, citation, requirements) fill the help page and the model listings. domain decides under which task group bmm_models() and the website list the model, so reuse the string of the models yours belongs with. Tests check that every constructor runs with no arguments and lists its fields in the order above.

The alias validates the user’s arguments and calls the constructor. Its help page is generated from the constructor’s fields by model_docs():

#' @title `r .model_my_model()$name`
#' @name my_model
#' @details `r model_docs(.model_my_model())`
#' @param resp_var1 A description of the response variable
#' @param required_arg1 A description of the required argument
#' @param ... used internally for testing, ignore it
#' @return An object of class `bmmodel`
#' @keywords bmmodel
#' @export
my_model <- function(resp_var1, required_arg1, links = NULL, ...) {
  call <- match.call()
  stop_missing_args()
  .model_my_model(resp_var1 = resp_var1, required_arg1 = required_arg1,
                  links = links, call = call, ...)
}

A model with versions adds version = c("v1", "v2") to the alias, validates it with match.arg() and appends the version to the class. use_model_template() generates either variant.

3.4 Class chains

The class vector of a model encodes where it sits in the taxonomy, and it decides which methods run:

class = c("bmmodel", <shared class>, <model>, <version>)

Task group	Models and class chains
Continuous reproduction	mixture2p: c("bmmodel", "circular", "mixture2p"); mixture3p: c("bmmodel", "circular", "non_targets", "mixture3p"); imm: c("bmmodel", "circular", "non_targets", "imm", "imm_full"); sdm: c("bmmodel", "circular", "sdm", "sdm_simple")
Categorical recall	m3: c("bmmodel", "m3", "m3_ss")
Processing trees	mpt: c("bmmodel", "mpt")
Signal detection	sdt_yn: c("bmmodel", "sdt", "sdt_yn"), likewise sdt_mafc, sdt_ranking, sdt_cdp; sdt_rating: c("bmmodel", "sdt", "sdt_rating", "sdt_rating_standard")
Choices and response times	ddm: c("bmmodel", "ddm"); ezdm: c("bmmodel", "ezdm", "ezdm_3par"); cswald: c("bmmodel", "cswald", "cswald_simple")

Shared classes such as "circular", "non_targets" and "sdt" carry methods that several models use: check_data.circular() checks that responses are in radians, and check_formula.sdt() refuses the stimulus column as a predictor. Dispatch walks the class vector from left to right, so check_data.bmmodel() runs first and each method hands on to the next with NextMethod(). The model’s own method runs last. Chapter 4 follows one such chain through the IMM.

3.5 S3 generics

Each pipeline step is a generic, such as

R/helpers-model.R

configure_model <- function(model, data, formula) {
   UseMethod('configure_model')
}

and a model adds a method for it, configure_model.my_model(). That is how new models are added without editing bmm(). The generics of the pipeline:

Generic	Signature	Defined in
check_model	(model, data = NULL, formula = NULL)	helpers-model.R
check_data	(model, data, formula)	helpers-data.R
check_formula	(model, data, formula)	bmmformula.R
bmf2bf	(model, formula)	bmmformula.R
configure_model	(model, data, formula)	helpers-model.R
configure_prior	(model, data, formula, user_prior, ...)	helpers-prior.R
create_initfun	(model, data, formula, prior = NULL, ...)	helpers-inits.R
postprocess_brm	(model, fit, ...)	helpers-postprocess.R

Match these signatures exactly, including ... and the = NULL defaults. R CMD check compares every method with its generic.

A model writes check_data() and configure_model() methods almost always, bmf2bf() when its response needs wrapping (DDM: rt | dec(response) ~ 1) or packs several columns (vreal(), vint(), trials()), and postprocess_brm() rarely. configure_model() returns nlist(formula, data, stanvars): the brmsformula with its family set, the data, and the Stan code for a custom family. A model built on a native brms family, such as brms::mixture() for the IMM, needs no stanvars.

Further generics let a model adjust behaviour around the pipeline. Each has a default that covers most models; Chapter 5 says when to write a method:

• before and after the fit: revert_check_data() (Section 3.7), settable_links() (Section A.1), data_check_findings();
• for post-processing: native_transform(), pp_observables() and pp_simulate().

A few more exist for single model families, such as generated_parameters() for MPT or summary_notes() for the SDT models. Their source files explain them.

3.6 Passing information from check_data() to configure_model()

check_data() and configure_model() are separate methods and cannot share local variables. When check_data() computes something that configure_model() needs, it stores it as an attribute of the data frame:

# in check_data.non_targets() (R/helpers-data.R)
attr(data, "max_set_size") <- max_set_size
attr(data, "lure_idx_vars") <- lure_idx_vars

# later, in configure_model.imm_full() (R/model_imm.R)
max_set_size <- attr(data, "max_set_size")
lure_idx <- attr(data, "lure_idx_vars")

Two cautions. The infrastructure uses the same channel for its own bookkeeping (check_data.bmmodel() sets data_name and checked), so not every attribute you find is yours. And a configure_model() method can be called directly, without check_data(), so the attribute may be missing. configure_model.sdm() recomputes the value in that case; if you do not guard, make the failure loud.

3.7 Re-running check_data() in update()

update(fit) without new data refits on the data stored in the fit. That is brms’s model frame, which holds only the variables the formula uses and is already the output of your check_data(). Running check_data() on it again fails for any model whose check adds, renames or packs columns.

So update() first calls revert_check_data(model, data) to bring the stored frame back into the shape the user’s data had. The default returns the data unchanged, which is right for a check_data() that only validates. If yours changes columns, write a method in R/update.R:

• End with NextMethod("revert_check_data"), so that a method on a shared class still runs.
• A column your check_data() needs may be missing from the stored frame, because brms keeps only the columns the formula uses. A method that has to rebuild such a column names it in the "rebuilt" attribute of the data, and update() drops it again after the check. Append to that attribute rather than overwriting it.

revert_check_data.m3() and the SDT methods are the examples. The stored_frame case your model registers for the tests (Chapter 5) is what checks that this works.

3.8 Default priors

configure_prior() reads each parameter’s entry in priors:

• main applies to the intercept. If the user drops the intercept (kappa ~ 0 + condition), it applies to each level of the factor instead. Keep it first in the list: without an intercept, it is found by position.
• effects applies to regression coefficients. Without it, brms’ defaults apply.
• sd applies to the random-effect standard deviations. Without it, brms falls back to student_t(3, 0, 2.5), which on a log link lets individual values vary by a factor of 12 around the group value. The shipped models use exponential() with a rate chosen by what the parameter measures (Section A.7).

The strings use the format of brms::set_prior(). Correlations among random effects get lkj(2) without anything in the model definition. For a parameter in both fixed_parameters and parameters, such as mu1 in the IMM, the fixed value holds unless the user gives it a formula.

3.9 File organization

File	Contents
R/bmm.R	bmm(). You do not edit it for a new model.
R/helpers-model.R	configure_fit(), check_model(), configure_model() defaults, use_model_template(), and the central tables a new model adds to: data_column_roles, response_annotations()
R/helpers-data.R, -prior.R, -inits.R, -postprocess.R, -parameters.R	The defaults of the other pipeline steps; also the RT data utilities (flag_contaminant_rts(), ezdm_summary_stats())
R/bmmformula.R	bmmformula()/bmf(), check_formula(), bmf2bf()
R/model_<name>.R	One model: specification, constructor, alias, its methods, and the log_lik, posterior_predict and posterior_epred functions of a custom family. MPT, the largest model, is split over several model_mpt_*.R files.
R/helpers-sdt.R, R/helpers-analysis.R	Code shared by the signal detection models, and their ROC/AUC analyses
R/distributions.R	The d*, r* (and where available p*, q*) functions of each model
R/update.R, R/restructure.R	update() with the revert_check_data() methods; bringing old fits up to date
R/data_check.R, R/bmm_setup.R	bmm_data_check() and bmm_setup()
post-processing files	conditional_effects.R, emmeans.R, pp_check.R, pp_observables.R, ppp.R, summary.R, plot.R, report-extract.R. They work on any fitted model.
R/brms-misc.R	A read-only copy of brms internals. Do not edit it.
inst/stan_chunks/	Stan code for custom families, named <model>_<block>.stan; code shared by several models gets its own file (sdt_dist_funs.stan)

Good model files to read first are R/model_ddm.R (flat specification), R/model_cswald.R (version table) and R/model_sdt_yn.R (aggregated counts).

When you write Stan code, remember that the sampler needs continuous gradients. Non-differentiable operations such as fabs() or fmax() cause divergent transitions; where they are mathematically needed, use a smooth stand-in such as sqrt(x^2 + epsilon) for |x|.