5 Adding a new model

If you have read Section 3 and Section 4, you should have a pretty good idea of how the bmm package functions. Now it’s time to add your new model.

The good news are, you don’t have to add any of the files manually. The bmm package includes a function use_model_template() that generates all the files with templates for all necessary functions. Thus, you can focus on the important work of filling these templates with the relevant information without worrying about missing something critical.

What `use_model_template()` does

Run it from the root of the package. It writes R/model_<model_name>.R containing the parameter specification, the model constructor, the user-facing alias, and stub methods for check_data, bmf2bf, configure_model and postprocess_brm. With custom_family = TRUE it also creates one empty inst/stan_chunks/<model_name>_<block>.stan file per entry in stanvar_blocks, and wires the code that reads them into configure_model.

The full argument list is on the reference page for use_model_template(), which is generated from the function itself and is therefore never out of date. Four things about it are worth knowing before you call it:

versions is the second positional argument. use_model_template("gcm", TRUE) means versions = TRUE, not custom_family = TRUE. Name the arguments.
The stan files are created before the R file is written. If the call fails partway through, empty .stan files are already on disk and you have to remove them by hand before you can call the function again.
An unversioned model gets version = "NA" — the literal string, not NA. That is what print() and the generated documentation will show unless you pass versions.
The template does not generate everything you may need. There are no stubs for check_formula(), configure_prior() or create_initfun(); add those yourself only if your model needs them. The generated alias also omits the @keywords bmmodel tag, which the package’s pkgdown test expects on a user-facing model alias — add it (see Chapter 6).

5.1 Example

Let’s add a new model called gcm. Let’s assume that you have tested the model in Stan and you have the Stan code ready. We want to define a custom family for the gcm model, and we want to define the following blocks: likelihood, functions (see ?brms::stanvar for an explanation of the blocks).

First you set up your system and git environment as described in Section 1. Then you can run the following code from R in the root directory of the bmm package:

use_model_template("gcm", custom_family = TRUE, stanvar_blocks = c("likelihood", "functions"))

This will create the file model_gcm.R in the R/ directory and the files gcm_likelihood.stan and gcm_functions.stan in the inst/stan_chunks/ directory. The function will also open the files in your editor, through usethis::edit_file(). You will see the following output in the console:

• Modify 'inst/stan_chunks/gcm_likelihood.stan'
• Modify 'inst/stan_chunks/gcm_functions.stan'
• Modify 'R/model_gcm.R'

Now you can fill the files with the appropriate code. The Stan files will be empty, but the R file will have the following structure:

#############################################################################!
# MODELS                                                                 ####
#############################################################################!
# see 'R/model_ddm.R' (flat defaults) or 'R/model_cswald.R' (versioned) for examples

.gcm_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)"),
    par2 = list(main = "normal(0, 0.5)", effects = "normal(0, 0.5)")
  ),
  init_ranges = list(
    par1 = c(-1, 1),
    par2 = c(0.5, 1.5)
  )
)


.model_gcm <- function(resp_var1 = NULL, required_arg1 = NULL, required_arg2 = NULL,
                                  links = NULL, call = NULL, ...) {
  out <- structure(
    list(
      resp_vars = nlist(resp_var1),
      other_vars = nlist(required_arg1, required_arg2),
      domain = "",
      task = "",
      name = "",
      citation = "",
      version = "NA",
      requirements = "",
      parameters = .gcm_defaults[["parameters"]],
      links = .gcm_defaults[["links"]],
      fixed_parameters = .gcm_defaults[["fixed_parameters"]],
      default_priors = .gcm_defaults[["priors"]],
      init_ranges = .gcm_defaults[["init_ranges"]]
    ),
    class = c("bmmodel", "gcm"),
    call = call
  )
  out$links[names(links)] <- links
  out
}
# user facing alias
# information in the title and details sections will be filled in
# automatically based on the information in the .model_gcm()$info

#' @title `r .model_gcm()$name`
#' @name gcm
#' @details `r model_info(.model_gcm())`
#' @param resp_var1 A description of the response variable
#' @param required_arg1 A description of the required argument
#' @param required_arg2 A description of the required argument
#' @param links A list of links for the model parameters.
#' @param ... used internally for testing, ignore it
#' @return An object of class `bmmodel`
#' @export
#' @examples
#' \dontrun{
#' # put a full example here (see 'R/model_ddm.R' for an example)
#' }
gcm <- function(resp_var1, required_arg1, required_arg2, links = NULL, ...) {
   call <- match.call()
   stop_missing_args()
   # uncomment if your model requires a specific backend:
   # stopif(!requireNamespace("cmdstanr", quietly = TRUE),
   #        'The "cmdstanr" package is required for this model.')
   .model_gcm(resp_var1 = resp_var1, required_arg1 = required_arg1, required_arg2 = required_arg2,
                links = links, call = call, ...)
}

#############################################################################!
# CHECK_DATA S3 methods                                                  ####
#############################################################################!
# A check_data.* function should be defined for each class of the model.
# If a model shares methods with other models, the shared methods should be
# defined in helpers-data.R. Put here only the methods that are specific to
# the model. See ?check_data for details.
# (YOU CAN DELETE THIS SECTION IF YOU DO NOT REQUIRE ADDITIONAL DATA CHECKS)

#' @export
check_data.gcm <- function(model, data, formula) {
   # retrieve required arguments
   required_arg1 <- model$other_vars$required_arg1
   required_arg2 <- model$other_vars$required_arg2

   # check the data (required)

   # compute any necessary transformations (optional)

   # save some variables as attributes of the data for later use (optional)

   NextMethod('check_data')
}

#############################################################################!
# Convert bmmformula to brmsformla methods                               ####
#############################################################################!
# A bmf2bf.* function should be defined if the default method for constructing
# the brmsformula from the bmmformula does not apply (e.g if aterms are required).
# The shared method for all `bmmodels` is defined in bmmformula.R.
# See ?bmf2bf for details.
# (YOU CAN DELETE THIS SECTION IF YOUR MODEL USES A STANDARD FORMULA WITH 1 RESPONSE VARIABLE)

#' @export
bmf2bf.gcm <- function(model, formula) {
   # retrieve required response arguments
   resp_var1 <- model$resp_vars$resp_var1
   resp_var2 <- model$resp_vars$resp_arg2

   # set the base brmsformula based
   brms_formula <- brms::bf(paste0(resp_var1, " | ", vreal(resp_var2), " ~ 1"))

   # return the brms_formula to add the remaining bmmformulas to it.
   brms_formula
}

#############################################################################!
# CONFIGURE_MODEL S3 METHODS                                             ####
#############################################################################!
# Each model should have a corresponding configure_model.* function. See
# ?configure_model for more information.

#' @export
configure_model.gcm <- function(model, data, formula) {
   # retrieve required arguments
   required_arg1 <- model$other_vars$required_arg1
   required_arg2 <- model$other_vars$required_arg2

   # retrieve arguments from the data check
   my_precomputed_var <- attr(data, 'my_precomputed_var')

   # construct brms formula from the bmm formula
   formula <- bmf2bf(model, formula)

   # construct the family & add to formula object
   gcm_family <- brms::custom_family(
     'gcm',
     dpars = c(),
     links = c(),
     lb = c(), # upper bounds for parameters
     ub = c(), # lower bounds for parameters
     type = '', # real for continous dv, int for discrete dv
     loop = TRUE, # is the likelihood vectorized
   )
   formula$family <- gcm_family

   # prepare initial stanvars to pass to brms, model formula and priors
   sc_path <- system.file('stan_chunks', package='bmm')
   stan_likelihood <- read_lines2(paste0(sc_path, '/gcm_likelihood.stan'))
   stan_functions <- read_lines2(paste0(sc_path, '/gcm_functions.stan'))

   stanvars <- stanvar(scode = stan_likelihood, block = 'likelihood') +
      stanvar(scode = stan_functions, block = 'functions')

   # return the list
   nlist(formula, data, stanvars)
}
#############################################################################!
# POSTPROCESS METHODS                                                    ####
#############################################################################!
# A postprocess_brm.* function should be defined for the model class. See
# ?postprocess_brm for details

#' @export
postprocess_brm.gcm <- function(model, fit) {
   # any required postprocessing (if none, delete this section)
   fit
}

Now you have to:

Fill the .gcm_defaults specification and the .model_gcm constructor. Most of your work is in .gcm_defaults: specify the model’s parameters (with descriptions), their link functions, which parameters are fixed (and to what value), the default priors, and the initialization ranges. It’s crucial that you set default priors for every parameter, informed by knowledge in the field — see .ddm_defaults in model_ddm.R for a worked example. In .model_gcm(), rename the response and required arguments (or delete the extras if you have none) and fill in the metadata (domain, task, name, citation, requirements). The constructor lists every field of the model object inline and reads the parameter specification from .gcm_defaults; its class must be c('bmmodel', 'gcm').

If your model has multiple versions with different parameter sets, generate the scaffold with use_model_template("gcm", versions = c("v1", "v2")). Instead of the flat .gcm_defaults, this produces a .gcm_version_table with one entry per version; a constructor that reads each field out of it (parameters = .gcm_version_table[[version]][["parameters"]], and so on) and appends paste0("gcm_", version) to the class vector; and an alias whose version argument is validated with match.arg().

The released versioned models do not all validate version this way, so do not read one of them as the norm. cswald() uses match.arg() as the template does (R/model_cswald.R:177), m3() checks the string with stopif() (R/model_m3.R:183-186), and ezdm() does not validate in the constructor at all — an unknown version reaches check_data.ezdm() and errors there (R/model_ezdm.R:200). Use the template’s match.arg() for a new model: it rejects the bad value where the user typed it.
Adjust the user-facing alias. Here you should only rename the required arguments and fill in the @examples section with a full example. Everything else will be filled in automatically based on the information in the .model_gcm function.
Fill the check_data.gcm function with the appropriate code. This function should check the data and return the data. You may or may not need to compute any transformations or save some variables as attributes of the data.
If necessary, define the bmf2bf.gcm method to convert the bmmformula to a brmsformula. You need a custom bmf2bf method if:
- Your model has multiple response variables (e.g., DDM needs rt | dec(response) ~ 1)
- Your model uses aggregate data with vreal(), vint(), trials() (e.g., EZDM)
- Your model needs special response wrapping
Keep in mind that brms automatically interprets this formula as the linear model formula for the mu parameter of your custom family. If your model does not have a natural mu parameter, suppress it: include mu = 0 in fixed_parameters, add a dummy mu (identity link) to the family’s dpars, and write the first formula line yourself (as DDM does with rt | dec(response) ~ 1). See the DDM example in Section 4 and the Suppressing mu section of the architecture chapter for details. If your model has a single response variable and no special wrapping, you can delete this section.
Fill the configure_model.gcm function with the appropriate code. This function should construct the formula, the family, and the stanvars, and return them as nlist(formula, data, stanvars). Note that configure_model does not return initial values — initialization is handled separately by create_initfun() using the init_ranges field in your model definition.

You need to fill information about your custom family, and then fill the STAN files with your STAN code. Conveniently, loading the STAN files and setting up the stanvars is set up automatically when calling the use_model_template function. Should you need to add more STAN files after you created the template, you can add the files in inst/stan_chunks/ manually and edit those lines to additionally load the manually added files.
Define init_ranges if needed. If your model uses a custom family (especially RT models), you should define init_ranges in your model definition to provide plausible starting ranges for MCMC sampling. This is particularly important for models where parameters have hard constraints relative to observed data – for example, non-decision time must be less than the fastest observed RT, otherwise the likelihood is undefined. The create_initfun() step will automatically use these to generate an initialization function. Each entry maps a parameter name to a c(lower, upper) range. If init_ranges is NULL, brms’ default initialization is used.
Fill the postprocess_brm.gcm function with the appropriate code. By post-processing, we mean changes to the fitted brms model - like renaming parameters, etc. If you don’t need any post-processing, you can delete this section.
Define log_lik and posterior_predict functions if your model uses a custom family. These R functions must be passed to brms::custom_family() and enable posterior predictive checks, model comparison (via LOO/WAIC), and posterior_predict(). They should use the d*() and r*() distribution functions from R/distributions.R. See the DDM example in Section 4.

5.2 Failure modes to expect

Five things in 1.3.2 will bite you while the model is still half-written. All of them are new or newly strict, so an older model file will not warn you about them.

configure_prior() errors if fixed_parameters names a parameter the formula never gets. Every name in fixed_parameters has to be wired into the formula by your configure_model(), as either a distributional parameter or a non-linear one. If it is neither, bmm stops with a model-definition error rather than letting a malformed b_Intercept ~ constant() prior reach brm() (#377, R/helpers-prior.R:165-174). This is the error you get when you add mu = 0 to fixed_parameters but forget to put a dummy mu in the family’s dpars.

Short parameter names are risky, and the rule is a word boundary. create_initfun() has to map Stan parameter names back to your model parameters to know which init_ranges entry applies. It matches with the regex (^|_)param(_|$) and takes the longest match when several apply (#354, #355, R/helpers-inits.R:75-87). So a parameter called s matches b_s_Intercept but not sd_id__sim, which is the point — before 1.3.2 this was a substring match and s collided with sim. Two consequences for naming: a parameter name that is a substring of another is now safe, but a name that appears as a whole underscore-delimited token somewhere else in the Stan parameter names is not. If nothing matches at all, the fallback is the first underscore-delimited token of the Stan name, which is usually wrong and silently so.

A parameter that is not a dpar is looked up in nlpars. create_initfun() resolves each parameter’s terms as bterms$dpars[[parameter]] %||% bterms$nlpars[[parameter]] (R/helpers-inits.R:56). Before #362 it only looked in dpars, and a model built as a non-linear brms formula — a native-multinomial model, for instance — failed with no applicable method for 'has_intercept' applied to an object of class "NULL". If you see that error on an older bmm, this is what it was.

A predictor that is also a parameter name and also a data column now warns. bmm decides whether a formula component is non-linear purely by whether one of its predictors is also a predicted parameter, and emits nlf() rather than lf() in that case. When the same name is also a column in the data, that choice changes the likelihood and the user has no way to see it. check_formula.bmmodel warns (#378, R/bmmformula.R:165-180). Short parameter names collide naturally with column names like accuracy or a condition coded c, which is another reason to avoid one-letter parameters in a new model.

Add your model to response_annotations(). Since #392, print.bmmodel() shows the required response variables with the coding they expect (“radians in [-pi, pi]”, “seconds”) and the default links, which answers “what does my data frame need to look like?” at the console. The annotations come from response_annotations() (R/helpers-model.R:184-205), a lookup keyed on the model’s class. A model whose class is not in that chain falls through to an empty list and prints its response variables bare. Add a branch for your class.

5.3 Testing

Unit testing is extremely important. You should test your model with the testthat package. You can use the use_test() function to create a test file for your model. See files like tests/testthat/test-model_ddm.R or tests/testthat/test-model_cswald.R for examples of how we test models. BRMS models take a long time to fit, so we don’t test the actual fitting process.

Instead we use brms’ mock backend. This is not a bmm feature — backend = "mock" and mock_fit are arguments of brms::brm(), and they reach it because bmm() passes ... through. The whole backend is two lines: whatever you pass as mock_fit is returned in place of the fitted Stan object (a function is called first). Everything before the fit — data checks, formula translation, family construction, stanvars, initial values — runs for real, which is exactly what you want to test. For example, here’s a mock test pattern:

test_that("gcm model runs without errors", {
  withr::local_options("bmm.silent" = 2)
  skip_on_cran()

  dat <- data.frame(
    resp = rnorm(10),
    exemplar = rep(c("A", "B"), 5),
    category = rep(c(1, 2), each = 5)
  )

  f <- bmmformula(param1 ~ 1, param2 ~ 1)
  mock_fit <- bmm(
    f, dat,
    gcm(resp_var1 = "resp", required_arg1 = "exemplar", required_arg2 = "category"),
    backend = "mock", mock_fit = 1, rename = FALSE
  )
  expect_equal(mock_fit$fit, 1)
})

Two things about that call are easy to get wrong. The alias generated by the template calls stop_missing_args(), so every argument without a default has to be supplied — a test that passes only the response variable errors before it reaches the pipeline. And the argument is mock_fit, not mock. You will see mock = 1 in the existing test files (48 occurrences against 21 of mock_fit = 1 in 1.3.2); it works only because R partial-matches it to mock_fit on the way through .... Write mock_fit.

The tests based on the testthat package are run every time you call the check() command. Before you submit your changes, make sure that all tests pass.

5.3.1 Parameter recovery

Additionally, you should perform a full test of the model by running it in a separate script and ensuring it gives meaningful results. At the very least, you should perform basic parameter recovery simulations for hyper-parameters (i.e., means and standard deviations) as well as subject-level parameters to give users an idea of how much data they need to adequately estimate the model. Chapter 6 covers what the study has to show — in particular why you must not simulate and fit with the same distribution function — and how to report it.

5.3.2 Distribution functions

If your model uses a custom family, you should also add distribution functions (d*, r*, and optionally p*, q*) to R/distributions.R. These are used for:

Data simulation (useful for parameter recovery and examples)
Posterior predictive checks via the log_lik and posterior_predict functions
Plotting model fits

5.4 Post-processing infrastructure

The bmm package provides several post-processing tools that work automatically for most models:

conditional_effects() — visualize marginal effects of predictors on model parameters. Works automatically via the bmmfit class.
emmeans integration — estimated marginal means via the emmeans package. Works automatically.
pp_check() — posterior predictive checks. Works automatically for standard families; custom support exists for multinomial families (like M3).

For these to work with your custom family, you need to have defined log_lik and posterior_predict functions and passed them to brms::custom_family(). No additional code is typically needed in your model file.

5.5 Add an example dataset

All new models should come with an example dataset, that can be loaded by users and can be used in the examples section. This should be either:

A new dataset that you add to the package
A dataset that already exists in the package but that can be used with the new model
A dataset that exists in another package that you can load with data() and use with the new model

For example, the vignettes for the mixture2p and mixture3p use an external dataset from the mixtur package that can be loaded with data('bays2009_full', package='mixtur'). The IMM models use a dataset included in the current package. For instructions on how to add a new dataset see here.

5.6 Add an article

All new models should come with an article that explains some basic information about the model and how to estimate it with bmm. You can use the use_article() function to create a new article. See here for more information. The articles will be published automatically on the package website under “Articles” when the pull request is approved. You can browse the source code for the existing articles in the vignettes/articles/ directory. You can see the published version of the existing vignettes here.

And that’s it! You have added a new model to the bmm package. You can now submit your changes to the bmm package repository.

What use_model_template() does