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:
.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)")
),
void_mu = FALSE
),
# 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.
void_mu: For models using a custom family that do not contain a location or mu parameter, for example the diffusion model, we recommend setting up a void_mu parameter. This avoids arbitrarily using one of the model parameters as the mu parameter.
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:
#' @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
warning("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:
#' @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 nt_distances
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 data = NextMethod("check_data") to call the check_data method of the more general class.
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)
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
set_size_var <- model$other_vars$set_size
nt_distances <- model$other_vars$nt_distances
# construct main brms formula from the bmm formula
formula <- bmf2bf(model, formula) +
brms::lf(kappa2 ~ 1) +
brms::lf(mu2 ~ 1) +
brms::nlf(theta1 ~ log(exp(c) + exp(a))) +
brms::nlf(kappa1 ~ kappa) +
brms::nlf(expS ~ exp(s))
# additional internal terms for the mixture model formula
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]} ~ {lure_idx[i]} * log(exp(c-expS*{nt_distances[i]}) + exp(a))",
"+ (1 - {lure_idx[i]}) * (-100)") +
glue_nlf("{mu_nts[i]} ~ {nt_features[i]}")
}
# define mixture family
formula$family <- brms::mixture(brms::von_mises("tan_half"),
brms::von_mises("identity"),
nmix = c(1, max_set_size),
order = "none")
nlist(formula, data)
}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.
# default method to paste the full brms formula for all bmmodels
#' @export
bmf2bf.bmmodel <- function(model, formula) {
# check if the model has only one response variable and extract if TRUE
brms_formula <- NextMethod("bmf2bf")
# for each dependent parameter, check if it is used as a non-linear predictor of
# another parameter and add the corresponding brms function
for (pform in formula) {
if (is_nl(pform)) {
brms_formula <- brms_formula + brms::nlf(pform)
} else {
brms_formula <- brms_formula + brms::lf(pform)
}
}
brms_formula
}
# paste first line of the brms formula for all bmmodels with 1 response variable
#' @export
bmf2bf.default <- function(model, formula){
# set base brms formula based on response
brms::bf(paste0(model$resp_vars[[1]], "~ 1"))
}The bmf2bf.bmmodel method initializes the conversion of the bmmformula to a brms_formula. The first step for this is to paste the first line of a brmsformula that includes the response as the dependent variable on the left-hand side. For models with a single response variable this is done in the bmf2bf.default method. For models with more than one response variable, you will have to provide a model specific method of bmf2bf.myModel to convert the bmmformula into the brmsformula . 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
.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)")
),
void_mu = FALSE
),
class = c('bmmodel', 'circular', 'sdm', paste0("sdm_", version)),
call = call
)
out$links[names(links)] <- links
out
}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:
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:
#' @export
configure_model.sdm <- function(model, data, formula) {
# construct the family
# note - c has a log link, but I've coded it manually for computational efficiency
sdm_simple <- brms::custom_family(
"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"))
stan_likelihood <- read_lines2(paste0(sc_path, "/sdm_simple_likelihood.stan"))
stanvars <- brms::stanvar(scode = stan_funs, block = "functions") +
brms::stanvar(scode = stan_tdata, block = "tdata") +
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
# set initial values to be sampled between [-1,1] to avoid extreme SDs that
# can cause the sampler to fail
init <- 1
# return the list
nlist(formula, data, stanvars, init)
}Lines 5-14 use the brms::custom_family function to define a custom family for the SDM model. The dpars argument specifies the parameters of the model, and the links argument specifies the link functions for the parameters. For more information, see here
Lines 17-23 read the custom STAN code from the inst/stan_chunks directory. This has to be specified with the system.file() command to ensure that the code is found when the package is installed. The stanvars object is used to pass custom STAN code to the brms package. The stanvars object is a list of brms::stanvar objects, each of which contains the STAN code for a specific part of the model. There is a separate .stan file for each part of the STAN code, and each file is read into a separate brms::stanvar object.
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:
#' @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
}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:
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)
}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.