Package 'bmm'

Description

Fit computational and measurement models using full Bayesian inference. The package provides a simple and accessible interface by translating complex domain-specific models into 'brms' syntax, a powerful and flexible framework for fitting Bayesian regression models using 'Stan'. The package is designed so that users can easily apply state-of-the-art models in various research fields, and so that researchers can use it as a new model development framework. References: Frischkorn and Popov (2025) doi:10.3758/s13428-025-02643-0.

Author(s)

Maintainer: Gidon T. Frischkorn [email protected] (ORCID) [copyright holder]

Authors:

• Gidon T. Frischkorn [email protected] (ORCID) [copyright holder]

• Vencislav Popov [email protected] (ORCID) [copyright holder]

Other contributors:

• Chenyu Li [contributor]

• Paul-Christian Bürkner [email protected] (Creator of 'brms', code portions of which are used in 'bmm'.) [copyright holder]

See Also

Useful links:

• https://github.com/popov-lab/bmm

• https://venpopov.com/bmm/

• Report bugs at https://github.com/popov-lab/bmm/issues

Description

Deprecated. ezdm_summary_stats() now returns contaminant-free n_upper and n_trials, so this second step is no longer needed; applying it corrects the same cell twice. Set the guess rate with the guess_rate argument of ezdm_summary_stats() instead.

Adjusts accuracy counts (n_upper, n_trials) by removing estimated contaminant trials using binomial sampling. Contaminant trials are assumed to produce correct responses at a fixed guess rate (e.g., 0.5 for 2AFC tasks).

Usage

adjust_ezdm_accuracy(n_upper, n_trials, contaminant_prop, guess_rate = 0.5)

Arguments

Details

The function adjusts one cell at a time. Vectors of length greater than one are an error; to correct several cells, use ezdm_summary_stats() on each group.

Uses binomial sampling to estimate the number of contaminant trials and contaminant correct responses, then subtracts these from the raw counts. Because of the stochastic sampling, results will vary across calls unless a seed is set by the user.

Value

A 1-row data.frame with columns n_upper_adj and n_trials_adj (integers). When contaminant_prop is NA or <= 0, returns the original counts unchanged.

See Also

ezdm_summary_stats() for computing the summary statistics and contamination proportions

Examples

# Deprecated. ezdm_summary_stats() returns contaminant-free counts already:
set.seed(42)
rt <- c(rnorm(80, 0.55, 0.05), runif(20, 0.1, 4))
response <- c(rbinom(80, 1, 0.85), rbinom(20, 1, 0.5))
ezdm_summary_stats(rt, response, contaminant_bound = c(0.1, 4))

Aggregate long-format Remember/Know data for sdt_cdp()

Description

Reshapes long-format Remember/Know recognition data – one row per trial, or one row per response category with a count column – into the aggregated wide format sdt_cdp() is fit to: one row per cell and one count column per response category, named ⁠new1 ... new<n_new>⁠, ⁠know<k>⁠ / ⁠remember<k>⁠ (and ⁠guess<k>⁠ when the data contain "guess" judgments), with k on the unified confidence scale ⁠n_new + 1 ... n_new + n_old⁠.

Usage

aggregate_sdt_cdp_data(data, judgment, confidence, count = NULL, response = "")

Arguments

Value

A data frame with one row per cell: the cell-defining columns followed by the response count columns, ready to pass to bmm() with an sdt_cdp() model. The numbers of "new" and "old" confidence levels are inferred from the data.

Examples

dat <- data.frame(
  stimulus = c(1L, 1L, 1L, 0L, 0L, 0L),
  judgment = c("remember", "know", "new", "new", "know", "new"),
  confidence = c(3L, 2L, 1L, 1L, 3L, 1L),
  count = c(40L, 25L, 15L, 55L, 10L, 20L)
)
aggregate_sdt_cdp_data(dat, "judgment", "confidence", "count")

Description

This function applies the specified link functions in the list of links to the formula or bmmformula that is passed to it. This function is mostly used internally for configuring bmmodels.

Usage

apply_links(formula, links = nlist())

Arguments

Value

The formula or bmmformula the links have been applied to

Examples

# specify a bmmformula
form <- bmf(x ~ a + c, kappa ~ 1, a ~ 1, c ~ 1)
links <- list(a = "log", c = "logit")

apply_links(form, links)

Description

Computes the posterior area under the ROC curve (AUC). For Gaussian and Gumbel (min or max) equal-variance binary SDT the AUC is available in closed form from the d draws; otherwise it is obtained by trapezoidal integration of the model-implied curve (for rating fits, the curve swept from the posterior of d and sdratio, lifted by Ro/Rn for version = "dpsdt", so the AUC is the area under the predicted mixture ROC and not $\Phi(d/\sqrt{2})$; for version = "metad" it is the type-1 curve). The returned AUC is always the area under the full curve, not the trapezoid of the discrete operating points or the K-1 rating thresholds; for binary multi-criteria fits it is one value per curve.

Usage

auc_sdt(
  fit,
  conditions = NULL,
  probs = c(0.025, 0.975),
  criterion_points = NULL,
  ...
)

Arguments

Details

The closed form is used only when sdratio is fixed at 0, where d (which is $d_a$) equals $d'$; every fit with sdratio estimated or fixed away from 0 takes the numerical route, so the AUC is invariant to the sensitivity parameterization.

Analytical formulas (equal variance, where $d_a = d'$): normal EV-SDT $AUC = \Phi(d'/\sqrt{2})$; Gumbel-min and Gumbel-max EV-SDT $AUC = \mathrm{logistic}(g')$.

Value

A data frame of class "bmm_sdt_auc" with columns AUC, .draw, and any condition columns, plus a summary attribute (AUC_mean, AUC_lower, AUC_upper), model_class, dist and conditions attributes.

See Also

roc_sdt(), plot.bmm_sdt_auc()

Examples

## Not run: 
# Three base-rate conditions shift the criterion, which identifies sdratio
dat <- expand.grid(id = 1:20, stimulus = c(0L, 1L),
                   condition = c("liberal", "neutral", "strict"))
dat$n_trials <- 100L
# simulated as bias, measured from the midpoint between the means (0 = boundary at the midpoint)
criteria <- c(liberal = -0.5, neutral = 0, strict = 0.5)
dat$n_old <- rsdt_yn(nrow(dat), dat$n_trials, dat$stimulus, d = 1.5,
                     bias = criteria[as.character(dat$condition)],
                     sdratio = 1.3)

fit <- bmm(
  formula = bmf(d ~ 1, criterion ~ 0 + condition, sdratio ~ 1),
  data = dat,
  model = sdt_yn(response = "n_old", stimulus = "stimulus",
                 n_trials = "n_trials"),
  cores = 4,
  backend = "cmdstanr"
)

auc <- auc_sdt(fit)
auc
plot(auc)

## End(Not run)

Description

Fit a Bayesian measurement model using brms as a backend interface to Stan.

Usage

bmm(
  formula,
  data,
  model,
  prior = NULL,
  sort_data = getOption("bmm.sort_data", "check"),
  silent = getOption("bmm.silent", 1),
  backend = getOption("brms.backend", NULL),
  file = NULL,
  file_compress = TRUE,
  file_refit = getOption("bmm.file_refit", FALSE),
  ...
)

fit_model(
  formula,
  data,
  model,
  prior = NULL,
  sort_data = getOption("bmm.sort_data", "check"),
  silent = getOption("bmm.silent", 1),
  backend = getOption("brms.backend", NULL),
  ...
)

Arguments

Value

An object of class brmsfit which contains the posterior draws along with many other useful information about the model. Use methods(class = "brmsfit") for an overview on available methods.

Supported Models

The following models are supported:

Continuous reproduction

• imm(): Interference measurement model by Oberauer and Lin (2017)

• mixture2p(): Two-parameter mixture model by Zhang and Luck (2008)

• mixture3p(): Three-parameter mixture model by Bays et al (2009)

• sdm(): Signal Discrimination Model (SDM) by Oberauer (2023)

Categorical recall and n-AFC decisions

• m3(): The Multinomial / Memory Measurement Model

Processing-tree models

• mpt(): Multinomial Processing Tree (MPT) models

Detection, recognition and confidence judgments

• sdt_cdp(): Continuous Dual-Process Signal Detection Theory (CDP)

• sdt_mafc(): Signal Detection Theory (m-AFC)

• sdt_ranking(): Signal Detection Theory (Ranking)

• sdt_rating(): Signal Detection Theory (Confidence Rating)

• sdt_yn(): Signal Detection Theory (Yes/No)

Choices and response times

• cswald(): Censored-Shifted Wald Model

• ddm(): Diffusion Decision Model

• ezdm(): EZ-Diffusion Model

Type ?modelname to get information about a specific model, e.g. ?imm

bmmformula syntax

see this online article for a detailed description of the syntax and how it differs from the syntax for brmsformula

Default priors, Stan code and Stan data

For more information about the default priors in bmm and about who to extract the Stan code and data generated by bmm and brms, see the online article.

Miscellaneous

Type help(package=bmm) for a full list of available help topics.

fit_model() is a deprecated alias for bmm().

References

Frischkorn, G. T., & Popov, V. (2023). A tutorial for estimating mixture models for visual working memory tasks in brms: Introducing the Bayesian Measurement Modeling (bmm) package for R. https://doi.org/10.31234/osf.io/umt57

See Also

bmm_models(), brms::brm(), default_prior(), bmmformula(), stancode(), standata()

Examples

# generate artificial data from the Signal Discrimination Model
dat <- data.frame(y = rsdm(2000))

# define formula
ff <- bmmformula(c ~ 1, kappa ~ 1)

# fit the model
fit <- bmm(
  formula = ff,
  data = dat,
  model = sdm(resp_error = "y"),
  cores = 4,
  backend = "cmdstanr"
)

Description

Inspect whether your data is coded the way a bmmodel expects before committing to Stan compilation and sampling. The function runs the same validation pipeline that bmm() runs internally (check_data() and check_formula()), but captures all errors, warnings and messages instead of throwing them, and combines them with a readable summary of:

• the response variables, their observed values and the coding the model expects

• the data columns each parameter formula uses, and how they are coded (factor levels, character vectors, numeric variables with few unique values)

• the number of observations per design cell (crossing the grouping variables from random effects terms with the categorical predictors), flagging sparse or empty cells

• model-specific diagnostics for common data mistakes (e.g. responses in degrees rather than radians for circular models, or misplaced NA values in the non-target features of set size varying designs). These are provided by data_check_findings() methods, which new models can extend

Problems are reported as findings rather than errors, so a single call shows everything that needs fixing at once.

The report covers the data and the formula only. It stops after check_formula() and takes no prior argument, so a malformed custom prior or a problem with the initial values is out of its reach and will only surface when you call bmm().

Usage

bmm_data_check(formula, data, model, min_trials = 10)

Arguments

Value

An object of class bmm_data_check with a print method. The object is a list containing the summarized response, predictors and cells information (including the full cell count table in ⁠$cells$counts⁠), all findings, and the captured results of the validation pipeline.

See Also

bmm(), check_data(), data_check_findings()

Examples

# three-parameter mixture model: radians check, NA padding of the
# non-target features across set sizes, trials per participant x set size
bmm_data_check(
  bmf(kappa ~ 0 + set_size + (0 + set_size | ID), thetat ~ 0 + set_size),
  data = oberauer_lin_2017,
  model = mixture3p(
    resp_error = "dev_rad",
    nt_features = paste0("col_nt", 1:7),
    set_size = "set_size"
  )
)

# interference measurement model: additionally checks the NA padding
# of the non-target distances
bmm_data_check(
  bmf(c ~ 1, a ~ 1, kappa ~ 1, s ~ 1),
  data = oberauer_lin_2017,
  model = imm(
    resp_error = "dev_rad",
    nt_features = paste0("col_nt", 1:7),
    nt_distances = paste0("dist_nt", 1:7),
    set_size = "set_size",
    version = "full"
  )
)

# signal discrimination model
bmm_data_check(
  bmf(c ~ 0 + setsize, kappa ~ 0 + setsize),
  data = zhang_luck_2008,
  model = sdm(resp_error = "response_error")
)

# memory measurement model: response category counts as response variables
bmm_data_check(
  bmf(c ~ 1 + (1 | ID), a ~ 1 + (1 | ID)),
  data = oberauer_lewandowsky_2019_e1,
  model = m3(
    resp_cats = c("corr", "other", "npl"),
    num_options = c("n_corr", "n_other", "n_npl"),
    choice_rule = "simple",
    version = "ss"
  )
)

# diffusion decision model: reaction time plausibility and response coding
bmm_data_check(
  bmf(drift ~ condition + (condition | ID)),
  data = data_color_judgement_task,
  model = ddm(rt = "rt", response = "response_correct")
)

Description

Measurement models available in bmm

Usage

bmm_models()

Details

Printed, the result lists the models grouped by the task they are meant for, one line per model with its constructor and full name. The groups are: continuous reproduction; categorical recall and n-AFC decisions; detection, recognition and confidence judgments; choices and response times; processing-tree models. Type ?modelname (for example ?imm) for the arguments of a model.

Value

A character vector of model names with class bmm_models, which prints as the grouped list. Use it as a character vector in base R or dplyr, e.g. "imm" %in% bmm_models(); as.character() gives a plain one where a function refuses the class.

Examples

bmm_models()
"imm" %in% bmm_models()

Description

View or change global bmm options

Usage

bmm_options(
  sort_data,
  parallel,
  default_priors,
  silent,
  color_summary,
  file_refit,
  step_size,
  reset_options = FALSE
)

Arguments

Details

The bmm_options function is used to view or change the current bmm options. If no arguments are provided, the function will return the current options. If arguments are provided, the function will change the options and return the old options invisibly. If you provide only some of the arguments, the other options will not be changed. The options are stored in the global options list and will be used by bmm() and other functions in the bmm package. Each of these options can also be set manually using the built-in options() function, by setting the bmm.sort_data, bmm.default_priors, and bmm.silent options.

Value

A message with the current bmm options and their values, and invisibly returns the old options for use with on.exit() and friends.

Examples

# view the current options
bmm_options()

# change the options to always sort the data and to use parallel sampling
bmm_options(sort_data = TRUE, parallel = TRUE)

# restore the default options
bmm_options(reset_options = TRUE)

# you can change the options using the options() function as well
options(bmm.sort_data = TRUE, bmm.parallel = TRUE)
bmm_options()

# reset the options to their default values
bmm_options(reset_options = TRUE)

# bmm_options(sort_data = TRUE, parallel = TRUE) will also return the old options
# so you can use it with on.exit()
old_op <- bmm_options(sort_data = TRUE, parallel = TRUE)
on.exit(bmm_options(old_op))

bmm_options(reset_options = TRUE)

Description

Fitting a model with bmm() needs a C++ toolchain (a compiler and make) and a Stan backend, the software that turns the model into a program and samples from it. There are two backends: the R package cmdstanr, which runs CmdStan, a separate program installed with cmdstanr::install_cmdstan(); and the R package rstan. rstan must load with either backend, because brms uses it to store every fit. bmm_setup() checks each of these, prints whether it passed, and gives one fix for every check that failed. It does not install or change anything.

Two checks concern the toolchain. "C++ toolchain" lets R compile a small test file; "CmdStan toolchain" is cmdstanr's own check, which only looks for make and a compiler on the PATH. The second can pass while the first fails, e.g. on a Mac whose command line tools are missing, and the first is the one that decides.

With smoke_test = TRUE, it finally compiles and samples a small two-parameter mixture model (mixture2p()) through bmm(), on the backend bmm() will use. This is the only check that shows the whole chain works. On an Apple Silicon Mac it takes about 10 seconds with cmdstanr and 45 seconds with rstan, and it is skipped when a check it needs has failed. Once every check passes, rerun the bmm() call that failed.

Usage

bmm_setup(smoke_test = TRUE, backend = getOption("brms.backend", NULL))

Arguments

Value

A data frame of class bmm_setup with one row per check and the columns check, status ("pass", "fail" or "skip"), detail and fix (NA unless the check failed). It prints as a report.

See Also

bmm(), bmm_data_check()

Examples

# checks without compiling anything
bmm_setup(smoke_test = FALSE)

# compiles and samples a small model as well
bmm_setup()

Description

This function is used to specify the formulas predicting the different parameters of a bmmodel.

Usage

bmmformula(...)

bmf(...)

Arguments

Value

A list of formulas for each parameters being predicted

General formula structure

The formula argument accepts formulas of the following syntax:

parameter ~ fixed_effects + (random_effects | grouping_variable)

bmm formulas are built on brms formulas and function in nearly the same way, so you can use most of the brms formula syntax. The main differences is that in bmm formulas, the response variable is not specified in the formula. Instead, each parameter of the model is explicitly specified as the left-hand side of the formula. In brms, the response variable is always specified as the left-hand side of the first formula, which implicitly means that any predictors in the first formula are predictors of the mu parameter of the model. In general, measurement models do not all have a mu parameter, therefore it is more straigthforward to explicitely predict each parameter of the model.

For example, in the following brms formula for the drift diffusion model, the first line corresponds to the drift rate parameter, but this is not explicitely stated.

      brmsformula(rt | dec(response) ~ condition + (condition | id),
                  bs ~ 1 + (1 | id),
                  ndt ~ 1 + (1 | id),
                  bias ~ 1 + (1 | id))

In bmm, the same formula would be written as:

      bmmformula(drift ~ condition + (condition | id),
                 bs ~ 1 + (1 | id),
                 ndt ~ 1 + (1 | id),
                 bias ~ 1 + (1 | id))

and the rt and response variables would be specified in the model argument of the bmm() function.

Aside from that, the bmm formula syntax is the same as the brms formula syntax. For more information on the brms formula syntax, see brms::brmsformula().

You can also use the bmf() function as a shorthand for bmmformula().

You can also set parameters to a constant value by using par = value syntax:

      bmf(drift ~ condition + (condition | id),
          bs ~ 1 + (1 | id),
          ndt ~ 1 + (1 | id),
          bias = 0.5)

in which case the bias parameter will not be estimated but it will be fixed to 0.5

Examples

imm_formula <- bmmformula(
  c ~ 0 + set_size + (0 + set_size | id),
  a ~ 1,
  kappa ~ 0 + set_size + (0 + set_size | id)
)

# or use the shorter alias 'bmf'
imm_formula2 <- bmf(
  c ~ 0 + set_size + (0 + set_size | id),
  a ~ 1,
  kappa ~ 0 + set_size + (0 + set_size | id)
)
identical(imm_formula, imm_formula2)

Description

Binary old/new recognition data from 40 subjects, aggregated to response counts. Each subject was tested under five base-rate conditions: the proportion of old items shifts the decision criterion from conservative (br1) to liberal (br5) while leaving sensitivity unchanged, tracing a five-point binary ROC per subject. That is the cleanest design for estimating the unequal-variance ratio (sdratio) of sdt_yn(), which needs more than one operating point: a single condition with no other varying predictor yields only one hit/false-alarm pair and cannot separate a wider signal distribution from a larger d'. Counts were digitised from the frequencies reported in the original article.

Usage

broeder_schuetz_2009_e3

Format

broeder_schuetz_2009_e3

A data frame with 400 rows (40 subjects x 5 conditions x 2 stimulus types) and 5 columns:

id

Integer uniquely identifying each subject

condition

Factor with five base-rate conditions, ordered from the most conservative (br1) to the most liberal (br5) induced criterion

stimulus

Integer stimulus type: 0 = new/lure, 1 = old/target

n_old

Integer count of "old" responses in that cell: hits for old items (stimulus == 1) and false alarms for new items (stimulus == 0)

n_trials

Integer number of items presented in that cell

Source

Broeder, A., & Schuetz, J. (2009). Recognition ROCs are curvilinear—or are they? On premature arguments against the two-high-threshold model of recognition. Journal of Experimental Psychology: Learning, Memory, and Cognition, 35(3), 587–606. doi:10.1037/a0015279

Examples

## Not run: 
# Unequal-variance yes/no SDT: the criterion varies across base-rate
# conditions, while sensitivity (d) and the signal/noise SD ratio (sdratio)
# are held constant across conditions.
model <- sdt_yn(
  response = "n_old", stimulus = "stimulus", n_trials = "n_trials"
)
fit <- bmm(
  formula = bmf(
    d ~ 1 + (1 | id),
    criterion ~ 0 + condition + (1 | id),
    sdratio ~ 1
  ),
  data = broeder_schuetz_2009_e3,
  model = model,
  backend = "cmdstanr"
)

## End(Not run)

Description

Convert between parametrizations of the c parameter of the SDM distribution

Usage

c_sqrtexp2bessel(c, kappa)

c_bessel2sqrtexp(c, kappa)

Arguments

Details

c_bessel2sqrtexp converts the memory strength parameter (c) from the bessel parametrization to the sqrtexp parametrization, c_sqrtexp2bessel converts from the sqrtexp parametrization to the bessel parametrization.

See the online article for details on the parameterization. The sqrtexp parametrization is the default in the bmm package.

Value

A numeric vector of the same length as c and kappa.

Examples

c_bessel <- c_sqrtexp2bessel(c = 4, kappa = 3)
c_sqrtexp <- c_bessel2sqrtexp(c = c_bessel, kappa = 3)

Description

Given a vector of responses, and the values of non-targets, this function computes the error relative to each of the non-targets.

Usage

calc_error_relative_to_nontargets(data, response, nt_features)

Arguments

Value

A data.frame with n*m rows, where n is the number of rows of data and m is the number of non-target variables. It preserves all other columns of data, except for the non-target locations, and adds a column y_nt, which contains the transformed response error relative to the non-targets

Examples

data <- oberauer_lin_2017
data <- calc_error_relative_to_nontargets(data, "dev_rad", paste0("col_nt", 1:7))
hist(data$y_nt, breaks = 100)

Description

The helper functions deg2rad and rad2deg should add convenience in transforming data from degrees to radians and from radians to degrees.

Usage

deg2rad(deg)

rad2deg(rad)

Arguments

Value

A numeric vector of the same length as deg or rad.

Examples

degrees <- runif(100, min = 0, max = 360)
radians <- deg2rad(degrees)
degrees_again <- rad2deg(radians)

Description

Compute conditional effects for parameters of a bmmfit object. This method provides a more intuitive interface than directly calling brms::conditional_effects() on bmmfit objects, by:

• Accepting model parameter names directly (e.g., "kappa", "thetat")

• Automatically determining whether parameters are distributional or non-linear

• Optionally applying inverse link transformations to show parameters on their natural scale

Usage

## S3 method for class 'bmmfit'
conditional_effects(x, par = NULL, scale = c("native", "sampling"), ...)

Arguments

Details

Parameter Types

bmm models use two types of parameters internally:

• Non-linear parameters (nlpar): Model parameters that enter the likelihood through a non-linear formula, like kappa and thetat of mixture3p() or c and a of m3()

• Distributional parameters (dpar): Parameters of the response distribution itself, like kappa and c of sdm()

Users should not need to know this distinction - conditional_effects.bmmfit() automatically routes to the correct parameter type.

Scale Transformations

By default (scale = "native"), parameters are shown on their natural scale by applying inverse link transformations:

• log link → exp transformation

• logit link → inverse logit (probability scale)

• tan_half link → 2*atan transformation (radians)

• identity link → no transformation

Use scale = "sampling" to see parameters on the scale used during MCMC sampling.

estimate__, lower__, upper__ and the spaghetti lines follow scale. se__ does not: it is the spread of the draws on the scale brms or bmm summarises them on, so for a link other than identity it can sit on a different scale than estimate__. For example, se__ of kappa is on the sampling scale in mixture3p() and on the native scale in sdm(), at either scale. Use lower__ and upper__ to describe uncertainty on the requested scale. With robust = FALSE, estimate__ is the mean and se__ the standard deviation of the draws on that scale, and estimate__ is then mapped to the requested scale: at scale = "native" the estimate of a non-linear parameter with a log link is the exponential of the mean log value, not the mean of the exponentiated draws.

The default of robust is TRUE (median and MAD), as in brms::conditional_effects(), with one exception: thetat and thetant of mixture3p() at scale = "native" default to robust = FALSE (mean and standard deviation).

Models with a multinomial family

In models fitted with a multinomial response (m3(), sdt_rating(), sdt_cdp(), sdt_ranking()), the parameters are latent quantities: they have no predictive distribution, no response categories and no response points. For these models:

• The effect is the posterior of the parameter itself, summarised by the median and MAD of the draws; use robust = FALSE for their mean and standard deviation.

• The grid of the plotted predictor and the values of all other predictors are set by the rules of brms::conditional_effects().

• Without effects, only the effects brms lists for the parameter's own formula are returned. A parameter without predictors has no effects.

• categorical, ordinal, select_points and transform are errors, and method must be "posterior_epred" or "posterior_linpred" (or an alias brms::conditional_effects() reads as one of them, such as "fitted"), which agree for a parameter.

Value

A brms_conditional_effects object (from brms), which can be:

• Plotted directly using plot()

• Converted to a data frame for custom plotting

• Combined with other conditional effects plots

See Also

brms::conditional_effects() for the underlying brms function

Examples

## Not run: 
# Fit a mixture model with set size effect on kappa
fit <- bmm(
  formula = bmf(kappa ~ 0 + setsize, thetat ~ 1),
  data = zhang_luck_2008,
  model = mixture3p(
    resp_error = "response_error",
    nt_features = paste0("col_lure", 1:5),
    set_size = "setsize"
  )
)

# Get conditional effects for kappa on natural scale (exp of log)
ce_kappa <- conditional_effects(fit, par = "kappa")
plot(ce_kappa)

# Get conditional effects for kappa on log scale (sampling scale)
ce_kappa_log <- conditional_effects(fit, par = "kappa", scale = "sampling")
plot(ce_kappa_log)

# Get effects for thetat (memory probability)
ce_thetat <- conditional_effects(fit, par = "thetat")
plot(ce_thetat)

# Specify which effects to plot
ce_specific <- conditional_effects(fit, par = "kappa", effects = "setsize")

# Combine with other brms options
ce_detailed <- conditional_effects(
  fit,
  par = "kappa",
  effects = "setsize",
  spaghetti = TRUE,
  ndraws = 100
)

## End(Not run)

Description

This function generates the activation functions for different versions of the Memory Measurement Model (m3) implemented in the bmm package. If no bmmodel object is passed then it will print the available model versions.

Usage

construct_m3_act_funs(model = NULL, warnings = TRUE)

Arguments

Value

A bmmformula object with the activation functions for the m3 version specified in the model object. The activation functions use the names of the response categories specified in the model object.

Examples

model <- m3(
 resp_cats = c("correct","other", "npl"),
 num_options = c(1, 4, 5),
 version = "ss"
)

construct_m3_act_funs(model, warnings = FALSE)

Description

Censored-Shifted Wald Model

Usage

cswald(rt, response, links = NULL, version = c("simple", "crisk"), ...)

Arguments

Details

• Domain: Decision Making / Response times

• Task: Choice Reaction Time tasks (with few errors)

• Name: Censored-Shifted Wald Model

• Citation:

  • Miller, R., Scherbaum, S., Heck, D. W., Goschke, T., & Enge, S. (2018). On the relation between the (censored) shifted Wald and the Wiener distribution as measurement models for choice response times. Applied Psychological Measurement, 42(2), 116-135. https://doi.org/10.1177/0146621617710465

• Version: simple

• Requirements:

  • Reaction times should be passed in seconds

• The response variable should be passed numerically: 0 = lower response, 1 = upper response

• Parameters:

  • drift: drift rate

  • bound: boundary (distance from starting point to correct boundary, half the boundary separation bound of version "crisk", ddm() and ezdm())

  • ndt: non-decision time

  • s: diffusion constant

• Fixed parameters:

  • mu = 0

  • s = 0

• Default parameter links:

  • drift = log; bound = log; ndt = log; s = log

• Default priors:

  • drift:

    • main: normal(0,1)

    • effects: normal(0,0.3)

    • sd: exponential(2)

  • bound:

    • main: normal(0,0.3)

    • effects: normal(0,0.3)

    • sd: exponential(2)

  • ndt:

    • main: normal(-2,0.3)

    • effects: normal(0,0.3)

    • sd: exponential(2)

  • s:

    • main: normal(0,0.3)

    • effects: normal(0,0.2)

    • sd: exponential(2)

Value

An object of class bmmodel

See Also

dcswald() and rcswald() for the density and random generation functions.

Examples

# generate simulated data from the diffusion model
dat <- rcswald(n = 500, drift = 2, bound = 1.5, ndt = 0.3, zr = 0.5, s = 1)

# specify the model
model <- cswald(rt = "rt", response = "response", version = "simple")

# specify the formula
formula <- bmf(
  drift ~ 1,
  bound ~ 1,
  ndt ~ 1
)

# fit the model
fit <- bmm(
  formula = formula,
  data = dat,
  model = model,
  cores = 4,
  backend = "cmdstanr"
)

Description

These functions provide the density, distribution, quantile, and random generation functions for the censored shifted Wald model: cswald. The random generation (rcswald) and distribution functions (pcswald, qcswald) use rtdists::rdiffusion(), rtdists::pdiffusion(), and rtdists::qdiffusion() internally for the "crisk" version, which is theoretically consistent since the censored shifted Wald model is an approximation to the Wiener diffusion model for tasks with high accuracy (few errors).

Usage

dcswald(
  rt,
  response,
  drift,
  bound,
  ndt,
  zr = 0.5,
  s = 1,
  version = c("simple", "crisk"),
  log = TRUE
)

rcswald(n, drift, bound, ndt, zr = 0.5, s = 1)

pcswald(
  q,
  response,
  drift,
  bound,
  ndt,
  zr = 0.5,
  s = 1,
  version = "simple",
  lower.tail = TRUE,
  log.p = FALSE
)

qcswald(
  p,
  response,
  drift,
  bound,
  ndt,
  zr = 0.5,
  s = 1,
  version = "simple",
  lower.tail = TRUE,
  log.p = FALSE
)

Arguments

Details

Cumulative Distribution Function (pcswald)

For the "simple" version, the CDF is only defined for response = 1 (correct responses), as errors are treated as censored observations. The CDF returns the probability that a correct response occurs by time q. For response = 0, NA is returned with a warning.

For the "crisk" version (competing risks), the CDF computes the defective cumulative distribution P(RT <= q, response = r), which is the probability of responding with the specified response by time q. This uses

rtdists::pdiffusion() internally for accurate computation.

Quantile Function (qcswald)

The quantile function returns the response time q such that P(RT <= q) = p. Similar to the CDF, for the "simple" version this is only defined for response = 1. For the "crisk" version, this uses rtdists::qdiffusion() internally.

Value

• dcswald() returns a numeric vector of (log-)likelihoods.

• rcswald() returns a data.frame with columns rt (response times) and response (1 = upper, 0 = lower).

• pcswald() returns a numeric vector of (log-)probabilities.

• qcswald() returns a numeric vector of quantiles (response times).

See Also

rtdists::rdiffusion(), rtdists::pdiffusion(), rtdists::qdiffusion() for the underlying functions

Examples

dat <- rcswald(n = 1000, drift = 2, bound = 1, ndt = 0.3)
head(dat)
hist(dat$rt)

Description

Raw data from 50 subjects that completed a color judgement task. Participants had to detect if there were more blue or red dots in a cloud of colored dots. The task included a difficulty manipulation coded by the condition variable. In the easy condition 60 percent of the dots were either blue or red, and in the difficult condition 52 percent were either blue or red. In addition the fixation cross preceding the onset of the dot matrix had a variable presentation time coded in the fix_duration variable

Usage

data_color_judgement_task

Format

data_color_judgement_task

A data frame with 9941 rows and 10 columns:

ID

Character String uniquely identifying each subject

trial

The trial number of the experimental trials

condition

A character variable coding the difficulty condition

fix_duration

The duration in ms the fixation cross was shown before onset of the dot matrix.

target_response

The color / response that was more frequent in the dot matrix

response

The response coded in red / blue given by the participant

response_correct

Boolean variable coding if the response was correct or not.

rt

The reaction time in seconds.

pressed_key

The actual key pressed by the participants

correct_key

The correct key that should have been pressed.

Description

Diffusion Decision Model

Usage

ddm(rt, response, links = NULL, ...)

Arguments

Details

• Domain: Decision Making / Response times

• Task: Two-Alternative Force Choice RT

• Name: Diffusion Decision Model

• Citation:

  • Ratcliff, R. (1978). A theory of memory retrieval. Psychological Review, 85(2), 59-108. https://doi.org/10.1037/0033-295X.85.2.59

• Requirements:

  • 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:

  • 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).

• Fixed parameters:

  • zr = 0

  • mu = 0

• Default parameter links:

  • drift = identity; bound = log; ndt = log; zr = logit

• Default priors:

  • drift:

    • main: cauchy(0,1)

    • effects: normal(0,0.5)

    • sd: exponential(1)

  • bound:

    • main: normal(0,0.5)

    • effects: normal(0,0.5)

    • sd: exponential(2)

  • ndt:

    • main: normal(-1.5,0.5)

    • effects: normal(0,0.3)

    • sd: exponential(2)

  • zr:

    • main: normal(0,0.5)

    • effects: normal(0,0.3)

    • sd: exponential(2)

Value

An object of class bmmodel

Default behavior

By default, zr is fixed at 0 on its logit link scale, that is at an unbiased starting point of zr = 0.5. If you want to estimate zr, add a formula for zr in your bmf() call. Estimates of zr are on the logit scale too, so use plogis() to read them as proportions.

With accuracy-coded responses (1 = correct, 0 = error) the unbiased default is the usual assumption. With response-coded data (1 and 0 for the two response options, e.g. "word" and "nonword"), a bias towards one option moves the starting point; with zr fixed, that bias is absorbed by the other parameters instead. Estimate zr for response-coded data, e.g. zr ~ 1, and code drift so that its sign follows the stimulus.

Examples

# Simulate data for one subject using rddm
set.seed(123)
n_trials <- 500

# Simulate DDM data with fixed parameters
sim_data <- rddm(
  n = n_trials,
  drift = 1.5,    # drift rate
  bound = 1.2,    # boundary separation
  ndt = 0.3,      # non-decision time
  zr = 0.5        # relative starting point
)

# Prepare data frame
dat <- data.frame(
  rt = sim_data$rt,
  response = sim_data$response
)

# Define formula (intercept-only model)
ff <- bmmformula(
  drift ~ 1,
  bound ~ 1,
  ndt ~ 1
)

# Specify the DDM model
model <- ddm(rt = "rt", response = "response")

# Fit the model
fit <- bmm(
  formula = ff,
  data = dat,
  model = model,
  cores = 4,
  iter = 1000,
  backend = "cmdstanr"
)

# Check parameter recovery
summary(fit)

Description

Density and random generation function for the Diffusion Decision Model.

Usage

dddm(rt, response, drift, bound, ndt, zr = 0.5, log = TRUE)

rddm(n, drift, bound, ndt, zr = 0.5)

Arguments

Description

Obtain the default priors for a Bayesian multilevel measurement model, as well as information for which parameters priors can be specified. Given the model, the data and the formula for the model, this function will return the default priors that would be used to estimate the model. Additionally, it will return all model parameters that have no prior specified (flat priors). This can help to get an idea about which priors need to be specified and also know which priors were used if no user-specified priors were passed to the bmm() function.

The default priors in bmm tend to be more informative than the default priors in brms, as we use domain knowledge to specify the priors.

Each model parameter carries three default priors, listed in the documentation of the model: main for the intercept (or for all levels of a factor when the intercept is suppressed), effects for the remaining regression coefficients, and sd for the standard deviations of its random effects. The sd prior is applied as a blanket prior to every random-effects standard deviation of that parameter (all grouping factors, intercepts and slopes alike) and replaces the student_t(3, 0, 2.5) default of brms. To override it, address the parameter with dpar or nlpar in brms::set_prior(), e.g. set_prior("exponential(2)", class = "sd", nlpar = "kappa").

Correlations among random effects belong to a grouping factor rather than to one parameter, so they get a single default for the whole model: lkj(2) instead of the lkj(1) of brms, set whenever the model estimates a correlation matrix. To override it, use e.g. set_prior("lkj(1)", class = "cor").

Usage

## S3 method for class 'bmmformula'
default_prior(object, data, model, formula = object, ...)

Arguments

Value

A data.frame with columns specifying the prior, the class, the coef and group for each of the priors specified. Separate rows contain the information on the parameters (or parameter classes) for which priors can be specified.

See Also

bmm_models(), brms::default_prior()

Examples

default_prior(bmf(c ~ 1, kappa ~ 1),
  data = oberauer_lin_2017,
  model = sdm(resp_error = "dev_rad")
)

Description

These S3 methods allow emmeans::emmeans() to work seamlessly with bmmfit objects by automatically resolving whether a parameter is a distributional (dpar) or non-linear (nlpar) parameter in brms. Users can always use dpar = "parameter_name" regardless of internal classification.

Usage

## S3 method for class 'bmmfit'
recover_data(object, ..., dpar = NULL, nlpar = NULL)

## S3 method for class 'bmmfit'
emm_basis(object, trms, xlev, grid, ..., dpar = NULL, nlpar = NULL)

Arguments

Details

bmm models use two types of parameters internally in brms:

• Distributional parameters (dpar): Used in models with custom_family(dpars = ...) (e.g., SDM, EZDM)

• Non-linear parameters (nlpar): Used in models with bmf2bf() + nlf() (e.g., mixture2p, mixture3p, IMM, M3)

Users should not need to know this distinction. These methods intercept the dpar argument and, if the parameter is actually an nlpar, silently re-route it to nlpar.

Value

See emmeans::recover_data() and emmeans::emm_basis() for return values.

See Also

emmeans::emmeans()

Examples

## Not run: 
# Works for any bmm model — no need to know dpar vs nlpar
em <- emmeans(fit, ~ condition, dpar = "kappa")
pairs(em)
confint(em)

## End(Not run)

Description

This functions extracts the names, dimensions, and types from a compiled STAN parameters blocks generated by bmm or brms. This function is used to specify initial values for bmm models.

Usage

extract_parameter_dimensions(parameters_block)

Arguments

Value

A list of all parameters, their types, and dimensions as as specified in the STAN data generated by bmm and brms. A declaration whose type this function does not model is returned with its name and NA for its type and dimensions, so that only that parameter is left without initial values

Examples

# generate simple stan code from brms
stan_code <- stancode(brms::bf(x ~ 1 + cond + (1 + cond | ID)),
  data = data.frame(
    x = rnorm(100),
    ID = rep(1:50, each = 2),
    cond = rep(1:2, times = 50)
  )
)

extracted_program_blocks <- extract_stan_blocks(stan_code)

par_dims <- extract_parameter_dimensions(extracted_program_blocks$parameters) #'

Description

This function extracts the code from the different program blocks of a STAN program. This can be used in combination with the stancode function to access information about the STAN code generated by brms and bmm.

Usage

extract_stan_blocks(
  stan_code,
  blocks = c("functions", "data", "transformed data", "parameters",
    "transformed parameters", "model", "generated quantities")
)

Arguments

Value

A named list with each element containing the code of one of the STAN program blocks. If a block

Examples

# generate simple stan code from brms
stan_code <- stancode(brms::bf(x ~ 1), data = data.frame(x = rnorm(100)))

extracted_program_blocks <- extract_stan_blocks(stan_code)

Description

EZ-Diffusion Model

Usage

ezdm(
  mean_rt,
  var_rt,
  n_upper,
  n_trials,
  links = NULL,
  version = c("3par", "4par"),
  ...
)

Arguments

Details

• Domain: Decision Making / Response times

• Task: Choice Reaction Time tasks

• Name: EZ-Diffusion Model

• Citation:

  • Wagenmakers, E.-J., Van Der Maas, H. L. J., & Grasman, R. P. P. P. (2007). An EZ-diffusion model for response time and accuracy. Psychonomic Bulletin & Review, 14(1), 3-22. https://doi.org/10.3758/BF03194023

  • Chávez De la Peña, A. F., & Vandekerckhove, J. (2025). An EZ Bayesian hierarchical drift diffusion model for response time and accuracy. Psychonomic Bulletin & Review, 32(6), 3067-3087. https://doi.org/10.3758/s13423-025-02729-y

• Version: 4par

• Requirements:

  Provide aggregated statistics for each subject and condition that model parameters should vary over:

• Mean reaction times (mean_rt) in seconds

• Variance of reaction times (var_rt) in seconds

• Number of responses to the upper decision threshold (n_upper)

• Total number of trials used to calculate aggregated statistics (n_trials)

• Parameters:

  • 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).

  • s: The diffusion constant, that is the standard deviation of the Gaussian noise during sampling

• Fixed parameters:

  • s = 0

  • mu = 0

• Default parameter links:

  • drift = identity; bound = log; ndt = log; zr = logit; s = log

• Default priors:

  • drift:

    • main: cauchy(0,1)

    • effects: normal(0,0.5)

    • sd: exponential(1)

  • bound:

    • main: normal(0,0.5)

    • effects: normal(0,0.5)

    • sd: exponential(2)

  • ndt:

    • main: normal(-1.5,0.5)

    • effects: normal(0,0.3)

    • sd: exponential(2)

  • zr:

    • main: normal(0,0.5)

    • effects: normal(0,0.3)

    • sd: exponential(2)

  • s:

    • main: normal(0,1)

    • effects: normal(0,0.3)

    • sd: exponential(2)

  In version "4par", a boundary reached fewer than twice, or without RT summaries (NA), enters the model through the response counts only. bmm() adds two columns to the data, rt_used_upper and rt_used_lower, which it sets to 0 for such a boundary and to 1 otherwise, and replaces the summaries of such a boundary with a placeholder the likelihood never reads. A 0 already in these columns is kept, so that update() does not read the placeholders in fit$data as data.

  newdata passed to log_lik(), predict() or posterior_predict() needs these columns too. Rows of fit$data have them. For raw data, set both to 1: a boundary reached fewer than twice or with NA summaries is then still left out of the likelihood. brms functions that compare the response with predictions, predictive_error() and residuals(method = "posterior_predict"), read the placeholders as data in cells whose rt_used_upper is 0, so leave those cells out first, for example with newdata = subset(fit$data, rt_used_upper == 1).

Value

An object of class bmmodel

Examples

## Not run: 
# Minimal parameter recovery example with 3-parameter EZDM

# Simulate data from known parameters
set.seed(123)
sim_data <- rezdm(
  n = 10,
  n_trials = 100,
  drift = 2,
  bound = 1.5,
  ndt = 0.3,
  version = "3par"
)

# Add subject ID
sim_data$id <- 1:10

# Specify model
model <- ezdm(
  mean_rt = "mean_rt",
  var_rt = "var_rt",
  n_upper = "n_upper",
  n_trials = "n_trials",
  version = "3par"
)

# Specify formula with random effects
formula <- bmf(
  drift ~ 1 + (1 | id),
  bound ~ 1 + (1 | id),
  ndt ~ 1
)

# Fit model (using cmdstanr backend)
fit <- bmm(
  formula = formula,
  data = sim_data,
  model = model,
  backend = "cmdstanr",
  cores = 4,
  chains = 4,
  iter = 2000,
  warmup = 1000
)

# Check parameter recovery
summary(fit)

# Extract population-level effects
# True values: drift = 2, bound = 1.5, ndt = 0.3
# drift has an identity link; bound and ndt have log links
fixed_effects <- brms::fixef(fit)
fixed_effects["drift_Intercept", "Estimate"]
exp(fixed_effects[c("bound_Intercept", "ndt_Intercept"), "Estimate"])

## End(Not run)

Description

Density and random generation functions for the EZ-Diffusion Model. The model operates on aggregated data: mean reaction time, variance of reaction time, and number of responses to the upper boundary.

Usage

dezdm(
  mean_rt,
  var_rt,
  n_upper,
  n_trials,
  drift,
  bound,
  ndt,
  zr = 0.5,
  s = 1,
  version = c("3par", "4par"),
  log = TRUE
)

rezdm(
  n,
  n_trials,
  drift,
  bound,
  ndt,
  zr = 0.5,
  s = 1,
  version = c("3par", "4par")
)

Arguments

Details

The number of upper-boundary responses is binomial. The two RT summaries follow the joint sampling distribution of the mean and the variance of n independent decision times, matched to the exact first four cumulants of the first-passage time. Writing $\mathrm{MDT}$ and $\mathrm{VRT}$ for its mean and variance and $\kappa_3$, $\kappa_4$ for its third and fourth cumulants, and $W = \kappa_4 / n + 2\,\mathrm{VRT}^2 / (n - 1)$ for the exact variance of the sample variance,

$\mathrm{var\_rt} \sim \mathrm{Gamma}(\mathrm{VRT}^2 / W, \mathrm{VRT} / W),$

$\mathrm{mean\_rt} \mid \mathrm{var\_rt} \sim N\left(\mathrm{ndt} + \mathrm{MDT} + \frac{\kappa_3 / n}{W}(\mathrm{var\_rt} - \mathrm{VRT}),\ \sqrt{\mathrm{VRT} / n - (\kappa_3 / n)^2 / W}\right),$

so that $\mathrm{Var}(\mathrm{mean\_rt}) = \mathrm{VRT} / n$ and $\mathrm{Cov}(\mathrm{mean\_rt}, \mathrm{var\_rt}) = \kappa_3 / n$ are exact. Decision times are right-skewed (with a symmetric start point their kurtosis is 8.8 at zero drift and falls towards 3 as drift grows), so the older form that assumes normal reaction times — independent normal and scaled chi-square $\mathrm{Gamma}((n - 1)/2, (n - 1)/(2\,\mathrm{VRT}))$ terms — understates the sampling variance of var_rt by a factor of $1 + (\mathrm{kurtosis} - 3)(n - 1)/(2n)$: about 3.8 with 100 trials per cell and 3.5 to 3.6 with 10, less at accuracies above .95. With an asymmetric start point in version "4par" it is more at the boundary nearer the start point (up to 5.7 for zr between .3 and .7) and less at the other. It also ignores a correlation of about 0.7 between the two statistics, making posteriors too narrow. The terms above reduce to it when $\kappa_3 = \kappa_4 = 0$.

For version "3par" the start point is symmetric, so the decision time is independent of which boundary is hit and all n_trials responses inform one set of cumulants. For version "4par" the two boundaries have different decision-time distributions, so each is given its own summaries and its own response count. A boundary reached fewer than twice has no sample variance, and one whose summaries are NA has nothing to evaluate; dezdm() lets either contribute only through the binomial term. rezdm() and ezdm_summary_stats() code such summaries as NA, and bmm() keeps these cells: it replaces the summaries of such a boundary with a placeholder that the likelihood never reads, so the response counts still inform the fit. The per-boundary formulas condition on the realised counts, which are themselves random.

The two additional cumulants cost sampling time. In two simulated designs (30 subjects with 200 or 250 trials, 3 seeds each, one machine) a gradient of the whole model took 1.3 times as long as with the older form for version "3par" and 1.9 times for version "4par". The likelihood alone took 1.5 to 2.0 and 1.7 to 2.5 times as long, most where accuracy is near chance, so the ratio grows with the number of cells.

Simulated mean_rt is not truncated at ndt. When a summary rests on few responses (a handful of trials in version "3par", or a rarely reached boundary in version "4par", whatever n_trials is), rezdm() can return mean_rt below ndt or even mean_rt <= 0, which bmm() rejects; drop those rows before fitting.

Value

dezdm gives the log-density of the observed summary statistics under the EZDM, and rezdm generates random summary statistics from the implied sampling distributions.

References

Wagenmakers, E.-J., Van Der Maas, H. L. J., & Grasman, R. P. P. P. (2007). An EZ-diffusion model for response time and accuracy. Psychonomic Bulletin & Review, 14(1), 3-22.

Chávez De la Peña, A. F., & Vandekerckhove, J. (2025). An EZ Bayesian hierarchical drift diffusion model for response time and accuracy. Psychonomic Bulletin & Review.

Examples

# 3-parameter version (single observation)
dezdm(
  mean_rt = 0.5, var_rt = 0.02, n_upper = 80, n_trials = 100,
  drift = 2, bound = 1.5, ndt = 0.3
)

# 3-parameter version (vectorized)
dezdm(
  mean_rt = c(0.5, 0.55), var_rt = c(0.02, 0.025),
  n_upper = c(80, 75), n_trials = c(100, 100),
  drift = 2, bound = 1.5, ndt = 0.3
)

# 4-parameter version (single observation)
dezdm(
  mean_rt = c(0.45, 0.55), var_rt = c(0.018, 0.025),
  n_upper = 80, n_trials = 100,
  drift = 2, bound = 1.5, ndt = 0.3, zr = 0.55, version = "4par"
)

# generate random summary statistics
rezdm(n = 100, n_trials = 100, drift = 2, bound = 1.5, ndt = 0.3)
rezdm(
  n = 100, n_trials = 100, drift = 2, bound = 1.5, ndt = 0.3,
  zr = 0.55, version = "4par"
)

Description

Computes robust summary statistics for the EZ-Diffusion Model by fitting mixture models to raw trial-level RT data, separating contaminant responses from true responses.

Usage

ezdm_summary_stats(
  rt,
  response,
  version = c("3par", "4par"),
  distribution = c("exgaussian", "lognormal", "invgaussian"),
  method = c("mixture", "simple", "robust"),
  robust_scale = c("iqr", "mad"),
  contaminant_bound = c("min", "max"),
  min_trials = 10,
  init_contaminant = 0.05,
  max_contaminant = 0.5,
  maxit = 500,
  tol = 1e-06,
  guess_rate = 0.5
)

Arguments

Details

RT outliers and contaminant responses (fast guesses, lapses of attention) can distort the mean and variance estimates used as input to the EZ-Diffusion equations. This function addresses this by fitting a mixture model with two components: a uniform distribution for contaminants and a parametric RT distribution for true responses. Robust moments are then extracted from the fitted parametric component.

The returned counts describe the same responses as the returned moments. For version = "3par" the estimated contaminants are shared between the boundaries at guess_rate; a cell whose observed accuracy lies outside ⁠[guess_rate * p, 1 - p * (1 - guess_rate)]⁠ for an estimated proportion p cannot have arisen that way, so the contaminants are removed proportionally from both boundaries instead, leaving the observed accuracy unchanged apart from rounding, and a warning is issued. For version = "4par" each boundary is corrected by its own estimate.

This function is designed to work with dplyr::group_by() and dplyr::reframe() for grouped operations.

Value

A 1-row data.frame. For version = "3par": mean_rt, var_rt, n_upper, n_trials, contaminant_prop. For version = "4par": mean_rt_upper, mean_rt_lower, var_rt_upper, var_rt_lower, n_upper, n_trials, contaminant_prop_upper, contaminant_prop_lower.

n_upper and n_trials count the responses the reported moments rest on. Whenever a contaminant proportion was estimated, the expected number of contaminants is removed from both counts and the result is rounded; with method = "simple", method = "robust", or an EM that did not converge, they are the raw counts.

See Also

flag_contaminant_rts() for trial-level contamination probabilities, ezdm() for fitting the EZ-Diffusion Model

Examples

# Generate example data
set.seed(123)
rt <- rgamma(100, shape = 5, rate = 10) + 0.3
response <- rbinom(100, 1, 0.8)

# 3par summary stats
ezdm_summary_stats(rt, response)

# With dplyr for grouped operations
# library(dplyr)
# mydata |>
#   group_by(subject) |>
#   reframe(ezdm_summary_stats(rt, response))

# 4par version with separate upper/lower moments
ezdm_summary_stats(rt, response, version = "4par")

Description

Extract information from a brmsfit object

Usage

fit_info(fit, what)

Arguments

Value

Depends on what and the class of fit. For brmsfit objects, information about the single fit is returned. For brmsfit_list objects, a list or data.frame with the information for each fit is returned.

• "time": A data.frame with the sampling time per chain

• "time_mean": A named numeric vector with the mean sampling time

Examples

fit <- bmm(
  formula = bmmformula(c ~ 1, kappa ~ 1),
  data = data.frame(y = rsdm(1000)),
  model = sdm(resp_error = "y")
)

fit_info(fit, "time")

Description

Identifies contaminant RTs (fast guesses, attention lapses) at the trial level using mixture modeling. For each trial, it computes the posterior probability of being a contaminant given a mixture of a uniform distribution (contaminants) and an RT distribution.

The function takes a numeric vector of RTs and returns a numeric vector of contamination probabilities, making it compatible with dplyr::mutate() and dplyr::group_by() workflows.

Usage

flag_contaminant_rts(
  rt,
  distribution = c("exgaussian", "lognormal", "invgaussian"),
  contaminant_bound = c("min", "max"),
  init_contaminant = 0.05,
  max_contaminant = 0.5,
  maxit = 500,
  tol = 1e-06
)

Arguments

Details

Mixture Model

The function fits: f(RT) = pi_c * Uniform(a,b) + (1-pi_c) * f_RT(RT|theta)

where pi_c is the contaminant proportion, Uniform(a,b) is the contaminant distribution over contaminant_bound, and f_RT is the specified RT distribution with parameters theta.

Grouping

To fit separate mixtures by condition or response boundary, use dplyr::group_by() before calling this function inside dplyr::mutate().

Diagnostics

Mixture fit diagnostics (parameters, convergence, log-likelihood) are attached as the "diagnostics" attribute of the returned vector. Access them with attr(result, "diagnostics").

Value

Numeric vector of posterior contamination probabilities P(contaminant | RT), NA for every trial of a group whose EM did not converge (a threshold such as filter(contam_prob <= 0.5) drops that whole group, so check for NA first), with a "diagnostics" attribute containing a one-row data.frame with columns: mixture_params (list), contaminant_prop, converged, iterations, loglik, n_trials, distribution, method.

See Also

ezdm_summary_stats() for aggregated RT statistics with contamination handling, validate_fast_guesses() for testing whether flagged contaminants show random guessing behavior

Examples

## Not run: 
# Simulate data with contaminants
library(bmm)
set.seed(123)
rt_clean <- rgamma(150, shape = 5, rate = 10)
rt_contam <- runif(50, 0.1, 0.2)

data <- data.frame(
  rt = c(rt_clean, rt_contam),
  subject = 1,
  response = sample(c("upper", "lower"), 200, replace = TRUE)
)

# Basic usage with mutate
library(dplyr)
data <- data |>
  mutate(contam_prob = flag_contaminant_rts(rt))

# Hard threshold: remove trials with P(contaminant) > 0.5
data_clean <- data |> filter(contam_prob <= 0.5)

# Separate fits by response boundary
data <- data |>
  group_by(subject, response) |>
  mutate(contam_prob = flag_contaminant_rts(rt))

# Access diagnostics
probs <- flag_contaminant_rts(data$rt)
attr(probs, "diagnostics")

## End(Not run)

Description

Three versions of the Interference measurement model by Oberauer and Lin (2017). - 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.

Usage

imm(
  resp_error,
  nt_features,
  nt_distances,
  set_size,
  regex = FALSE,
  version = c("full", "bsc", "abc"),
  ...
)

IMMfull(resp_error, nt_features, nt_distances, set_size, regex = FALSE, ...)

IMMbsc(resp_error, nt_features, nt_distances, set_size, regex = FALSE, ...)

IMMabc(resp_error, nt_features, set_size, regex = FALSE, ...)

Arguments

Details

• Domain: Visual working memory

• Task: Continuous reproduction

• Name: Interference measurement model by Oberauer and Lin (2017).

• Citation:

  • 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: full

• Requirements:

  • 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:

  • mu1: 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

• Fixed parameters:

  • mu1 = 0

  • mu2 = 0

  • kappa2 = -100

• Default parameter links:

  • mu1 = tan_half; kappa = log; a = log; c = log; s = log

• Default priors:

  • mu1:

    • main: normal(0, 0.5)

    • effects: normal(0, 0.25)

    • sd: exponential(4)

  • kappa:

    • main: normal(2, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • a:

    • main: normal(0, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • c:

    • main: normal(0, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • s:

    • main: normal(0, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

Version: bsc

• Requirements:

  • 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:

  • mu1: 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

  • c: Context activation

  • s: Spatial similarity gradient

• Fixed parameters:

  • mu1 = 0

  • mu2 = 0

  • kappa2 = -100

• Default parameter links:

  • mu1 = tan_half; kappa = log; c = log; s = log

• Default priors:

  • mu1:

    • main: normal(0, 0.5)

    • effects: normal(0, 0.25)

    • sd: exponential(4)

  • kappa:

    • main: normal(2, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • c:

    • main: normal(0, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • s:

    • main: normal(0, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

Version: abc

• Requirements:

  • 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:

  • mu1: 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

• Fixed parameters:

  • mu1 = 0

  • mu2 = 0

  • kappa2 = -100

• Default parameter links:

  • mu1 = tan_half; kappa = log; a = log; c = log

• Default priors:

  • mu1:

    • main: normal(0, 0.5)

    • effects: normal(0, 0.25)

    • sd: exponential(4)

  • kappa:

    • main: normal(2, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • a:

    • main: normal(0, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • c:

    • main: normal(0, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

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

Value

An object of class bmmodel

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

Description

Density, distribution, and random generation functions for the interference measurement model with the location of mu, strength of cue- dependent activation c, strength of cue-independent activation a, the generalization gradient s, and the precision of memory representations kappa.

Usage

dimm(
  x,
  mu = c(0, 2, -1.5),
  dist = c(0, 0.5, 2),
  c = 5,
  a = 2,
  b = 1,
  s = 2,
  kappa = 5,
  log = FALSE
)

pimm(
  q,
  mu = c(0, 2, -1.5),
  dist = c(0, 0.5, 2),
  c = 1,
  a = 0.2,
  b = 0,
  s = 2,
  kappa = 5
)

qimm(
  p,
  mu = c(0, 2, -1.5),
  dist = c(0, 0.5, 2),
  c = 1,
  a = 0.2,
  b = 0,
  s = 2,
  kappa = 5
)

rimm(
  n,
  mu = c(0, 2, -1.5),
  dist = c(0, 0.5, 2),
  c = 1,
  a = 0.2,
  b = 1,
  s = 2,
  kappa = 5
)

Arguments

Value

dimm gives the density of the interference measurement model, pimm gives the cumulative distribution function of the interference measurement model, qimm gives the quantile function of the interference measurement model, and rimm gives the random generation function for the interference measurement model.

References

Oberauer, K., Stoneking, C., Wabersich, D., & Lin, H.-Y. (2017). Hierarchical Bayesian measurement models for continuous reproduction of visual features from working memory. Journal of Vision, 17(5), 11.

Examples

# generate random samples from the imm and overlay the density
r <- rimm(10000,
  mu = c(0, 2, -1.5), dist = c(0, 0.5, 2),
  c = 5, a = 2, s = 2, b = 1, kappa = 4
)
x <- seq(-pi, pi, length.out = 10000)
d <- dimm(x,
  mu = c(0, 2, -1.5), dist = c(0, 0.5, 2),
  c = 5, a = 2, s = 2, b = 1, kappa = 4
)
hist(r, breaks = 60, freq = FALSE)
lines(x, d, type = "l", col = "red")

Description

This function transforms the precision parameter kappa of the von Mises distribution to the circular standard deviation. Adapted from Matlab code by Paul Bays (https://www.paulbays.com/code.php)

Usage

k2sd(K)

Arguments

Value

A vector of sd values.

Examples

kappas <- runif(1000, 0.01, 100)

# calcualte SD (in radians)
SDs <- k2sd(kappas)

# transform SDs from radians to degrees
SDs_degress <- SDs * 180 / pi

# plot the relationship between kappa and circular SD
plot(kappas, SDs)
plot(kappas, SDs_degress)

Description

Computes the noise and signal latent decision-variable densities implied by a fitted sdt_yn() or sdt_rating() model, together with the response criterion (binary) or the K-1 confidence thresholds (rating). This is the canonical signal detection picture: two evidence distributions separated by the sensitivity and cut by one or more decision boundaries.

Usage

latent_sdt(
  fit,
  conditions = NULL,
  n_grid = 200,
  probs = c(0.025, 0.975),
  collapse = NULL,
  show_competitors = FALSE,
  ...
)

Arguments

Details

Densities are evaluated at the posterior-mean parameters of each condition, on the evidence axis the model's criterion is measured on: noise at 0 with unit SD and signal at sep with SD exp(sdratio), where sep = d * sqrt((1 + exp(sdratio)^2) / 2) is the separation in noise-SD units. A fit that estimated bias has its boundaries drawn on the same axis, at bias + sep / 2. d is $d_a$, measured in root-mean-square SD units, so sep equals d whenever sdratio is fixed at 0. The decision boundaries additionally carry a posterior credible band on their location. The area beyond a boundary equals the corresponding hit/false-alarm rate for every noise distribution: the likelihood evaluates the survivor function of the evidence distribution at the boundary, which is exactly that area.

Available for all four SDT models. sdt_yn() draws the criterion and sdt_rating() the K-1 confidence thresholds as boundary lines. sdt_mafc() and sdt_ranking() have no response criterion – the decision is a max/rank rule over the m alternatives – so only the densities are drawn (no boundary lines); there the noise density represents each of the m - 1 distractor alternatives and the signal density the target.

Because the densities depend only on d/sdratio, predictors that vary the criterion/thresholds only (e.g. a base-rate manipulation, as in broeder_schuetz_2009_e3) leave the densities unchanged: by default they are collapsed into a single panel with their several boundaries overlaid and colour-coded, rather than shown as repeated identical panels. Predictors that vary d/sdratio produce distinct density panels (faceted). For sdt_mafc/sdt_ranking the set size is likewise collapsed unless it predicts d. Use collapse to override the auto-detection.

Value

A data frame of class "bmm_sdt_latent" with columns x (the evidence axis), density, distribution ("noise" or "signal"), and any faceting (density) condition columns. A lines attribute holds the decision-boundary positions (columns position, lower, upper, marker, level, plus condition columns), or NULL for the criterion-free sdt_mafc/sdt_ranking models. When show_competitors = TRUE, a competitors attribute holds the max-of-distractors densities. An extra attribute carries the response-process parameters of the sdt_rating() versions as a posterior summary (parameter, mean, lower, upper, plus condition columns): for version = "dpsdt", Ro and Rn on the probability scale; for version = "metad", the M-ratio (mratio) and meta-d' (metad); NULL for version = "standard". The object also carries probs, model_class, dist, is_rating and conditions (the conditions of the density panels).

See Also

roc_sdt(), plot.bmm_sdt_latent()

Examples

## Not run: 
# Three base-rate conditions shift the criterion, which identifies sdratio
dat <- expand.grid(id = 1:20, stimulus = c(0L, 1L),
                   condition = c("liberal", "neutral", "strict"))
dat$n_trials <- 100L
# simulated as bias, measured from the midpoint between the means (0 = boundary at the midpoint)
criteria <- c(liberal = -0.5, neutral = 0, strict = 0.5)
dat$n_old <- rsdt_yn(nrow(dat), dat$n_trials, dat$stimulus, d = 1.5,
                     bias = criteria[as.character(dat$condition)],
                     sdratio = 1.3)

fit <- bmm(
  formula = bmf(d ~ 1, criterion ~ 0 + condition, sdratio ~ 1),
  data = dat,
  model = sdt_yn(response = "n_old", stimulus = "stimulus",
                 n_trials = "n_trials"),
  cores = 4,
  backend = "cmdstanr"
)

latent <- latent_sdt(fit)
latent
plot(latent)

## End(Not run)

Description

The Multinomial / Memory Measurement Model (M3) is a measurement model that was originally introduced for working memory tasks with categorical responses. It assumes that each candidate in each response category is activated by a combination of sources of activation. The probability of choosing a response category is determined by the activation of the candidates. The model can be used for any n-AFC categorical decision task.

Usage

m3(
  resp_cats,
  num_options,
  choice_rule = "softmax",
  version = c("custom", "ss", "cs"),
  ...
)

Arguments

Details

• Domain: Working Memory (categorical), Categorical Decision Making

• Task: n-AFC retrieval

• Name: The Multinomial / Memory Measurement Model

• Citation:

  • Oberauer, K., & Lewandowsky, S. (2019). Simple measurement models for complex working-memory tasks. Psychological Review, 126(6), 880-932. https://doi.org/10.1037/rev0000159

Version: ss

• Requirements:

  • Provide names for variables specifying the number of responses in a set of response categories.

  • Specify activation sources for each response categories

  • Include at least an activation source "b" for all response categories

  • Predict the specified activation at least by a fixed intercept and any additional predictors from your data

• Parameters:

  • b: Background activation. Added to each response category. Fixed for scaling, necessary in all models.

  • c: Context activation. Added to the item cued to be recalled, that is the correct item.

  • a: General activation. Added to all items that were presented during the current trial.

• Fixed parameters:

  • b = 0

• Default parameter links:

  • c = identity; a = identity

• Default priors:

  • a:

    • main: normal(3,1)

    • effects: normal(0,0.5)

    • sd: exponential(1)

  • c:

    • main: normal(3,1)

    • effects: normal(0,0.5)

    • sd: exponential(1)

Version: cs

• Requirements:

  • Provide names for variables specifying the number of responses in a set of response categories.

  • Specify activation sources for each response categories

  • Include at least an activation source "b" for all response categories

  • Predict the specified activation at least by a fixed intercept and any additional predictors from your data

• Parameters:

  • b: Background activation. Added to each response category. Fixed for scaling, necessary in all models.

  • c: Context activation. Added to the item cued to be recalled, that is the correct item.

  • a: General activation. Added to all items that were presented during the current trial.

  • f: Filtering. This parameter captures the extent to which distractors remained in working memory.

• Fixed parameters:

  • b = 0

• Default parameter links:

  • c = identity; a = identity; f = logit

• Default priors:

  • a:

    • main: normal(3,1)

    • effects: normal(0,0.5)

    • sd: exponential(1)

  • c:

    • main: normal(3,1)

    • effects: normal(0,0.5)

    • sd: exponential(1)

  • f:

    • main: logistic(0,1)

    • effects: normal(0,1)

    • sd: exponential(1)

Version: custom

• Requirements:

  • Provide names for variables specifying the number of responses in a set of response categories.

  • Specify activation sources for each response categories

  • Include at least an activation source "b" for all response categories

  • Predict the specified activation at least by a fixed intercept and any additional predictors from your data

• Parameters:

  • b: Background activation. Added to each response category. Fixed for scaling, necessary in all models.

• Fixed parameters:

  • b = 0

• Default parameter links:

  • =

• Default priors:

Missing values and reserved names

A missing response count (NA) is counted as 0. If the category has options in that row, bmm() warns and says how many counts were replaced; if it has none (num_options is 0 in that row), the 0 is true and there is no warning, as for dist in oberauer_lewandowsky_2019_e1. A missing value in a column named in num_options is an error: enter 0 where there were no options. Data columns named Y, nTrials or ⁠Idx_<category>⁠ are refused because bmm() creates columns with these names; response categories may be called Y or nTrials.

Value

An object of class bmmodel

Examples

data <- oberauer_lewandowsky_2019_e1

# initiate the model object
m3_model <- m3(
  resp_cats = c("corr", "other", "dist", "npl"),
  num_options = c("n_corr", "n_other", "n_dist", "n_npl"),
  choice_rule = "simple"
)

# specify the model formula including the activation formulas for each response category
m3_formula <- bmf(
  corr ~ b + a + c,
  other ~ b + a,
  dist ~ b + d,
  npl ~ b,
  c ~ 1 + cond + (1 + cond | ID),
  a ~ 1 + cond + (1 + cond | ID),
  d ~ 1 + (1 | ID)
)

# specify links for the model parameters
m3_model$links <- list(
  c = "log",
  a = "log",
  d = "log"
)

# check if the default priors are applied correctly
default_prior(m3_formula, data = data, model = m3_model)

# fit the model
m3_fit <- bmm(
  formula = m3_formula,
  data = data,
  model = m3_model,
  cores = 4
)

# print summary of the model
summary(m3_fit)

Description

Density and random generation functions for the memory measurement model. Please note that these functions are currently not vectorized.

Usage

dm3(x, pars, m3_model, act_funs = NULL, log = TRUE, ...)

rm3(n, size, pars, m3_model, act_funs = NULL, unpack = FALSE, ...)

Arguments

Value

dm3 gives the density of the memory measurement model, and rm3 gives the random generation function for the memory measurement model.

References

Oberauer, K., & Lewandowsky, S. (2019). Simple measurement models for complex working-memory tasks. Psychological Review, 126(6), 880–932. https://doi.org/10.1037/rev0000159

Examples

# Basic usage - b parameter is added automatically
model <- m3(
  resp_cats = c("corr", "other", "npl"),
  num_options = c(1, 4, 5),
  choice_rule = "simple",
  version = "ss"
)

# No need to provide b parameter
dm3(x = c(20, 10, 10), pars = c(a = 1, c = 2), m3_model = model)
rm3(n = 10, size = 100, pars = c(a = 1, c = 2), m3_model = model)

# Can also use full formula (activation formulas are extracted automatically)
full_formula <- bmf(
  corr ~ b + a + c,
  other ~ b + a,
  npl ~ b,
  a ~ 1,
  c ~ 1
)
rm3(
  n = 10, size = 100, pars = c(a = 1, c = 2),
  m3_model = model, act_funs = full_formula
)

## Not run: 
# Use with dplyr::reframe() for automatic unpacking into columns
library(dplyr)
library(tibble)
param_grid <- expand.grid(a = c(0.5, 1, 1.5), c = c(1, 2, 3))

simulated_data <- param_grid |>
  rowwise() |>
  reframe(
    a = a,
    c = c,
    # unpack=TRUE returns named vector; wrap in as_tibble_row for auto-unpacking
    as_tibble_row(rm3(
      n = 1, size = 100, pars = c(a = a, c = c),
      m3_model = model, unpack = TRUE
    ))
  )
# Result has columns: a, c, corr, other, npl

## End(Not run)

Description

Ranking signal-detection data from 60 subjects. On each trial, participants saw 3, 4, or 5 face images—one a studied target—and ranked them by perceived oldness; the rank assigned to the target is recorded. The set size varies across trials, so the data are aggregated to target rank-frequency counts per subject and set size, in the wide format sdt_ranking() consumes: one count column per rank position, with structural zeros where the rank exceeds the trial's set size. Fitting all set sizes jointly relies on the per-row set-size feature of sdt_ranking(): pass the set-size column to m.

Usage

meyer_grant_jakob_2025

Format

meyer_grant_jakob_2025

A data frame with 180 rows (60 subjects x set sizes 3, 4, 5) and 7 columns:

id

Factor with sequential codes s01–s60 assigned by bmm; the original participant identifiers are not shipped

set_size

Integer number of ranked items on the trial (3, 4, or 5)

rank1, rank2, rank3, rank4, rank5

Integer number of trials in which the target received that rank (rank1 = most likely target). Columns beyond set_size are structural zeros. The counts sum to 56 within each row.

Source

Meyer-Grant, C. G., & Jakob, M. (2025). Ranking tasks in recognition memory: A direct test of the two-high-threshold contrast model. Journal of Experimental Psychology: General, 154(5), 1445–1455. doi:10.1037/xge0001700. Data on OSF: https://osf.io/gtzu7/.

Examples

## Not run: 
# Ranking SDT with set size varying per row: pass the set-size column to `m`
# so trials with 3, 4, and 5 alternatives are fit jointly.
model <- sdt_ranking(
  response = c("rank1", "rank2", "rank3", "rank4", "rank5"), m = "set_size"
)
fit <- bmm(
  formula = bmf(d ~ 1 + (1 | id)),
  data = meyer_grant_jakob_2025,
  model = model,
  backend = "cmdstanr"
)

## End(Not run)

Description

Two-parameter mixture model by Zhang and Luck (2008).

Usage

mixture2p(resp_error, ...)

Arguments

Details

• Domain: Visual working memory

• Task: Continuous reproduction

• Name: Two-parameter mixture model by Zhang and Luck (2008).

• Citation:

  • Zhang, W., & Luck, S. J. (2008). Discrete fixed-resolution representations in visual working memory. Nature, 453(7192), 233-235. https://doi.org/10.1038/nature06860

• Requirements:

  • The response vairable should be in radians and represent the angular error relative to the target

• Parameters:

  • mu1: 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

  • thetat: Mixture weight for target responses

• Fixed parameters:

  • mu1 = 0

  • mu2 = 0

  • kappa2 = -100

• Default parameter links:

  • mu1 = tan_half; kappa = log; thetat = logit

• Default priors:

  • mu1:

    • main: normal(0, 0.5)

    • effects: normal(0, 0.25)

    • sd: exponential(4)

  • kappa:

    • main: normal(2, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • thetat:

    • main: logistic(0, 1)

    • effects: normal(0, 2.5)

    • sd: exponential(1)

Value

An object of class bmmodel

Examples

# generate artificial data
dat <- data.frame(y = rmixture2p(n = 2000))

# define formula
ff <- bmmformula(kappa ~ 1, thetat ~ 1)

model <- mixture2p(resp_error = "y")

# fit the model
fit <- bmm(
  formula = ff,
  data = dat,
  model = model,
  cores = 4,
  iter = 500,
  backend = "cmdstanr"
)

Description

Density, distribution, and random generation functions for the two-parameter mixture model with the location of mu, precision of memory representations kappa and probability of recalling items from memory p_mem.

Usage

dmixture2p(x, mu = 0, kappa = 5, p_mem = 0.6, log = FALSE)

pmixture2p(q, mu = 0, kappa = 7, p_mem = 0.8)

qmixture2p(p, mu = 0, kappa = 5, p_mem = 0.6)

rmixture2p(n, mu = 0, kappa = 5, p_mem = 0.6)

Arguments

Value

dmixture2p gives the density of the two-parameter mixture model, pmixture2p gives the cumulative distribution function of the two-parameter mixture model, qmixture2p gives the quantile function of the two-parameter mixture model, and rmixture2p gives the random generation function for the two-parameter mixture model.

References

Zhang, W., & Luck, S. J. (2008). Discrete fixed-resolution representations in visual working memory. Nature, 453.

Examples

# generate random samples from the mixture2p model and overlay the density
r <- rmixture2p(10000, mu = 0, kappa = 4, p_mem = 0.8)
x <- seq(-pi, pi, length.out = 10000)
d <- dmixture2p(x, mu = 0, kappa = 4, p_mem = 0.8)
hist(r, breaks = 60, freq = FALSE)
lines(x, d, type = "l", col = "red")

Description

Three-parameter mixture model by Bays et al (2009).

Usage

mixture3p(resp_error, nt_features, set_size, regex = FALSE, ...)

Arguments

Details

• Domain: Visual working memory

• Task: Continuous reproduction

• Name: Three-parameter mixture model by Bays et al (2009).

• Citation:

  • Bays, P. M., Catalao, R. F. G., & Husain, M. (2009). The precision of visual working memory is set by allocation of a shared resource. Journal of Vision, 9(10), Article 7. https://doi.org/10.1167/9.10.7

• Requirements:

  • 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:

  • mu1: 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

  • thetat: Mixture weight for target responses

  • thetant: Mixture weight for non-target responses

• Fixed parameters:

  • mu1 = 0

  • mu2 = 0

  • kappa2 = -100

• Default parameter links:

  • mu1 = tan_half; kappa = log; thetat = softmax; thetant = softmax

• Default priors:

  • mu1:

    • main: normal(0, 0.5)

    • effects: normal(0, 0.25)

    • sd: exponential(4)

  • kappa:

    • main: normal(2, 1)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • thetat:

    • main: logistic(0, 1)

    • effects: normal(0, 2.5)

    • sd: exponential(1)

  • thetant:

    • main: logistic(0, 1)

    • effects: normal(0, 2.5)

    • sd: exponential(1)

Value

An object of class bmmodel

Examples

# generate artificial data from the Bays et al (2009) 3-parameter mixture model
dat <- data.frame(
  y = rmixture3p(n = 2000, mu = c(0, 1, -1.5, 2)),
  nt1_loc = 1,
  nt2_loc = -1.5,
  nt3_loc = 2
)

# define formula
ff <- bmmformula(
  kappa ~ 1,
  thetat ~ 1,
  thetant ~ 1
)

# specify the 3-parameter model with explicit column names for non-target features
model1 <- mixture3p(resp_error = "y", nt_features = paste0("nt", 1:3, "_loc"), set_size = 4)

# fit the model
fit <- bmm(
  formula = ff,
  data = dat,
  model = model1,
  cores = 4,
  iter = 500,
  backend = "cmdstanr"
)

# alternatively specify the 3-parameter model with a regular expression to match non-target features
# this is equivalent to the previous call, but more concise
model2 <- mixture3p(resp_error = "y", nt_features = "nt.*_loc", set_size = 4, regex = TRUE)

# fit the model
fit <- bmm(
  formula = ff,
  data = dat,
  model = model2,
  cores = 4,
  iter = 500,
  backend = "cmdstanr"
)