Package {bmm}

Title:	Easy and Accessible Bayesian Measurement Models Using 'brms'
Version:	1.4.1
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>.
License:	GPL-2
URL:	https://github.com/popov-lab/bmm, https://venpopov.com/bmm/
BugReports:	https://github.com/popov-lab/bmm/issues
Depends:	R (≥ 4.1.0)
Imports:	bayesplot, brms (≥ 2.21.0), crayon, fs, glue, graphics, matrixStats, methods, mvtnorm, parallel, posterior (≥ 1.6.0), rtdists, rlang, stats, withr
Suggests:	bookdown, cmdstanr (≥ 0.7.0), cowplot, dplyr, emmeans (≥ 1.8.0), fansi, ggplot2, ggthemes, knitr, loo (≥ 2.10.0), mixtur, MPTinR, pkgbuild, remotes, rmarkdown, roxygen2, rstan, stringr, testthat (≥ 3.3.2), tidybayes, tidyr, usethis, vctrs, waldo, yaml
Additional_repositories:	https://stan-dev.r-universe.dev
Config/testthat/edition:	3
Config/testthat/parallel:	true
Config/testthat/start-first:	conditional_effects, update, model_sdt_rating, model_sdt_cdp, default-priors, model_m3, model_cswald, helpers-init, model_ezdm, pp_check
Config/Needs/website:	rmarkdown
Encoding:	UTF-8
LazyData:	true
Config/roxygen2/version:	8.1.0
NeedsCompilation:	no
Packaged:	2026-10-16 06:22:35 UTC; gidonfrischkorn
Author:	Vencislav Popov [aut, cph], Gidon T. Frischkorn [aut, cre, cph], Chenyu Li [ctb], Paul-Christian Bürkner [cph] (Creator of 'brms', code portions of which are used in 'bmm'.)
Maintainer:	Gidon T. Frischkorn <gfrischkorn@icloud.com>
Repository:	CRAN
Date/Publication:	2026-10-16 06:50:02 UTC

Easy and Accesible Bayesian Measurement Models Using 'brms'

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 gfrischkorn@icloud.com (ORCID) [copyright holder]

Authors:

• Gidon T. Frischkorn gfrischkorn@icloud.com (ORCID) [copyright holder]

• Vencislav Popov vencislav.popov@gmail.com (ORCID) [copyright holder]

Other contributors:

• Chenyu Li [contributor]

• Paul-Christian Bürkner paul.buerkner@gmail.com (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

Distribution functions for the Interference Measurement Model (IMM)

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

Distribution functions for the Signal Discrimination Model (SDM)

Description

Density, distribution function, and random generation for the Signal Discrimination Model (SDM) Distribution with location mu, memory strength c, and precision kappa. Currently only a single activation source is supported.

Usage

dsdm(x, mu = 0, c = 3, kappa = 3.5, log = FALSE, parametrization = "sqrtexp")

psdm(
  q,
  mu = 0,
  c = 3,
  kappa = 3.5,
  lower.tail = TRUE,
  log.p = FALSE,
  lower.bound = -pi,
  parametrization = "sqrtexp"
)

qsdm(p, mu = 0, c = 3, kappa = 3.5, parametrization = "sqrtexp")

rsdm(n, mu = 0, c = 3, kappa = 3.5, parametrization = "sqrtexp")

Arguments

Details

Parametrization

See the online article for details on the parameterization. Oberauer (2023) introduced the SDM with the bessel parametrization. The sqrtexp parametrization is the default in the bmm package for numerical stability and efficiency. The two parametrizations are related by the functions c_bessel2sqrtexp() and c_sqrtexp2bessel().

The cumulative distribution function

Since responses are on the circle, the cumulative distribution function requires you to choose a lower bound of integration. The default is -\pi, as for the brms::pvon_mises() function but you can choose any value in the argument lower_bound of psdm. Another useful choice is the mean of the response distribution minus \pi, e.g. lower_bound = mu-pi. This is the default in circular::pvonmises(), and it ensures that 50% of the cumulative probability mass is below the mean of the response distribution.

Value

dsdm gives the density, psdm gives the distribution function, qsdm gives the quantile function, rsdm generates random deviates, and .dsdm_integrate is a helper function for calculating the density of the SDM distribution.

References

Oberauer, K. (2023). Measurement models for visual working memory - A factorial model comparison. Psychological Review, 130(3), 841–852

Examples

# plot the density of the SDM distribution
x <- seq(-pi, pi, length.out = 10000)
plot(x, dsdm(x, 0, 2, 3),
  type = "l", xlim = c(-pi, pi), ylim = c(0, 1),
  xlab = "Angle error (radians)",
  ylab = "density",
  main = "SDM density"
)
lines(x, dsdm(x, 0, 9, 1), col = "red")
lines(x, dsdm(x, 0, 2, 8), col = "green")
legend("topright", c(
  "c=2, kappa=3.0, mu=0",
  "c=9, kappa=1.0, mu=0",
  "c=2, kappa=8, mu=1"
),
col = c("black", "red", "green"), lty = 1, cex = 0.8
)

# plot the cumulative distribution function of the SDM distribution
p <- psdm(x, mu = 0, c = 3.1, kappa = 5)
plot(x, p, type = "l")

# generate random deviates from the SDM distribution and overlay the density
r <- rsdm(10000, mu = 0, c = 3.1, kappa = 5)
d <- dsdm(x, mu = 0, c = 3.1, kappa = 5)
hist(r, breaks = 60, freq = FALSE)
lines(x, d, type = "l", col = "red")

Utility functions for Signal Detection Theory

Description

Compute sensitivity, criterion and bias from hit and false alarm rates for different SDT distribution families. A single (hit, false alarm) pair cannot identify the signal-to-noise SD ratio, so all three are the equal-variance values; see sdt_yn() for the unequal-variance model.

Usage

sdt_d(
  hit_rate,
  fa_rate,
  dist = c("normal", "gumbel_min", "gumbel_max", "logistic")
)

sdt_criterion(
  hit_rate,
  fa_rate,
  dist = c("normal", "gumbel_min", "gumbel_max", "logistic")
)

sdt_bias(
  hit_rate,
  fa_rate,
  dist = c("normal", "gumbel_min", "gumbel_max", "logistic")
)

Arguments

Value

sdt_d returns the distance between the signal and noise distributions on the latent evidence axis, obtained by inverting the decision rule "respond old when the evidence exceeds the criterion": Q(1 - FA) - Q(1 - H), where Q is the quantile function of dist. For dist = "normal" this reduces to the familiar d' = \Phi^{-1}(H) - \Phi^{-1}(FA). Because one operating point implies equal variance, this matches the d parameter of sdt_yn() whenever sdratio is at its default.

sdt_criterion returns the location of the decision criterion measured from the mean of the noise distribution, in noise-SD units, as the criterion parameter of sdt_yn() is: Q(1 - FA). It depends on the false-alarm rate alone. For dist = "normal" this is -\Phi^{-1}(FA).

sdt_bias returns the criterion measured from the midpoint between the noise and signal means instead, as the bias parameter of sdt_yn() is: (Q(1 - FA) + Q(1 - H)) / 2, which is sdt_criterion() - sdt_d() / 2. For dist = "normal" this is the familiar c = -(\Phi^{-1}(H) + \Phi^{-1}(FA)) / 2. Positive values are a conservative bias, toward "new"/"noise" responses.

References

Green, D. M., & Swets, J. A. (1966). Signal detection theory and psychophysics. Wiley.

See Also

sdt_yn(), whose d, criterion and bias parameters these functions compute in closed form from observed rates: sdt_d(), sdt_criterion() and sdt_bias() return the same quantities as the parameters of the same name whenever sdratio is at its default.

Examples

# Compute d from hit and false alarm rates (Gaussian SDT)
sdt_d(hit_rate = 0.8, fa_rate = 0.2, dist = "normal")

# The extreme-value analogue
sdt_d(hit_rate = 0.75, fa_rate = 0.25, dist = "gumbel_min")
# Compute criterion from hit and false alarm rates
sdt_criterion(hit_rate = 0.8, fa_rate = 0.2, dist = "normal")
# Compute the bias relative to the midpoint
sdt_bias(hit_rate = 0.8, fa_rate = 0.2, dist = "normal")

Adjust Accuracy Counts for Contamination (deprecated)

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

Apply link functions for parameters in a formula or bmmformula

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)

Area under the ROC curve from a fitted SDT model

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)

Convert bmmformula objects to brmsformula objects

Description

Called by configure_model() inside bmm() to convert the bmmformula into a brmsformula based on information in the model object. It will call the appropriate bmf2bf.\* methods based on the classes defined in the model_\* function.

Usage

bmf2bf(model, formula)

Arguments

Value

A brmsformula defining the response variables and the additional parameter formulas for the specified bmmodel

Examples

model <- mixture2p(resp_error = "error")

formula <- bmmformula(
  thetat ~ 0 + set_size + (0 + set_size | id),
  kappa ~ 1 + (1 | id)
)

brms_formula <- bmf2bf(model, formula)

Fit Bayesian Measurement Models

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

Human-readable pre-fit data report for bmm models

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

Measurement models available in bmm

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

View or change global bmm options

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)

Check whether this machine can fit bmm models

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

Create formula for predicting parameters of a bmmodel

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)

Recognition ROC data from Broeder & Schuetz (2009, Experiment 3)

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)

Convert between parametrizations of the c parameter of the SDM distribution

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)

Calculate response error relative to non-target values

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)

Generic S3 method for checking data based on model type

Description

Called by bmm() to automatically perform checks on the data depending on the model type. It will call the appropriate check_data methods based on the list of classes defined in the .model_* functions. For models with several classes listed, it will call the functions in the order they are listed. Thus, any operations that are common to a group of models should be defined in the appropriate check_data.* function, where \* corresponds to the shared class. For example, for the .model_imm_abc model, this corresponds to the following order of check_data.* functions: check_data() -> check_data.circular(), check_data.non_targets() the output of the final function is returned to bmm().

Usage

check_data(model, data, formula)

Arguments

Value

A data.frame with the same number of rows as the input data, but with additional columns added as necessary, any necessary transformations applied, and attributes added to the data.frame for later use. If you need to reuse variables created by the check_data.* functions in subsequent stages (e.g. in configure_model()), you can store and access them using the attr() function.

Examples

data <- oberauer_lin_2017
model <- sdmSimple(resp_error = "dev_rad")
formula <- bmf(c ~ 1, kappa ~ 1)
checked_data <- check_data(model, data, formula)

Generic S3 method for checking if the formula is valid for the specified model

Description

Generic S3 method for checking if the formula is valid for the specified model

Usage

check_formula(model, data, formula)

Arguments

Value

the formula object

Generic S3 method for checking if the model is supported and model preprocessing

Description

In addition for validating the model, specific methods might add information to the model object based on the provided data and formula

Usage

check_model(model, data = NULL, formula = NULL)

Arguments

Value

An object of type 'bmmodel'

Convert degrees to radians or radians to degrees.

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)

Conditional Effects for BMM Models

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)

Generic S3 method for configuring the model to be fit by brms

Description

Called by bmm() to automatically construct the model formula, family objects and default priors for the model specified by the user. It will call the appropriate configure_model.* functions based on the list of classes defined in the .model_* functions. Currently, we have a method only for the last class listed in the .model_* functions. This is to keep model configuration as simple as possible. In the future we may add shared methods for classes of models that share the same configuration.

Usage

configure_model(model, data, formula)

Arguments

Details

A bare bones configure_model.* method should look like this:

configure_model.newmodel <- function(model, data, formula) {

   # preprocessing - e.g. extract arguments from data check, construct new variables
   <preprocessing code>

   # construct the formula
   formula <- bmf2bf(formula, model)

   # construct the family
   family <- <code for new family>

   # construct the default prior
   prior <- <code for new prior>

   # return the list
   nlist(formula, data, family, prior)
}

Value

A named list containing at minimum the following elements:

• formula: An object of class brmsformula. The constructed model formula

• data: the user supplied data.frame, preprocessed by check_data

• family: the brms family object

• prior: the brms prior object

• stanvars: (optional) An object of class stanvars (for custom families). See brms::custom_family() for more details.

Examples

configure_model.mixture3p <- 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

  # construct initial brms formula
  formula <- bmf2bf(model, formula) +
    brms::lf(kappa2 ~ 1) +
    brms::lf(mu2 ~ 1) +
    brms::nlf(theta1 ~ thetat) +
    brms::nlf(kappa1 ~ kappa)

  # 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]} * (thetant + log(inv_ss)) + ",
        "(1 - {lure_idx[i]}) * (-100)"
      ) +
      glue_nlf("{mu_nts[i]} ~ {nt_features[i]}")
  }

  # define mixture family
  vm_list <- lapply(1:(max_set_size + 1), function(x) brms::von_mises(link = "identity"))
  vm_list$order <- "none"
  formula$family <- brms::do_call(brms::mixture, vm_list)

  nlist(formula, data)
}

Generic S3 method for configuring the default prior for a bmmodel

Description

Called by bmm() to automatically construct the priors for a given model, data and formula, and combine it with the prior given by the user. The first method executed is configure_prior.bmmodel, which will build the prior based on information from the model object such as fixed_parameters, default_priors, etc. Thus it is important to define these values in the model object. The function will also recognize if the user has specified that some parameters should be fixed to a constant and put the appropriate constant priors. Any additional priors that a developer wants to specify, which are not based on information in the model object, can be defined in the configure_prior.* method for the model. See configure_prior.imm_full for an example.

Usage

configure_prior(model, data, formula, user_prior, ...)

Arguments

Value

A brmsprior object containing the default priors for the model

Examples

configure_prior.mixture3p <- function(model, data, formula, user_prior, ...) {
  # if there is set_size 1 in the data, set constant prior over thetant for set_size1
  prior <- brms::empty_prior()
  set_size_var <- model$other_vars$set_size
  prior_cond <- any(data$ss_numeric == 1) && !is.numeric(data[[set_size_var]])

  thetant_preds <- rhs_vars(formula$pforms$thetant)
  if (prior_cond && set_size_var %in% thetant_preds) {
    prior <- prior + brms::prior_("constant(-100)",
      class = "b",
      coef = paste0(set_size_var, 1),
      nlpar = "thetant"
    )
  }
  # check if there is a random effect on theetant that include set_size as predictor
  bterms <- brms::brmsterms(formula$pforms$thetant)
  re_terms <- bterms$dpars$mu$re
  if (!is.null(re_terms)) {
    for (i in 1:nrow(re_terms)) {
      group <- re_terms$group[[i]]
      form <- re_terms$form[[i]]
      thetant_preds <- rhs_vars(form)

      if (prior_cond && set_size_var %in% thetant_preds) {
        prior <- prior + brms::prior_("constant(1e-8)",
          class = "sd",
          coef = paste0(set_size_var, 1),
          group = group,
          nlpar = "thetant"
        )
      }
    }
  }

  prior
}

Get Activation Functions for different M3 versions

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)

Generic S3 method for creating an initial values function

Description

Called by bmm() to create an initfun for models that require initial values to properly start sampling

Usage

create_initfun(model, data, formula, prior = NULL, ...)

Arguments

Value

An initfun with no arguments that generates inital values suitable for the STAN parameters generated by the respective model call

Censored-Shifted Wald Model

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

Distribution functions for the censored shifted Wald model (cswald)

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)

Construct a single finding for a pre-fit data report

Description

Builds one finding in the shape bmm_data_check() expects. Use it inside data_check_findings() methods rather than assembling the list by hand, so that an unrecognized severity is caught where it is written instead of being silently printed as a note.

Usage

data_check_finding(severity, message)

Arguments

Value

A list with elements severity and message

See Also

data_check_findings(), bmm_data_check()

Generic S3 method for model-specific pre-fit data findings

Description

Called by bmm_data_check() to collect model-specific diagnostics of common data coding mistakes. Like check_data(), methods dispatch over the model class chain from general to specific, and every method should concatenate its own findings with NextMethod(). Each finding is a list with elements severity ("warning" or "note") and message. Unlike check_data(), these diagnostics never throw - they only describe likely problems, so a report can always be produced.

Defining a method is entirely optional: the default method returns an empty list, and bmm_data_check() produces its full generic report for models without any method. Add one only when a model has a common data mistake that the hard checks in check_data() deliberately tolerate - a mistake check_data() already warns about needs no method, because bmm_data_check() captures that warning and reports it under "Hard checks".

Usage

data_check_findings(model, data, formula)

Arguments

Value

A list of findings (possibly empty), each built with data_check_finding()

See Also

data_check_finding(), bmm_data_check()

Examples

# a method for a model whose check_data() tolerates unnormalized weights
data_check_findings.my_model <- function(model, data, formula) {
  weights <- data[[model$other_vars$weights]]
  findings <- list()
  if (!isTRUE(all.equal(sum(weights), 1))) {
    findings <- list(data_check_finding(
      "warning",
      glue::glue("The weights sum to {sum(weights)} rather than 1.")
    ))
  }
  c(findings, NextMethod())
}

Example Data from a Color Judgement Task

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.

Diffusion Decision Model

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)

Distribution function for the Diffusion Decision Model (ddm)

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

Get Default priors for Measurement Models specified in BMM

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

emmeans support for bmmfit objects

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)

Extract dimension from parameters in STAN parameter block

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) #'

Extract code from different STAN program blocks

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)

EZ-Diffusion Model

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)

Distribution functions for the EZ-Diffusion Model (ezdm)

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

Compute Robust Summary Statistics for EZ-Diffusion Model

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

Extract information from a brmsfit object

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

Flag contaminant reaction times using mixture modeling

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)

Interference measurement model by Oberauer and Lin (2017).

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

Transform kappa of the von Mises distribution to the circular standard deviation

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)

Model-implied latent decision-variable distributions

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)

The Multinomial / Memory Measurement Model

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)

Category log-probability under the Gaussian choice rule of m3()

Description

R companion to the Stan function m3_gauss_logp that m3(choice_rule = "gaussian") places in each category's activation formula. brms evaluates the non-linear formula in R for log_lik(), posterior_predict() and posterior_epred(), looking the function up on the search path; it is exported for that reason and is not meant to be called directly.

Usage

m3_gauss_logp(k, ...)

Arguments

Value

The log probability of category k, with the shape of the first activation.

Distribution functions for the Memory Measurement Model (M3)

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)

Ranking recognition data from Meyer-Grant & Jakob (2025)

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)

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

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

Distribution functions for the two-parameter mixture model (mixture2p)

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

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

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

Distribution functions for the three-parameter mixture model (mixture3p)

Description

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

Usage

dmixture3p(
  x,
  mu = c(0, 2, -1.5),
  kappa = 5,
  p_mem = 0.6,
  p_nt = 0.2,
  log = FALSE
)

pmixture3p(q, mu = c(0, 2, -1.5), kappa = 5, p_mem = 0.6, p_nt = 0.2)

qmixture3p(p, mu = c(0, 2, -1.5), kappa = 5, p_mem = 0.6, p_nt = 0.2)

rmixture3p(n, mu = c(0, 2, -1.5), kappa = 5, p_mem = 0.6, p_nt = 0.2)

Arguments

Value

dmixture3p gives the density of the three-parameter mixture model, pmixture3p gives the cumulative distribution function of the two-parameter mixture model, qmixture3p gives the quantile function of the two-parameter mixture model, and rmixture3p gives the random generation function for the two-parameter mixture model.

References

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

Examples

# generate random samples from the mixture3p model and overlay the density
r <- rmixture3p(10000, mu = c(0, 2, -1.5), kappa = 4, p_mem = 0.6, p_nt = 0.2)
x <- seq(-pi, pi, length.out = 10000)
d <- dmixture3p(x, mu = c(0, 2, -1.5), kappa = 4, p_mem = 0.6, p_nt = 0.2)
hist(r, breaks = 60, freq = FALSE)
lines(x, d, type = "l", col = "red")

References for a measurement model

Description

Returns the published source(s) of a model, for the reference list of a paper that uses it.

Usage

model_citation(x, ...)

Arguments

Details

The references come from the model constructor of the installed version of bmm, so a fit made with an older version gets the current, corrected references. When the installed version no longer has a constructor for the model, or its references are empty, the copy stored in the fit or model object is returned instead.

Some versions of a model cite the papers that introduced them instead, for example sdt_rating(version = "dpsdt").

A custom m3() model is cited with the paper that introduced the M3 framework. The model structure itself is defined by the user, so its source cannot be known to bmm. The returned vector then has the attribute uncited, naming the part of the model you have to cite yourself.

bmm itself is cited separately, see citation("bmm").

Value

A character vector with one reference in APA style per element, or character(0) if no reference is available.

See Also

citation()

Examples

model_citation(sdm(resp_error = "y"))
model_citation(
  sdt_rating(response = paste0("r", 1:4), stimulus = "old", version = "dpsdt")
)

Multinomial Processing Tree (MPT) models

Description

Multinomial Processing Tree (MPT) models measure the probabilities of discrete latent cognitive states from categorical responses. A model is specified as a set of trees created with mpt_tree(), one per distinct branch structure, such as one tree for targets and one for lures in a recognition test. Branch expressions that terminate in the same response category are summed within each tree, and parameters shared across trees are equated by giving them the same name. All latent probability parameters are estimated on an unconstrained latent scale (logit or probit); predictor formulas supplied via bmmformula() apply on that latent scale.

Usage

mpt(
  trees,
  tree_id = NULL,
  covariates = NULL,
  simplex = NULL,
  restrictions = NULL,
  links = "logit",
  ...
)

Arguments

Details

• Domain: Categorical decision making, memory, and reasoning

• Task: Any task with categorical responses generated by a discrete processing-tree structure

• Name: Multinomial Processing Tree (MPT) models

• Citation:

  • Batchelder, W. H., & Riefer, D. M. (1999). Theoretical and empirical review of multinomial process tree modeling. Psychonomic Bulletin & Review, 6(1), 57-86. https://doi.org/10.3758/BF03210812

  The data hold aggregated counts: one row per participant and tree (and per condition, if conditions are crossed with the trees), with one count column per response category. The count columns must be named exactly like the branch names given to mpt_tree() (old and new in the example below); bmm() finds them by these names, so no response variable is passed. The tree_id column names the tree of each row, and the counts of a row must come from trials of that tree. Wide tables with one column per category of an EQN file, as MPTinR and TreeBUGS take them, are reshaped with mpt_data_reshape().

  A separate tree is needed only when the branch expressions differ, that is, when the trial determines which latent processes apply or which response categories are reachable (targets versus lures, a response set with or without distractors, inclusion versus exclusion instructions). Experimental factors whose effect on the parameters you want to estimate belong in the parameter formulas, exactly as in every other bmm model: bmf(D ~ 0 + condition). The two roles are independent, and one data column can serve both when a factor happens to change the structure as well. When several levels of a factor share a branch structure, add a column naming the tree and keep the factor for the formulas.

  Equating and fixing parameters. Parameters are equated across trees by giving them the same name. Beyond that, there are two mechanisms, and they differ in when the constraint is decided. A restriction (restrictions = c("Dn = Do", "g = 0.5")) is applied when the model is built: it is substituted into the branch expressions before the parameters are identified, so the restricted parameter is not part of the model, gets no formula and no prior, cannot be a simplex member, and is listed by print(model) under ⁠Restrictions:⁠ but not by parameter_info(). A constant in the formula (bmf(D ~ 1, g = 0.5)) is applied at fit time: the parameter stays part of the model, is listed by print() under ⁠Fixed:⁠ and by parameter_info() with its value, and is freed again by giving it a formula. The value is mapped to the latent scale when the constant prior is built, so default_prior() shows constant(0) for a guessing rate of 0.5 under the logit link, and summary() lists it as ⁠g_Intercept 0.00⁠ under "Constant Parameters"; prior_info() and parameter_info() show the probability, 0.5. Members of a simplex group and their praw components cannot be fixed this way; fix parameters outside the group instead. As a rule of thumb, a constant that defines the model goes into restrictions, a constant you compare against a free version goes into the formula; the section "Constraining parameters" of the MPT article walks through both.

  Default priors. Under the logit link, intercepts and regression coefficients get a logistic(0, 1) prior and random-effect SDs get exponential(1). Under the probit link, intercepts and coefficients get normal(0, 1) and SDs get exponential(2). Group-level correlation matrices get lkj(2). A parameter switched to the other link after construction (model$links$D <- "probit") takes that link's priors unless its default prior was replaced. default_prior() lists the rows a given formula produces.

  The default prior of a simplex group of K members is the uniform (symmetric Dirichlet(1)) distribution on the simplex: each member has prior mean 1/K, whatever the listing order. Under stick breaking this is a Beta(1, K - k) prior on stick k, which default_prior() shows as logitbeta(1, K - k) (or probitbeta(1, K - k)) on the intercept of the k-th praw; the last stick keeps logistic(0, 1) (or normal(0, 1)), which is the same distribution. The result is exact under both links, and under 0 + cond it holds in every cell. Under 1 + cond it holds in the reference cell only, because the effects keep their own prior. TreeBUGS instead puts a uniform prior on each branch parameter, under which the first member has prior mean 0.5 (0.50/0.25/0.25 for three members). To get it for a group gA, gB, gC, gD, pass prior = brms::set_prior("logistic(0, 1)", nlpar = c("gAraw", "gBraw")), with "normal(0, 1)" under the probit link; without coef, it covers every cell of a 0 + cond formula. A predictor on member p acts on the logit (or probit) of p's share of the probability left over by the members listed before it. A predictor on the first member therefore scales all later members by the same factor, leaving their ratios unchanged. Read such effects on the probability scale with native_parameters().

  Identifiability. Printing the model ends with an identifiability check for intercept-only formulas. It compares the number of free parameters with the degrees of freedom (response categories minus 1, summed over trees) and computes the rank of the Jacobian of the category probabilities with respect to the free parameters, from exact derivatives at five interior test values (and five test values of each covariate); parameters fixed in the formula enter at their values. A rank below the number of free parameters means that a combination of the named parameters is not identified by the branch expressions, even when the count passes. The check is local: it holds at the test values, not at the boundaries of the parameter space. A test value with a non-finite derivative is left out, and the printout says how many; if stats::D() cannot differentiate a branch expression, the rank is not computed and the count decides. The rank sees the tree parameters only: after a non-linear formula, the rank line names that formula's parameter as not analysed.

  bmm() warns when the rank is deficient and every formula is intercept-only. When a formula has a predictor or is non-linear, it only says that these may identify the parameters and that it did not check.

  Example: separate detection parameters Do and Dn in a single old/new test give three free parameters and rank 2 of 3, so one combination of Do, Dn and g cannot be estimated; a second condition that varies the response bias identifies them.

  Reading the estimates. summary() reports intercepts and regression coefficients on the latent (logit or probit) scale. native_parameters() returns the posterior draws of every parameter on the probability scale for each observed combination of the predictors, e.g. native_parameters(fit, re_formula = NA, summary = TRUE); a difference between conditions is the difference of these draws. With re_formula = NA these are the probabilities of the median participant (the quantity TreeBUGS reports as "mean"), not the population mean of the participants' probabilities; see the section on group-level effects in native_parameters().

  Person effects written as (1 | id) in each parameter's formula are independent across parameters. A shared label, (1 |p| id) in every formula, estimates their correlations (prior lkj(2)), which corresponds to the latent-trait MPT model.

  For posterior predictive checks, brms::posterior_predict() returns simulated counts (draws x rows x categories) and brms::posterior_epred() the expected counts; ppp() computes the T1 and T2 statistics of Klauer (2010) and their posterior predictive p-values from them. With newdata, brms::posterior_predict() needs the columns nTrials, ⁠Idx_<tree>⁠ and ⁠Poss_<category>⁠ that check_data() builds; fit$data contains them.

  Non-linear formulas. Parameters that receive a non-linear predictor formula (a formula whose right-hand side references other formula parameters, e.g. D ~ Dmax * (1 - exp(-rate * ptime))) are not transformed by the link function: the user-supplied expression must keep the parameter within (0, 1). Sub-parameters of such formulas (e.g., Dmax and rate) are estimated on the identity scale with normal(0, 1) default priors on their intercepts and normal(0, 0.5) on their regression coefficients. Their random-effect SDs keep brms's default prior, student_t(3, 0, 2.5), because no single rate fits every scale such a formula can give them.

  Order constraints between parameters (Do > Dn) are expressed by reparameterizing the larger parameter in the model formula, e.g. bmf(Do ~ Dn + (1 - Dn) * inv_logit(phi), Dn ~ 1, phi ~ 1); the section "Order constraints" of the MPT article walks through the recipe.

  Names. Parameter and response category names must start with a letter and may contain only letters and digits, because brms does not allow underscores or dots in non-linear parameter names.

  All latent parameters are non-linear parameters of the underlying brms model, and brms has no separate Intercept prior class for those: their intercept prior is stored as ⁠class = "b", coef = "Intercept"⁠. A prior given as ⁠class = "b", nlpar = "D"⁠ replaces the defaults of all coefficients of D, the intercept included; add coef = "Intercept" to set the intercept alone. Call default_prior() to see which rows a given formula produces.

Value

An object of class bmmodel

Examples

## Not run: 
# two-high-threshold (2HTM) model of recognition memory
# one tree per probe type: targets (studied items) and lures (new items),
# each ending in the responses "old" and "new"
tree_target <- mpt_tree("target",
  old ~ D + (1 - D) * g,
  new ~ (1 - D) * (1 - g)
)
tree_lure <- mpt_tree("lure",
  old ~ (1 - D) * g,
  new ~ D + (1 - D) * (1 - g)
)

model <- mpt(
  trees = list(tree_target, tree_lure),
  tree_id = "item_type"
)

# simulate data for 20 participants with D = 0.7, g = 0.5:
# one row per participant x tree; count columns named after the branches
data <- data.frame(
  id = rep(1:20, each = 2),
  item_type = rep(c("target", "lure"), 20)
)
counts <- t(sapply(data$item_type, function(cond) {
  rmpt(
    n = 1, size = 50, pars = c(D = 0.7, g = 0.5),
    mpt_model = model, tree = cond, unpack = TRUE
  )
}))
data <- cbind(data, counts)

# predict both parameters by a fixed intercept and correlated random
# intercepts (the shared label |p| estimates the correlation)
formula <- bmf(
  D ~ 1 + (1 |p| id),
  g ~ 1 + (1 |p| id)
)

fit <- bmm(
  formula = formula,
  data = data,
  model = model,
  cores = 4
)

summary(fit)
# population-level D and g on the probability scale
native_parameters(fit, re_formula = NA, summary = TRUE)

## End(Not run)

Reshape wide MPT frequency tables to the long format of bmm()

Description

MPTinR and TreeBUGS take one row per participant and one column per response category of the EQN file. bmm() takes one row per participant and tree, a column naming the tree, and one count column per shared response category (see mpt()). mpt_data_reshape() reshapes the first layout into the second for a given model.

Usage

mpt_data_reshape(data, model, id = NULL, columns = c("names", "mptinr"))

Arguments

Details

Named columns (columns = "names"). For a model from mpt_from_eqn(), a column is matched to a category by the EQN label (hit), or, when trees share labels, by the tree name and label (target.hit, with the bmm name of the tree); the qualified name wins where both exist. A label used in several trees must be qualified. Labels that the categories argument of mpt_from_eqn() merged onto one response category are summed. All columns of data that are neither count columns nor id are carried unchanged to every long row, e.g. between-participant covariates.

Positional columns (columns = "mptinr"). MPTinR data for an EQN model have no column names. Every column of data except id is read as a count, in the order of check.mpt() in MPTinR (version 1.14.1), elements eqn.order.trees and eqn.order.categories: trees sorted by their label in the EQN file and, within each tree, the category labels of the file sorted. No other column can be carried; add covariates afterwards, e.g. with merge() on the participant column. The help page of check.mpt() states only that these elements give the order; the sorting rule is that of its internal .get.EQN.model.order(), which sorts the columns that read.table() returns for the file: numerically if all tree labels (all category labels) are numbers, else as character strings in the collation of the R session. mpt_data_reshape() applies the same rule, so run it in the session where MPTinR's order was determined if the labels mix upper and lower case. This mode needs a model from mpt_from_eqn(): MPTinR's own model format orders the columns by the lines of the model text, and mpt() reorders the branches of every tree to the categories of the first, so the line order cannot be recovered from the model. For other models, name the columns ⁠<tree>.<category>⁠ and use columns = "names".

Impossible categories. A response category that a tree cannot produce (impossible in mpt_tree(), mpt_from_string() and mpt_from_eqn()) has no column in data and gets 0 in the long rows of that tree, which bmm() accepts.

Counts. The count columns must be numeric, without missing values, and hold non-negative whole numbers; the long data hold them as integers.

Checks. The function cannot tell a correct reshape from one that swaps columns of equal meaning. Run bmm_data_check() on the result, and compare one participant's long rows with the wide row (see the Examples of mpt_from_eqn()).

Value

A data frame with one row per participant and tree, participants first. Its columns are the participant column (id), the tree column named by the tree_id of the model (absent in a model without it), the carried columns, and one integer count column per response category of the model.

See Also

mpt_from_eqn() and mpt_from_string() for importing the model itself; mpt() for the data layout bmm() expects; bmm_data_check() to check the result before fitting.

Examples

eqn_file <- tempfile(fileext = ".eqn")
writeLines(c(
  "6",
  "target  hit   D_o",
  "target  hit   (1-D_o)*g",
  "target  miss  (1-D_o)*(1-g)",
  "lure    fa    (1-D_n)*g",
  "lure    cr    D_n",
  "lure    cr    (1-D_n)*(1-g)"
), eqn_file)
model <- mpt_from_eqn(
  eqn_file,
  restrictions = "D_n = D_o",
  categories = c(hit = "yes", fa = "yes", miss = "no", cr = "no"),
  tree_id = "item_type"
)

# TreeBUGS data: one named column per category label of the EQN file.
# Columns that are no counts, here a between-participant factor, are
# carried to the long data
wide <- data.frame(
  subject = c("s1", "s2"), group = c("a", "b"),
  hit = c(30, 25), miss = c(10, 15), fa = c(8, 12), cr = c(32, 28)
)
long <- mpt_data_reshape(wide, model, id = "subject")
long
bmm_data_check(bmf(Do ~ 1, g ~ 1), long, model)

# MPTinR data: unnamed columns, sorted by tree, then by label (cr, fa for
# tree lure, then hit, miss for tree target)
wide_mptinr <- rbind(c(32, 8, 30, 10), c(28, 12, 25, 15))
identical(
  mpt_data_reshape(wide_mptinr, model, columns = "mptinr"),
  mpt_data_reshape(wide[c("cr", "fa", "hit", "miss")], model)
)

# trees with response categories that only one tree can produce
model_pc <- mpt(
  list(
    mpt_tree("pairs", list(E1 = "c * r", E2 = "1 - c * r"),
             impossible = "F1"),
    mpt_tree("singles", list(F1 = "u", E2 = "1 - u"),
             impossible = "E1")
  ),
  tree_id = "item_type"
)
wide_pc <- data.frame(
  pairs.E1 = c(20, 22), pairs.E2 = c(20, 18),
  singles.F1 = c(25, 24), singles.E2 = c(15, 16)
)
mpt_data_reshape(wide_pc, model_pc)

Create an MPT model from an EQN model file

Description

Reads an EQN model file, the format used by multiTree, MPTinR, and TreeBUGS, and returns the corresponding mpt() model object; model files in MPTinR's own "easy" format are read by mpt_from_string(). Each line of an EQN file contains three whitespace-separated fields – tree, response category, branch probability expression – and branch lines that terminate in the same response category of a tree are summed.

Usage

mpt_from_eqn(
  file,
  restrictions = NULL,
  categories = NULL,
  tree_id = NULL,
  covariates = NULL,
  simplex = NULL,
  links = "logit",
  impossible = NULL,
  tree_names = NULL
)

Arguments

Details

File format. Blank lines and lines that start with ⁠#⁠ are skipped. If the first remaining line holds a single integer and nothing else, it is read as the equation count of classical EQN files and skipped, with a warning when it differs from the number of equation lines that follow. Any other first line is read as an equation. Note that TreeBUGS always ignores the first line of an EQN file, whatever it holds, and so does MPTinR with model.type = "eqn" (its default for .eqn files; with "eqn2" it reads the first line as the equation count). bmm therefore reads a file without a count line in full, while those programs drop its first equation; a title line produces an error that says to delete it or start it with ⁠#⁠. A file that starts with the equation count, or with ⁠#⁠ lines followed by equations, is read the same way by all three. A line with fewer than three fields, or an expression that cannot be parsed, is an error that names the line; a tree whose branches do not sum to 1 at numeric test values is an error that names the tree.

Response categories. In bmm, every observation (row) belongs to one tree, and all trees share the same response categories. Classical EQN files often label the response categories per tree instead (e.g., hit/miss in the old-item tree and fa/cr in the new-item tree of a recognition model). Use the categories argument to map such tree-specific labels onto the shared response categories (e.g., the response options "yes" and "no"). If a category cannot occur in some trees (e.g., the pair and singleton trees of the pair-clustering model), declare it with the impossible argument instead.

Names. Because brms does not allow underscores or dots in non-linear parameter names, parameter and response category names are sanitized: a leading ⁠<tree>_⁠ prefix is removed from category names that are not mapped with categories, and all remaining underscores and dots are stripped. Two labels of the file, or two different values within one tree, that become the same name are an error; mapped values that become the same name in different trees are that one shared category. Names are case-sensitive, so P_G1 and P_g1 stay distinct (PG1, Pg1). Tree labels are kept, unless tree_names maps them. The formulas passed to bmm() use the bmm names (bmf(Do ~ 1) for a file parameter D_o; bmf() is short for bmmformula()); the parameter rows of the renaming map translate them, and print(model) lists them. Of the arguments that name parameters, only restrictions and covariates use the file's names.

Renaming map. The mpt_renaming attribute of the returned model is a data frame with one row per parameter, per tree, and per category of each tree, renamed or not, and the columns kind ("parameter", "tree", or "category"), tree (the bmm tree name for categories, NA otherwise), file_name, bmm_name, and renamed (TRUE where the two names differ). The names that changed are also reported in a message.

Data for an imported model. bmm() expects one row per participant and tree: a column named by tree_id holding the bmm tree names, and one count column per shared response category, named by bmm_name in the renaming map (see mpt()). A category that is impossible in a tree gets 0 (or NA) in the rows of that tree. TreeBUGS data are wide: one row per participant, one column per category label of the EQN file. MPTinR data for an EQN model are unnamed, and their columns do not follow the file: the trees and, within each tree, the category labels are sorted (alphabetically, or numerically when all labels in the file are numbers); MPTinR::check.mpt("model.eqn") returns the order in its elements eqn.order.trees and eqn.order.categories. (For MPTinR's own model format the columns follow the lines instead; see mpt_from_string().) mpt_data_reshape() reshapes both layouts with the renaming map, sums the labels that categories merged onto one response category, and gives 0 to the categories a tree cannot produce (see Examples).

bmm cannot tell a reshape that swaps tree labels or category columns from a correct one: such data fit without an error. The row sums of the count columns are the trials per row; they cannot show swapped columns, and they show swapped tree labels only when the trees differ in their number of trials. Comparing the response proportions per tree with the design (targets must get more "yes" responses than lures) catches swapped tree labels and response columns swapped in every tree, but not always a column slipped or a map entry wrong in one tree, nor a swap between two trees of the same kind; so also compare one participant's long rows with the wide row. bmm_data_check() runs the remaining checks (columns, tree labels, counts in impossible cells) without a fit and prints ⁠Hard checks: passed⁠ when they hold.

Provenance. The returned model carries an mpt_source attribute, list(file = , md5 = , lines = ): the full path of the EQN file, its MD5 checksum, and the lines as read, so that a saved fit records which file it was built from.

Value

An object of class bmmodel (see mpt()), with the attributes mpt_renaming and mpt_source (see Details).

See Also

mpt_from_string() for model text in MPTinR's own format; mpt() and mpt_tree() for building the trees in R, including the data layout bmm() expects.

Examples

eqn_file <- tempfile(fileext = ".eqn")
writeLines(c(
  "6",
  "target  hit   D_o",
  "target  hit   (1-D_o)*g",
  "target  miss  (1-D_o)*(1-g)",
  "lure    fa    (1-D_n)*g",
  "lure    cr    D_n",
  "lure    cr    (1-D_n)*(1-g)"
), eqn_file)
# restrictions use the file's names; the formulas use the bmm names
model <- mpt_from_eqn(
  eqn_file,
  restrictions = "D_n = D_o",
  categories = c(hit = "yes", fa = "yes", miss = "no", cr = "no"),
  tree_id = "item_type"
)
attr(model, "mpt_renaming")

# TreeBUGS data: one row per participant, columns named by the category
# labels of the EQN file, no participant column. mpt_data_reshape() reshapes
# them to bmm's one row per participant and tree, through the map; the
# categories a tree cannot produce get 0
wide <- data.frame(
  hit = c(30, 25), miss = c(10, 15), fa = c(8, 12), cr = c(32, 28)
)
long <- mpt_data_reshape(wide, model)
long
# the row sums are the trials per row (40 here). The proportions catch
# swapped tree labels and swapped response columns (targets must get
# more "yes" than lures), but not always a column slipped in one
# tree: also compare one participant's long rows with the wide row
rowSums(long[c("yes", "no")])
tapply(long$yes / rowSums(long[c("yes", "no")]), long$item_type, mean)
long[long$id == 1, ]
wide[1, ]
bmm_data_check(bmf(Do ~ 1, g ~ 1), long, model)

# MPTinR data for an EQN model are unnamed and sorted by tree, then by
# label: MPTinR::check.mpt(eqn_file)$eqn.order.categories lists the
# labels (here cr, fa for tree new, then hit, miss for tree old)
wide_mptinr <- rbind(c(32, 8, 30, 10), c(28, 12, 25, 15))
identical(mpt_data_reshape(wide_mptinr, model, columns = "mptinr"), long)

Create an MPT model from an MPTinR-style model string

Description

Reads a model definition in the "easy" format of MPTinR model files and returns the corresponding mpt() model object; EQN files (multiTree, MPTinR, TreeBUGS) are read by mpt_from_eqn(). Each line is one branch probability expression; blank lines separate trees. The response category of each line is given either by an inline comment (expression # category) or by the categories argument. Lines that terminate in the same response category are summed.

Usage

mpt_from_string(
  text,
  restrictions = NULL,
  tree_names,
  categories = NULL,
  tree_id = NULL,
  covariates = NULL,
  simplex = NULL,
  links = "logit",
  impossible = NULL
)

Arguments

Details

Comments. A line whose first non-blank character is ⁠#⁠ is a comment and is skipped. It does not separate trees; only blank lines do. On a line with an expression, the text after ⁠#⁠ is the response category of that line. This differs from MPTinR, which ignores everything after a ⁠#⁠ and treats a line that starts with ⁠#⁠ as a blank line, so that there comment lines can separate trees. An MPTinR model file whose trees are separated only by comment lines, or whose inline comments are not category labels, has to be edited before it is imported.

Trees with different response categories. In bmm all trees share the same response categories. If a category cannot occur in some trees (e.g., the pair and singleton trees of the pair-clustering model), declare it with the impossible argument. If the trees only label the same responses differently, give those lines the same category label.

Names. Parameter names are used as written in the model text and are case-sensitive, so the formulas passed to bmm() name them the same way (bmf(D ~ 1); bmf() is short for bmmformula()). They must be valid bmm names (letters and digits, see mpt()); unlike mpt_from_eqn(), this function does not rename them. The import fails when the branches of a tree do not sum to 1 at numeric test values.

Data for an imported model. The Details of mpt_from_eqn() describe the layout bmm() expects (one row per participant and tree, a tree column, one count column per response category, 0 for impossible categories) and the checks to run before fitting; they apply to models from this function as well, with the tree names of tree_names and the category labels of the lines. The data MPTinR fits to a model file in this format are wide and positional: one row per participant (or one vector for aggregated counts) and one unnamed column per line of the model file, in file order, because each line is one response category. (For EQN files MPTinR sorts the columns instead; see mpt_from_eqn().) Name the columns ⁠<tree>.<category>⁠; mpt_data_reshape() then stacks the trees and gives 0 to the categories a tree cannot produce (see Examples).

Provenance. The returned model carries an mpt_source attribute, list(text = ), with the model definition as one string.

Value

An object of class bmmodel (see mpt()), with an mpt_source attribute holding the imported text.

See Also

mpt_from_eqn() for EQN files (multiTree, MPTinR, TreeBUGS); mpt() and mpt_tree() for building the trees in R, including the data layout bmm() expects.

Examples

model_2htm <- "
# two-high-threshold model
D + (1 - D) * g        # old
(1 - D) * (1 - g)      # new

(1 - D) * g            # old
D + (1 - D) * (1 - g)  # new
"
model <- mpt_from_string(
  model_2htm,
  tree_names = c("target", "lure"),
  tree_id = "item_type"
)
model

# fix the guessing rate as in MPTinR / TreeBUGS
mpt_from_string(
  model_2htm,
  tree_names = c("target", "lure"),
  tree_id = "item_type",
  restrictions = "g = 0.5"
)

# pair clustering, written as in MPTinR without category labels: pairs and
# singletons have different response categories
model_pc <- "
c * r
(1 - c) * u * u
2 * (1 - c) * u * (1 - u)
c * (1 - r) + (1 - c) * (1 - u) * (1 - u)

u
1 - u
"
mpt_from_string(
  model_pc,
  tree_names = c("pairs", "singles"),
  categories = list(
    pairs = c("E1", "E2", "E3", "E4"),
    singles = c("F1", "F2")
  ),
  tree_id = "item_type",
  impossible = list(pairs = c("F1", "F2"), singles = c("E1", "E2", "E3", "E4"))
)

# an MPTinR model file read with readLines() keeps the blank lines that
# separate the trees. A restrictions file is not read by bmm: drop its
# comments and blank lines and pass the remaining lines.
model_file <- tempfile(fileext = ".model")
writeLines(model_2htm, model_file)
restr_file <- tempfile(fileext = ".restr")
writeLines(c("# fix the guessing rate", "g = 0.5"), restr_file)
restr <- trimws(sub("#.*", "", readLines(restr_file)))
model <- mpt_from_string(
  readLines(model_file),
  tree_names = c("target", "lure"),
  tree_id = "item_type",
  restrictions = restr[nzchar(restr)]
)

# MPTinR data: one unnamed column per line of the model file, in file
# order. Name the columns by tree and category, then let mpt_data_reshape()
# stack the trees to get bmm's one row per participant and tree.
wide <- rbind(c(30, 10, 8, 32), c(25, 15, 12, 28))
colnames(wide) <- c("target.old", "target.new", "lure.old", "lure.new")
long <- mpt_data_reshape(wide, model)
long
# the row sums are the trials per row (40 here), but they cannot show
# swapped columns or tree labels: compare one participant with the wide
# row, then run the data checks (bmm_data_check() prints
# "Hard checks: passed" when the columns, tree labels and counts fit)
rowSums(long[c("old", "new")])
long[long$id == 1, ]
wide[1, ]
bmm_data_check(bmf(D ~ 1), long, model)

Create a tree for a Multinomial Processing Tree (MPT) model

Description

Specifies a single tree of a multinomial processing tree model by naming the tree and providing one branch-probability expression per response category. Trees created with mpt_tree() are combined into a full model specification with mpt().

Usage

mpt_tree(name, branches, ..., impossible = NULL)

Arguments

Details

Nothing in the formulas is evaluated: the left-hand side must be a single response category name, the right-hand side is kept as an expression (with constant arithmetic folded, see below), and the formula environment is dropped. A one-sided formula, a left-hand side that is not a name, and a mix of formulas and strings are errors.

Constant arithmetic such as 1/4 or 1/(2*2) is folded into a decimal literal (0.25) when the tree is created. Stan compiles a bare integer fraction as integer division (1/4 == 0), which would silently corrupt the likelihood.

Categories listed in impossible keep their column in the data and their place in the response matrix, but their probability is set to zero for the observations belonging to this tree. Because the model is estimated on the log scale, "zero" is implemented as a large negative linear predictor rather than a literal zero; no probability has to be taken away from the remaining branches to compensate.

Value

An object of class mpt_tree

Examples

# the target tree of the two-high-threshold model: targets (studied items)
# end in the responses "old" and "new"
tree_target <- mpt_tree(
  "target",
  old ~ D + (1 - D) * g,
  new ~ (1 - D) * (1 - g)
)
tree_target

# the same tree with the branches as a named list of strings
tree_target_strings <- mpt_tree(
  name = "target",
  branches = list(
    old = "D + (1 - D) * g",
    new = "(1 - D) * (1 - g)"
  )
)
identical(tree_target, tree_target_strings)

# a three-alternative recognition task in which lures are shown on some
# trials only: `lure` responses cannot occur on the others
tree_lures <- mpt_tree(
  "lures",
  target ~ D * b + (1 - D) * (1/3),
  lure ~ D * (1 - b) + (1 - D) * (1/3),
  new ~ (1 - D) * (1/3)
)

# the branches of the lure-free tree sum to 1 without the `lure` category;
# no probability is reserved for it
tree_nolures <- mpt_tree(
  "nolures",
  target ~ D + (1 - D) * (1/2),
  new ~ (1 - D) * (1/2),
  impossible = "lure"
)
tree_nolures

# a `tree` column in the data holds "lures" or "nolures" per observation;
# effects of experimental conditions go in the parameter formulas instead,
# e.g. bmf(D ~ 0 + cond)
model <- mpt(list(tree_lures, tree_nolures), tree_id = "tree")

Distribution functions for Multinomial Processing Tree (MPT) models

Description

Density and random generation functions for multinomial processing tree models specified with mpt(). Please note that these functions are currently not vectorized.

Usage

dmpt(x, pars, mpt_model, tree = NULL, log = TRUE, covariates = NULL, ...)

rmpt(
  n,
  size,
  pars,
  mpt_model,
  tree = NULL,
  unpack = FALSE,
  covariates = NULL,
  ...
)

Arguments

Value

dmpt gives the multinomial density of a vector of response counts, and rmpt generates random response counts for the response categories.

Note

Unlike the densities of the circular models in this package (dsdm(), dmixture2p(), dmixture3p(), dimm()) and unlike stats::dmultinom(), dmpt() returns the log density by default, matching dm3().

References

Batchelder, W. H., & Riefer, D. M. (1999). Theoretical and empirical review of multinomial process tree modeling. Psychonomic Bulletin & Review, 6(1), 57-86. https://doi.org/10.3758/BF03210812

Examples

tree_target <- mpt_tree("target", list(
  old = "D + (1 - D) * g",
  new = "(1 - D) * (1 - g)"
))
tree_lure <- mpt_tree("lure", list(
  old = "(1 - D) * g",
  new = "D + (1 - D) * (1 - g)"
))
model <- mpt(list(tree_target, tree_lure), tree_id = "item_type")

dmpt(
  x = c(35, 15), pars = c(D = 0.7, g = 0.5),
  mpt_model = model, tree = "target"
)
rmpt(
  n = 10, size = 50, pars = c(D = 0.7, g = 0.5),
  mpt_model = model, tree = "lure"
)

Metacognitive efficiency (M-ratio) from a fitted meta-d' SDT model

Description

Extracts the posterior of the M-ratio (\mathrm{meta\text{-}d'}/d') and of meta-d' from a sdt_rating() model fit with version = "metad". That model estimates the log M-ratio (logmratio) directly, so the M-ratio is exp(logmratio) and meta-d' is exp(logmratio) * d. The M-ratio is the standard measure of metacognitive efficiency (Maniscalco & Lau, 2012; Fleming, 2017): 1 is ideal metacognition, below 1 is inefficiency, and above 1 is hyper-efficiency. Type-1 and type-2 sensitivity are reported on one scale, d' or, once sdratio is estimated, d_a (see sdt_sensitivity()), so the ratio is unaffected by sdratio. The type-1 criterion keeps its location in the meta-d' space, as in the HMeta-d model equations (Fleming, 2017, Appendix); maximum-likelihood meta-d' constrains meta-c' = c' (Maniscalco & Lau, 2014), so the two differ when the criterion is far from the midpoint.

Usage

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

Arguments

Value

A data frame of class "bmm_sdt_mratio" summarising the posterior: columns parameter ("mratio" or "metad"), mean, median, lower, upper (the credible-interval bounds at probs), and any condition columns. The full per-draw posteriors are kept in the draws attribute (long format: parameter, value, .draw, plus condition columns) for downstream computation such as plotting or condition contrasts.

References

Maniscalco, B., & Lau, H. (2012). A signal detection theoretic approach for estimating metacognitive sensitivity from confidence ratings. Consciousness and Cognition, 21(1), 422–430. doi:10.1016/j.concog.2011.09.021

Maniscalco, B., & Lau, H. (2014). Signal detection theory analysis of type 1 and type 2 data: meta-d', response-specific meta-d', and the unequal variance SDT model. In S. M. Fleming & C. D. Frith (Eds.), The cognitive neuroscience of metacognition (pp. 25–66). Springer.

Fleming, S. M. (2017). HMeta-d: hierarchical Bayesian estimation of metacognitive efficiency from confidence ratings. Neuroscience of Consciousness, 2017(1), nix007. doi:10.1093/nc/nix007

See Also

sdt_rating(), latent_sdt()

Posterior draws of model parameters on the native scale

Description

bmm samples all parameters on their link scale (log for precision and boundary parameters, logit or softmax for mixture weights, tan_half for circular locations). native_parameters() returns posterior draws of those parameters on the scale they are interpreted and reported on, evaluated over a grid of predictor values.

Because draws are returned rather than summaries, any contrast is ordinary arithmetic on the draws, with the baseline made explicit by the grid.

Usage

native_parameters(
  x,
  newdata = NULL,
  pars = NULL,
  re_formula = NULL,
  scale = c("native", "sampling"),
  ndraws = NULL,
  draw_ids = NULL,
  summary = FALSE,
  prob = 0.95,
  robust = FALSE,
  ...
)

Arguments

Value

If summary = FALSE, a data.frame with one row per draw, grid cell and parameter, with columns .chain, .iteration, .draw, the grid variables, parameter and value. The draw indices refer to the draws of x, so they can be joined with the output of posterior::as_draws_df(x). If summary = TRUE, a data.frame with one row per grid cell and parameter, with columns Estimate, Est.Error and the interval bounds.

Transform first, then summarise

Quantiles are preserved by monotone inverse links but means are not, so the median and the credible interval of a native-scale parameter are exact while its mean is not the inverse link of the linear predictor's mean. This function transforms the draws and only then summarises, which is why summary = TRUE is not the same as transforming the output of summary.bmmfit().

The same order is what makes summary = TRUE correct for the mixture weights, whose softmax is not an elementwise map at all: summarising the transformed draws reports quantiles of the weight's own marginal posterior, whereas transforming the summaries would report the softmax of three separate quantiles, which is not a quantile of anything.

Transforming coefficients is not the same thing

The inverse link applies to the linear predictor of a grid cell, not to an individual regression coefficient. For a model with kappa ~ condition, exp(b_kappa_conditionB) is a multiplicative factor, not kappa in condition B; kappa in condition B is exp(b_kappa_Intercept + b_kappa_conditionB). Building the grid is exactly what removes this step, which is why contrasts are taken between rows of the output rather than read off the coefficients.

Group-level effects

re_formula = NULL (the default) returns subject-specific parameters and adds the grouping variables to the grid. re_formula = NA sets all group-level effects to zero, which under a non-identity link gives the median subject, not the population mean. The population mean requires marginalising over the distribution of group-level effects, which is done by predicting for a new level:

native_parameters(fit, newdata = transform(nd, id = "new"),
                  allow_new_levels = TRUE, sample_new_levels = "gaussian")

These three quantities differ substantially at between-subject standard deviations typical for working memory data.

What is returned

Parameters fixed to a constant are returned at that constant, transformed to the native scale. Because bmm fixes parameters on the link scale, the native value can differ from the value shown by parameter_info(): the ddm relative starting point zr is fixed at 0 under a logit link and is therefore reported as 0.5, and the ezdm and cswald diffusion constant s is fixed at 0 under a log link and is reported as 1.

Trial-level derived quantities, such as the response probabilities of the mixture models or the category probabilities of m3, are not returned; use brms::posterior_epred() for those.

The default grid contains the observed combinations of the predictors, not their full crossing, so cells that were never presented do not appear. A continuous predictor therefore produces one grid cell per observed value; supply newdata for such models.

Circular location parameters

The circular models sample their location parameter (mu for sdm, mu1 for mixture2p, mixture3p and imm) through a tan_half link, and it is returned in radians in ⁠(-pi, pi)⁠ — a response bias relative to the target, since the response variable is the angular error. It is fixed to 0 unless the bmmformula predicts it explicitly, so an all-zero mu1 means the model never estimated one.

The inverse link is 2 * atan(), and the caveat above applies to it with particular force: 2 * atan(b_mu1_conditionB) is neither the bias in condition B nor the difference between conditions. The bias in condition B is 2 * atan(b_mu1_Intercept + b_mu1_conditionB), and the difference between the conditions is that value minus 2 * atan(b_mu1_Intercept), taken draw by draw.

Mixture weights

For mixture3p the returned thetat and thetant are probabilities but they do not sum to 1: brms holds the linear predictor of one mixture component at zero, and that component is the guessing distribution, so the remaining mass 1 - thetat - thetant is the probability of a guess. thetant is the total probability of a non-target response, summed over the lures, not the probability per lure.

At set sizes where no non-target was presented the model switches its non-target components off, so thetant is reported as exactly 0 and thetat as plogis() of its linear predictor, rather than as the value the set-size regression extrapolates to. This correction reads the lure indicators out of the prediction grid; if you supply newdata, columns you do not supply are filled from the first row of the model data rather than recomputed, and native_parameters() warns when the two disagree.

See Also

parameter_info(), native_transform(), conditional_effects.bmmfit()

Examples

fit <- bmm(
  bmf(c ~ 0 + set_size, kappa ~ 1),
  data = oberauer_lin_2017,
  model = sdm(resp_error = "dev_rad")
)

# draws of c and kappa for every set size, on the native scale
np <- native_parameters(fit, re_formula = NA)
head(np)

# a contrast is plain arithmetic on the draws
c_draws <- subset(np, parameter == "c")
quantile(
  c_draws$value[c_draws$set_size == 1] - c_draws$value[c_draws$set_size == 4],
  probs = c(0.025, 0.5, 0.975)
)

# posterior summaries of the transformed draws
native_parameters(fit, re_formula = NA, summary = TRUE)

Transform parameter draws from the sampling scale to the native scale

Description

Applies a bmmodel's inverse link functions to posterior draws of its parameters. This is the extension point used by native_parameters(): the transformation is a property of the model, not of an individual parameter, because some models map several parameters jointly (e.g. mixture weights through a softmax).

A model whose parameters transform elementwise, or through one softmax group that is active on every row of the data, needs no method — declaring links in its ⁠.model_*()⁠ constructor is sufficient, and that covers every model bmm currently ships. Write a method when either condition fails. Two cases that are known to fail it: a model with several independent multinomial branches, because a "softmax" link carries no group identity and all such parameters are softmaxed together as one group; and a transformation that depends on the design rather than only on the parameter's own value.

Usage

native_transform(model, linpred, data, ...)

## Default S3 method:
native_transform(model, linpred, data, ...)

## S3 method for class 'non_targets'
native_transform(model, linpred, data, ...)

Arguments

Details

Methods for specific models follow the bmm S3 chain, which runs from the general classes to the specific ones. A method should therefore transform the parameters it owns, remove them from linpred, and call NextMethod() so the remaining parameters reach the default method:

native_transform.mymodel <- function(model, linpred, data, ...) {
  own <- .my_joint_transform(linpred[c("p1", "p2")])
  linpred <- linpred[not_in(names(linpred), c("p1", "p2"))]
  c(NextMethod(), own)
}

Removing them before NextMethod() is what keeps a method idempotent with respect to the parameters it does not own: a parameter left in linpred is transformed a second time by the next method in the chain, which does not error and does not warn.

Value

A named list of matrices on the native scale. A method must return one element for every name it was given and preserve each matrix's dimensions; native_parameters() checks both and errors, naming the model, if either is broken. The order of the list is free.

See Also

native_parameters()

Data from Experiment 1 reported by Oberauer & Lewandowsky (2019)

Description

Raw data of 40 subjects that completed a verbal memory recall task in three different conditions using different types of distractor words.

Usage

oberauer_lewandowsky_2019_e1

Format

oberauer_lewandowsky_2019_e1

A data frame with 120 rows and 10 columns:

ID

Integer uniquely identifying each subject

cond

Factor sperating the three experimental conditions: ⁠new distractors⁠ refers to new words being used as distractors, ⁠old reordered⁠ refers to the to be remembered words being the distractors, but reordered relative to the serial position, ⁠old same⁠ refers to the to be remebered words being the distractors, and appearing in the same order as the to be remembered words.

corr

The frequency a subject recalled the correct item

other

The frequency a subject recalled one of the other to be remebered words

dist

The frequency a subject recalled one of the distractors

npl

The frequency a subject recalled a not-presented lure (NPL), that is a word that was not presented during a trial

n_corr, n_other, n_dist, n_npl

The number of candidataes in each of the response categories

Data from Experiment 1 reported by Oberauer & Lin (2017)

Description

Raw data of 19 subjects that completed a continuous reproduction task with set size 1 to 8 reported by Oberauer & Lin (2017).

Usage

oberauer_lin_2017

Format

oberauer_lin_2017

A data frame with 15,200 rows and 19 columns:

ID

Integer uniquely identifying different subjects

session

Session number

trial

Trial number within each session

set_size

The set_size of the data in this row

dev_rad

The response error, that is the difference between the response given and the target color in radians.

col_nt1, col_nt2, col_nt3, col_nt4, col_nt5, col_nt6, col_nt7

The non-target items' color value relative to the target.

dist_nt1, dist_nt2, dist_nt3, dist_nt4, dist_nt5, dist_nt6, dist_nt7, dist_nt8

The spatial distance between all non-target items and the target item in radians.

Source

https://osf.io/m4shu

Get parameter information for a bmm model

Description

Returns a data frame with information about the model parameters: their descriptions, whether they are fixed (and to which value), and their link functions. Use prior_info() for the priors of a fitted model.

Usage

parameter_info(x, ...)

## S3 method for class 'bmmodel'
parameter_info(x, formula = NULL, ...)

## S3 method for class 'bmmfit'
parameter_info(x, ...)

Arguments

Value

A data frame of class bmm_parameters with one row per parameter and columns: parameter, description, fixed, value, and link.

Examples

# For an unfitted model
parameter_info(sdm(resp_error = "y"))

# For an M3 model
parameter_info(m3(
  resp_cats = c("corr", "other", "npl"),
  num_options = c(1, 4, 5),
  version = "ss"
))

Deprecated: use parameter_info()

Description

parameters() is deprecated as of bmm 1.4.0 and will be removed in bmm 1.6.0. It shares its name with parameters::parameters(), so whichever of the two packages is attached last decides what parameters() returns. Use parameter_info() instead; it takes the same arguments and returns the same table.

Usage

parameters(x, ...)

Arguments

Value

The output of parameter_info().

Plot the posterior AUC distribution from a fitted SDT model

Description

Plot the posterior AUC distribution from a fitted SDT model

Usage

## S3 method for class 'bmm_sdt_auc'
plot(x, condition_col = NULL, ...)

Arguments

Value

A ggplot2 object.

See Also

auc_sdt(), roc_sdt()

Plot the latent decision-variable distributions of an SDT model

Description

Draws the model-implied noise and signal evidence densities on the latent decision axis. For sdt_yn()/sdt_rating() the response criterion or the K-1 confidence thresholds are added as dashed vertical lines with a shaded credible band; when several boundaries are shown together (e.g. base-rate criteria or the rating thresholds) they are colour-coded. sdt_mafc() and sdt_ranking() have no boundary; with show_competitors = TRUE in latent_sdt() the max-of-distractors densities are overlaid as dashed lines, one per set size. Distinct d/sdratio conditions are faceted.

Usage

## S3 method for class 'bmm_sdt_latent'
plot(x, condition_col = NULL, line_alpha = 0.9, ...)

Arguments

Value

A ggplot2 object.

See Also

latent_sdt(), roc_sdt()

Plot a model-implied SDT ROC curve

Description

Shows the smooth model-implied ROC as a credible-band ribbon, with the model-implied operating points overlaid as colour-coded markers with crosshair error bars (uncertainty in both directions). For binary multi-criteria fits the points are the several criterion levels; for rating models they are the K-1 confidence thresholds (labelled t1..t(K-1)). Pass observed to overlay empirical points from roc_observed().

Usage

## S3 method for class 'bmm_sdt_roc'
plot(
  x,
  observed = NULL,
  condition_col = NULL,
  add_diagonal = TRUE,
  scale = c("probability", "quantile", "z"),
  ribbon_alpha = 0.25,
  point_size = 2.5,
  ...
)

Arguments

Details

With scale = "quantile" (or its alias scale = "z") the rates are read on the distribution's quantile axis via the transform -qf(1 - rate), where qf is the inverse CDF of the fitted noise distribution, and the axis label names that transform (z for normal, logit for logistic, loglog/cloglog for the Gumbel distributions). On this axis the binary model ROC is a straight line with slope 1 / exp(sdratio) and intercept d * sqrt((1 + exp(sdratio)^2) / 2) / exp(sdratio) – the separation in noise-SD units over the signal SD, since d is d_a – so a slope below 1 is the unequal-variance signature (signal SD > noise SD), and departures of the observed points from a straight line diagnose misfit. This linearity holds for all four distributions. The (0,0) and (1,1) endpoints map to infinity and are dropped on the transformed scale.

Value

A ggplot2 object. The credible band (ribbon) is the pointwise posterior interval of the hit rate at each criterion value, plotted at that criterion's posterior-mean false-alarm rate; it does not include uncertainty in the false-alarm rate.

See Also

roc_sdt(), roc_observed(), auc_sdt()

Plot the structure of MPT trees

Description

Draws the processing-tree diagram implied by the branch probability expressions: each branch expression is expanded into its root-to-leaf paths (a sum of products, with the factors in writing order), paths with a common prefix share the corresponding tree nodes, and the leaves show the response categories. Calling plot() on a model created with mpt() draws all trees of the model.

Usage

## S3 method for class 'mpt_tree'
plot(x, cex = 0.9, ...)

## S3 method for class 'mpt'
plot(x, cex = 0.9, ...)

Arguments

Details

Edges are drawn in the order the branch expressions first mention them, bare constants included. A complement edge 1 - X is placed directly below its sibling X; the complement of a product, 1 - PG1 * PG2, is one atomic edge. A numeric constant factor above 1 is not a separate edge: the constant and all factors after it form one edge, as in Batchelder & Riefer (1999, Fig. 1). For example, B * 2 * u * (1 - u) draws one edge B, then one edge 2 * u * (1 - u). If the constant is at the start (2 * (1 - c) * u * (1 - u)), the entire expression is one root edge beside its siblings, drawn without a warning, so write the constant after the process steps it follows, as in (1 - c) * 2 * u * (1 - u). If it is in the middle (u * 2 * (1 - u)), the constant and the factors after it join into one edge, and a trailing constant joins the factor before it. The plot warns when this merging hides a process step: a factor after the constant that has no complement among the merged factors (D * 2 * g * (1 - g) * r), a factor separated from its complement by the constant (u * 2 * (1 - u)), or a quotient whose denominator does not contain the complemented parameter ((1 - Pi) * 1 / 15 is one edge; write (1 - Pi) * (1/15)). A quotient normalised by a sum that contains the parameter, as in x * (1 - d) / (x * (1 - d) + y), does not warn. Expressions are expanded distributively, so Pm * (Pb + (1 - Pb) * 0.25) displays two paths sharing the Pm edge.

The plot warns when edges leaving an internal node do not sum to 1, but only for nodes whose edge labels are all constants, plain parameters, or 1 - parameter; product edges, quotient edges, function calls, and covariate fans that lack complements are not checked, nor are nodes with a single edge or, when plotting a bare tree, fans of plain parameters (which a model may declare as a simplex). Write factors in the same order across all categories to avoid mispairing due to different factor orders. A bare mpt_tree() whose constants were already substituted (such as the stored restricted trees of a model) does not recognise a complement written as a number (0.8 beside 0.2), so those edges keep the order of first mention; plot() of a model draws the unrestricted trees and is not affected.

For a model created with mpt(), the plot draws the trees as written and marks restricted parameters: the edge of a fixed parameter reads g = 0.5, the edge of an equated parameter ⁠Dn (= Do)⁠, and an edge that contains a restricted parameter shows the restriction in parentheses, as in 1 - g (g = 0.5). Impossible categories are omitted. A literal 0 in an expression (as in D + 0 * x) is drawn as an edge labelled 0. The ... argument is ignored; the title cannot be changed. Every product of sums is expanded into all its paths, so each multiplied sum such as (a + b) doubles the number of paths: trees with more than about ten multiplied sums draw slowly. Leaf labels wider than about 90% of a panel are clipped, and long edge labels near the root can be clipped earlier, when the leaf labels leave little room.

Value

The object x, invisibly.

References

Batchelder, W. H., & Riefer, D. M. (1999). Theoretical and empirical review of multinomial process tree modeling. Psychonomic Bulletin & Review, 6(1), 57-86. https://doi.org/10.3758/BF03210812

Examples

tree_target <- mpt_tree("target", list(
  old = "D + (1 - D) * g",
  new = "(1 - D) * (1 - g)"
))
plot(tree_target)

model <- mpt(
  list(
    tree_target,
    mpt_tree("lure", list(
      old = "(1 - D) * g",
      new = "D + (1 - D) * (1 - g)"
    ))
  ),
  tree_id = "item_type"
)
plot(model)

Expected response of a fitted bmm model

Description

posterior_epred() returns draws of the expected value of the response for each observation: the mean of what brms::posterior_predict() simulates for it. What that is depends on the model:

With dpar or nlpar, posterior_epred() returns draws of that model parameter for every model, as in brms: a dpar on its native scale, an nlpar on the scale of its link. native_parameters() returns the model parameters on their native scale over a grid of predictor values.

Usage

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

## S3 method for class 'bmmfit'
fitted(object, ..., scale = c("response", "linear"), dpar = NULL, nlpar = NULL)

Arguments

Value

A draws by observations matrix (an array with a third dimension for the response categories of the multinomial models).

See Also

brms::posterior_epred(), native_parameters()

Examples

## Not run: 
fit <- bmm(
  bmf(drift ~ 1, bound ~ 1, ndt ~ 1),
  data = rddm(200, drift = 1.5, bound = 1.2, ndt = 0.3),
  model = ddm(rt = "rt", response = "response"),
  backend = "cmdstanr"
)
# expected response time of each observation, one row per draw
epred <- posterior_epred(fit)

## End(Not run)

Generic S3 method for postprocessing the fitted brm model

Description

Called by bmm() to automatically perform some type of postprocessing depending on the model type. It will call the appropriate postprocess_brm.* methods based on the list of classes defined in the .model_* functions. For models with several classes listed, it will call the functions in the order they are listed. Thus, any operations that are common to a group of models should be defined in the appropriate postprocess_brm.* function, where \* corresponds to the shared class. For example, for the sdm model, the postprocessing involves setting the link function for the c parameter to "log", because it was coded manually in the stan code, but it was specified as "identity" in the brms custom family. If your model requires no postprocessing, you can skip this method, and the default method will be used (which returns the same brmsfit object that was passed to it).

Usage

postprocess_brm(model, fit, ...)

Arguments

Value

An object of class brmsfit, with any necessary postprocessing applied

Examples

fit <- readRDS("my_saved_fit.rds")
postprocessed_fit <- prostprocess_brm(fit)

Posterior predictive check for bmmfit objects

Description

For models where brms provides pp_check support, this method delegates to brms::pp_check(). For models with multinomial families (e.g., the m3 model), brms's pp_check is unavailable; this method dispatches to a model-specific visualisation instead.

Usage

## S3 method for class 'bmmfit'
pp_check(
  object,
  type = NULL,
  ndraws = NULL,
  group = NULL,
  resp_var = NULL,
  ...
)

Arguments

Details

For multinomial models, the plot mirrors the bayesplot ppc_bars style: observed proportions are shown as bars and posterior predictive medians with credible intervals are shown as point-ranges, using the bayesplot default colour scheme and theme.

Some models describe several observables jointly (e.g. ddm(): response times and responses; ezdm(): mean RT, RT variance and accuracy), but brms only ever checks the primary response. For these models the resp_var argument selects which observable to check; pp_check_vars() lists the available checks. All selected observables are drawn from one joint posterior predictive simulation, so resp_var = "all" panels are mutually consistent.

Some observables are undefined for some cells — an ezdm() boundary has no mean response time when fewer than two responses reach it. Observations whose observed value is undefined are dropped from the check; undefined values in the simulated replicates are absorbed by dropping those posterior draws instead, so the number of observations checked does not depend on ndraws. Dropped draws are not missing at random — they are draws whose parameters made a boundary sparse — so the retained predictive is mildly conditioned; both reductions are reported with a warning. With resp_var = "all" one reduction is shared by every panel, so the panels are computed on the same observations and draws.

Value

For multinomial models, for a 4-parameter ezdm() fit with placeholders in mean_rt_upper, or when resp_var is specified, a ggplot2 object (a bayesplot_grid for resp_var = "all"). For other models, the result of brms::pp_check().

See Also

brms::pp_check(), pp_check_vars()

List the posterior predictive checks available for a model fit

Description

Lists the values resp_var accepts in pp_check.bmmfit() for a fitted model whose likelihood involves several observables.

Usage

pp_check_vars(fit)

Arguments

Value

A data.frame with one row per available check (columns resp_var, label, default_type, slot, listing the brms standata slots the check reads, and default, flagging the observable that pp_check() plots when resp_var is not specified), or NULL invisibly for models without multi-observable support.

See Also

pp_check.bmmfit()

Examples

## Not run: 
fit <- bmm(bmmformula(drift ~ condition), data, ddm(rt = "rt", response = "response"))
pp_check_vars(fit)
pp_check(fit, resp_var = "response")

## End(Not run)

Declare the observables of a bmm model for posterior predictive checks

Description

pp_observables() returns the model's observable declaration used by pp_check.bmmfit() when resp_var is specified, or NULL for models that delegate fully to brms::pp_check(). pp_simulate() draws all observables jointly from the posterior predictive distribution.

Usage

pp_observables(model)

pp_simulate(model, prep)

Arguments

Details

A pp_observables() method returns list(observed, checks):

• observed: named character vector mapping observable names to slots of the brms standata ("Y", "vreal1", "vint1", "trials", "dec"). The observable mapped to "Y" is the default check.

• checks: named list of entries built by the internal .pp_observable() constructor. Each compute closure receives a named list keyed by names(observed) and must be elementwise, so the identical closure produces y from length-N vectors and yrep from ndraws x N matrices.

Two optional elements serve observed data that holds placeholders rather than observations, such as the summaries of an unused ezdm() boundary:

• defaults: named vector giving, for observables whose slot a fit saved by an older bmm version lacks, the value to use for every observation.

• y_placeholders: a function of the fit's data that returns TRUE if the "Y" slot holds placeholders. brms would plot them as data, so pp_check.bmmfit() without resp_var then checks the observable mapped to "Y" itself.

A pp_simulate() method returns a named list of ndraws x nobs matrices drawn jointly, typically through the internal .pp_simulate_joint() helper around the model's ⁠r*()⁠ function. Simulating observables independently would break their joint distribution (e.g. rt and response under the DDM). Names not in observed are ignored; declared observables that are not simulated (design quantities such as trial counts) are filled in from the data.

Register exactly one method per model at the most general class level where the declaration is identical across versions.

These two generics are exported so that model methods defined outside bmm can be registered against them, but they are an internal developer interface documented for bmm's own model authors and carry no stability guarantee across releases.

Value

pp_observables() returns NULL for a model that delegates fully to brms::pp_check(), or a list with elements observed (a named character vector mapping observable names to brms standata slots), checks (a named list of check definitions, each with a compute closure, a label and a default plot type: a bayesplot ⁠ppc_*⁠ type or bmm's "bars_binned") and, optionally, defaults and y_placeholders (see Details). pp_simulate() returns a named list of ndraws x nobs matrices, one per simulated observable.

Posterior predictive p-values for count models

Description

Computes the goodness-of-fit statistics T_1 and T_2 of Klauer (2010) and their posterior predictive p-values (PPP) for models that fit aggregated response counts: mpt(), m3(), sdt_yn(), sdt_rating(), sdt_mafc(), sdt_ranking() and sdt_cdp(). TreeBUGS' PPP() reports the same two statistics; Details says where their values differ.

Usage

ppp(fit, ...)

## S3 method for class 'bmmfit'
ppp(
  fit,
  statistic = c("T1", "T2"),
  by = NULL,
  ndraws = NULL,
  level = c("group", "participant"),
  participant = NULL,
  ...
)

Arguments

Details

Both statistics are computed once per posterior draw, for the observed counts and for counts simulated from that draw with brms::posterior_predict(), each against the expected counts of the same draw from brms::posterior_epred(). The simulated data are new data from the same participants. The PPP is the proportion of draws in which the simulated counts deviate more from the expectation than the observed counts do. A PPP near 0 flags misfit; a PPP away from 0 means that the model reproduces the aspect of the data the statistic measures.

T_1: mean structure. The rows of the data are summed into cells, and T_1 is Pearson's X^2 over all cells and response categories,

T_1 = \sum_{c,k} \frac{(n_{ck} - \hat{n}_{ck})^2}{\hat{n}_{ck}},

where n_{ck} is the observed and \hat{n}_{ck} the expected count of category k in cell c. Klauer (2010, Eq. 17) and TreeBUGS average the counts over participants instead of summing them. With the same number T of participants in every cell, the two give the same PPP and T_1 here is T times theirs; otherwise the mean gives the cells with fewer participants more weight.

T_2: covariance structure. For every participant, the counts are summed per cell, giving one vector per participant over all cells and categories. T_2 compares the covariance matrix s of these vectors across participants with the covariance matrix \sigma the model expects (Klauer, 2010, Eq. 18):

T_2 = \sum_{i,j} \frac{(s_{ij} - \sigma_{ij})^2}{\sqrt{\sigma_{ii}\sigma_{jj}}}.

\sigma is the covariance of the expected counts across participants plus the average multinomial covariance within a participant, which is the expected value of s given the draw. Both covariances divide by T - 1, where T is the number of participants. Klauer (2010) multiplies the multinomial part by (T - 1) / T, the expected value for covariances that divide by T; that version of T_2 is (T - 1) / T times this one and has the same PPP. TreeBUGS combines the factor with covariances that divide by T - 1, so its \sigma falls short of the expected s by 1 / T of the multinomial part. The sum runs over the full matrix, so the covariances between trees or conditions, which carry the correlations between the person parameters, are part of the check. When participants differ in the cells they have, as in a between-participants design, they are grouped by their set of cells, T_2 is computed per group, and the results are summed. A participant who shares the set of cells with nobody else has no covariance and is left out, with a message.

Cells. By default, a cell is a combination of the tree (mpt()), the stimulus (the SDT models) or the set size (sdt_ranking()), and the categorical population-level predictors of the model formula. Numeric predictors are not part of the default, because every distinct value would form a cell of its own. Rows of a cell are summed, so the trees or stimuli must be part of by for the statistic to compare each tree's own categories.

Impossible categories. Categories with an expected count of about 0 in every draw, such as the impossible categories of an mpt() or m3() model or the ranks beyond the set size of an sdt_ranking() model, are left out of both statistics. A cell with a single possible category, such as an m3() set size without other candidates, is left out of T_2: the number of trials fixes its count.

Value

A data frame of class bmm_ppp with one row per statistic (and per participant for level = "participant") and the columns statistic, observed and replicated (the means of the statistic over draws, for the observed and the simulated counts, as reported by Klauer, 2010, and TreeBUGS), ppp and ndraws. The draws of the statistics are stored in the attribute "draws".

References

Klauer, K. C. (2010). Hierarchical multinomial processing tree models: A latent-trait approach. Psychometrika, 75(1), 70–98. https://doi.org/10.1007/s11336-009-9141-0

See Also

pp_check.bmmfit() for the graphical check.

Examples

trees <- list(
  mpt_tree("target", list(old = "D + (1 - D) * g", new = "(1 - D) * (1 - g)")),
  mpt_tree("lure", list(old = "(1 - D) * g", new = "D + (1 - D) * (1 - g)"))
)
dat <- expand.grid(id = factor(1:20), item_type = c("target", "lure"),
                   stringsAsFactors = FALSE)
dat$old <- rbinom(nrow(dat), 50, ifelse(dat$item_type == "target", 0.8, 0.2))
dat$new <- 50 - dat$old

fit <- bmm(
  bmf(D ~ 1 + (1 | id), g ~ 1 + (1 | id)), dat, mpt(trees, "item_type"),
  cores = 4, backend = "cmdstanr"
)
ppp(fit)
ppp(fit, level = "participant")

Generic S3 method for model-specific lines in print() of a model

Description

Called by print() on a bmmodel after the lines all models share and before the closing pointer to parameter_info(). The default method prints nothing.

Usage

print_model_details(model, ...)

Arguments

Value

Called for its printed output; returns NULL invisibly

Generate a markdown list of the measurement models available in bmm

Description

Used internally to populate the README and the "Get started" article. Models are grouped as in bmm_models(), and every model links to its reference page on the website.

Usage

print_pretty_models_md(group = NULL)

Arguments

Value

Markdown code for printing the list of measurement models available in bmm

Examples

print_pretty_models_md()
print_pretty_models_md(group = "Continuous reproduction")

Priors of a fitted bmm model and where they came from

Description

Lists the priors a fitted bmm model was actually sampled with and where each one came from: a bmm default, a brms default, or a user-specified prior. Each row also gives the link function of the parameter, because the prior applies on that sampling scale. Parameters without a proper prior are flagged as flat.

Usage

prior_info(fit, format = "table")

Arguments

Details

All priors and constants apply on the sampling scale set by each parameter's link function: a prior for a log-link parameter describes the log of that parameter, and a constant fixes the parameter on that scale (e.g. constant(0) with a log link fixes the parameter to 1 on the native scale). Constants are therefore printed together with their exact native-scale value whenever the two differ.

The provenance of each prior is determined by re-deriving the default priors from the model, formula and data stored in the fit. A user-specified prior that is identical to the bmm default is therefore reported as a default. Coefficients that inherit their prior from a more general class are collapsed into the row of the prior they inherit from.

The one exception is a parameter fixed to a constant in the formula: bmm folds such a constant into the model object at fit time, overwriting the default it replaced, so the original default is no longer recoverable from the fit. Constants named in the formula are therefore always reported as "user", even when the value restates the bmm default.

Re-derivation uses the default_priors stored inside the fit, so later changes to a model's default values do not affect the report of an older fit. It does use the currently installed code that turns those values into prior rows, so if that construction changes between versions, the reported provenance of a fit made with an earlier version may differ from what was actually applied. No check against fit$version$bmm is performed.

Flat priors are flagged because they are improper: Bayes factors via bridge sampling are undefined when any parameter has an improper prior.

Parameters that exist only because the family machinery requires them are omitted: a fixed parameter the model does not declare (the mu brms forces on the custom families of the response-time models, or the mu2 and kappa2 of a two-component brms::mixture() family), and the mixture-weight reference component theta2, which brms lists with a default prior although the sampler holds it at zero. Fixed parameters the model does declare are reported, since the user can estimate them: mu for sdm(), mu1 for the circular mixture models, zr for ddm().

Value

A data.frame of class bmm_prior_info with columns parameter, link, class, coef, group, prior and source. Subsetting the report with [ returns a plain data.frame, as the report-specific printing depends on columns and attributes that subsetting drops.

See Also

default_prior(), parameter_info()

Examples

fit <- bmm(
  bmf(c ~ 0 + set_size, kappa ~ 1),
  data = oberauer_lin_2017,
  model = sdm(resp_error = "dev_rad")
)
prior_info(fit)
prior_info(fit, format = "text")

Rejection Sampling

Description

Performs rejection sampling to generate samples from a target distribution. Each draw can come from its own target: draw i is sampled from f evaluated at the i-th element of every per-draw argument in ..., under the envelope max_f[i].

Usage

rejection_sampling(n, f, max_f, proposal_fun, ...)

Arguments

Value

A numeric vector of length n containing samples from the target distribution.

Examples

target_density <- function(x) brms::dvon_mises(x, mu = 0, kappa = 10)
proposal <- function(n) runif(n, min = -pi, max = pi)
samples <- rejection_sampling(10000, target_density, max_f = target_density(0), proposal)
hist(samples, freq = FALSE)
curve(target_density, col = "red", add = TRUE)

# one location per draw
mu <- rep(c(0, 2), 5000)
samples <- rejection_sampling(
  10000, brms::dvon_mises, max_f = brms::dvon_mises(0, 0, 10), proposal,
  mu = mu, kappa = 10
)
tapply(samples, mu, mean)

Restructure Old bmmfit Objects

Description

Restructure old bmmfit objects to work with the latest bmm version. This function is called internally when applying post-processing methods.

Usage

## S3 method for class 'bmmfit'
restructure(x, ...)

Arguments

Value

A bmmfit object compatible with the latest version of bmm and brms.

Examples

# Load an old bmmfit object
old_fit <- readRDS("bmmfit_old.rds")
new_fit <- restructure(old_fit)

Generic S3 method for reverting any postprocessing of the fitted brm model

Description

Called by update.bmmfit() to automatically revert some of the postprocessing depending on the model type. It will call the appropriate revert_postprocess_brm.* methods based on the list of classes defined in the .model_* functions. For models with several classes listed, it will call the functions in the order they are listed. For example, for the sdm model, the postprocessing involves setting the link function for the c parameter to "log", because it was coded manually in the stan code, but it was specified as "identity" in the brms custom family. However, during the update process, the link function should be set back to "identity". Only use this if you have a specific reason to revert the postprocessing (if otherwise the update method would produce incorrect results).

Usage

revert_postprocess_brm(model, fit, ...)

Arguments

Value

An object of class brmsfit, with any necessary postprocessing applied

Examples

fit <- readRDS("my_saved_fit.rds")
postprocessed_fit <- prostprocess_brm(fit)
reverted_fit <- revert_postprocess_brm(postprocessed_fit)

Observed ROC points from a fitted SDT model

Description

Computes empirical (observed) ROC points from the response-count data used to fit a sdt_rating() or sdt_yn() model. Rating models pool counts across observations within each stimulus type (and optional condition) to produce cumulative hit/false-alarm rates at each confidence threshold; binary models produce one operating point per criterion level (the empirical counterpart of the model-implied points from roc_sdt()).

Usage

roc_observed(fit, conditions = NULL)

Arguments

Value

A data frame of class "bmm_sdt_roc_observed" with columns FA, Hit, and any condition columns. Rating models additionally include the (0,0) and (1,1) endpoints. The attribute model_type is "rating" or "binary"; rating results also carry n_ratings.

See Also

roc_sdt(), plot.bmm_sdt_roc()

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

# One point per base-rate condition, from the response counts
obs <- roc_observed(fit)
obs
plot(roc_sdt(fit), observed = obs)

## End(Not run)

Model-implied ROC curve from a fitted SDT model

Description

Computes the receiver operating characteristic (ROC) curve implied by the posterior distribution of a signal detection model fit with bmm(). ROC curves require a response criterion, so they are defined for sdt_yn() and sdt_rating() only (not the criterion-free sdt_mafc() and sdt_ranking()).

Usage

roc_sdt(
  fit,
  conditions = NULL,
  n_points = 100,
  probs = c(0.025, 0.975),
  criterion_points = NULL,
  ...
)

Arguments

Details

For binary models the ROC is traced analytically from the posterior of d (and sdratio for unequal-variance SDT) over a grid of criterion values. When the criterion varies across conditions but d/sdratio do not (e.g. a base-rate manipulation, as in broeder_schuetz_2009_e3), the several criteria are operating points on a single curve: roc_sdt() returns one smooth curve and attaches the model-implied points (one per criterion level) as a points attribute.

For rating models the K-1 confidence thresholds define K-1 empirical ROC points per posterior draw (returned as the data frame). The smooth model-implied curve is traced over a virtual cut from the posterior of d (and sdratio) and attached as the summary attribute, with the K-1 thresholds attached as the points attribute (labelled t1..t(K-1)) so they fall on the curve. For version = "metad" the smooth curve is the type-1 (d, sdratio) ROC and the threshold points are the model's confidence operating points; they lie inside the curve when the M-ratio is below 1, by the amount of the metacognitive inefficiency.

Value

A data frame of class "bmm_sdt_roc" with columns FA, Hit, .draw, and any condition columns. The object carries a summary attribute (FA, Hit_mean, Hit_lower, Hit_upper) with the smooth model-implied curve, and a points attribute with the model-implied operating points: one per criterion level for binary multi-criteria fits, or the K-1 confidence thresholds (labelled t1..t(K-1)) for rating fits. It also carries the attributes probs, model_class, dist, is_rating and conditions, which the print() and plot() methods read. The Hit_lower and Hit_upper columns in the summary attribute are the pointwise posterior quantiles of the hit rate at each criterion value, plotted at that criterion's posterior-mean false-alarm rate; this band does not include uncertainty in the false-alarm rate and is narrower than a credible band at a fixed false-alarm rate.

See Also

auc_sdt(), roc_observed(), plot.bmm_sdt_roc()

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

roc <- roc_sdt(fit)
roc
plot(roc, observed = roc_observed(fit))
plot(roc, observed = roc_observed(fit), scale = "quantile")

## End(Not run)

Signal Discrimination Model (SDM) by Oberauer (2023)

Description

Signal Discrimination Model (SDM) by Oberauer (2023)

Usage

sdm(resp_error, version = "simple", ...)

sdmSimple(resp_error, version = "simple", ...)

Arguments

Details

see the online article for a detailed description of the model and how to use it. * Domain: Visual working memory

• Task: Continuous reproduction

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

• Citation:

  • Oberauer, K. (2023). Measurement models for visual working memory—A factorial model comparison. Psychological Review, 130(3), 841-852. https://doi.org/10.1037/rev0000328

• Version: simple

• Requirements:

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

• Parameters:

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

• Fixed parameters:

  • mu = 0

• Default parameter links:

  • mu = tan_half; c = log; kappa = log

• Default priors:

  • mu:

    • main: normal(0, 0.5)

    • effects: normal(0, 0.25)

    • sd: exponential(4)

  • kappa:

    • main: student_t(5, 1.75, 0.75)

    • effects: normal(0, 1)

    • sd: exponential(1)

  • c:

    • main: student_t(5, 2, 0.75)

    • effects: normal(0, 1)

    • sd: exponential(1)

Value

An object of class bmmodel

Examples

# simulate data from the model
dat <- data.frame(y = rsdm(n = 1000, c = 4, kappa = 3))

# specify formula
ff <- bmf(
  c ~ 1,
  kappa ~ 1
)

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

Continuous Dual-Process Signal Detection Theory Model

Description

Continuous Dual-Process Signal Detection Theory Model

Usage

sdt_cdp(
  response = "",
  stimulus,
  n_new,
  n_old,
  dist = "normal",
  threshold_type = c("parsimonious", "equidistant", "log_distance"),
  links = NULL,
  ...
)

Arguments

Details

• Domain: Recognition Memory

• Task: Old/New Recognition with Remember/Know Judgments

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

• Citation:

  • Wixted, J. T., & Mickes, L. (2010). A continuous dual-process model of remember/know judgments. Psychological Review, 117(4), 1025-1054. https://doi.org/10.1037/a0020874

• Requirements:

  Provide aggregated data with one row per cell (unique combination of predictors and stimulus class) and one integer count column per response category:

  • new1 ... ⁠new<n_new>⁠ for 'new' judgments

  • ⁠know<k>⁠ and ⁠remember<k>⁠ (and optionally ⁠guess<k>⁠) for 'old' judgments, with k on the unified confidence scale n_new+1 ... n_new+n_old

  • an optional common column prefix is set via response

  • stimulus: 0 = new/lure, 1 = old/target

Use aggregate_sdt_cdp_data() to build these columns from long-format (trial-level) data

• Parameters:

  • dfam: Familiarity sensitivity: target mean on the familiarity axis, in lure-SD units (d' on that axis, where targets and lures share SD 1)

  • drec: Recollection sensitivity: target mean on the recollection axis, in lure-SD units. The target SD is exp(sigmar), so once sigmar is estimated drec is not the d_a the other SDT models report as d

  • criterion: Response bias: old/new boundary on the strength (F+R) axis

  • spacing: Threshold spacing: controls distance between adjacent thresholds (exp(spacing) ensures positive spacing)

  • rcrit: Remember criterion: threshold on the recollection axis

  • sigmar: Log SD of the recollection target distribution, so exp(sigmar) is the SD and 0 means SD = 1

  • rho: Familiarity-recollection correlation on an unconstrained scale; tanh(rho) is the correlation and 0 means independent processes

  • kcrit: Know criterion: threshold on the familiarity axis that splits Know from Guess. Active only when the data include 'guess' counts

• Fixed parameters:

  • sigmar = 0

  • rho = 0

  • kcrit = -100

• Default parameter links:

  • dfam = identity; drec = identity; criterion = identity; spacing = identity; rcrit = identity; sigmar = identity; rho = identity; kcrit = identity

• Default priors:

  • dfam:

    • main: normal(1, 1)

    • effects: normal(0, 0.5)

    • sd: exponential(1)

  • drec:

    • main: normal(1, 1)

    • effects: normal(0, 0.5)

    • sd: exponential(1)

  • criterion:

    • main: normal(0, 1.5)

    • effects: normal(0, 0.5)

    • sd: exponential(2)

  • spacing:

    • main: normal(0, 0.5)

    • effects: normal(0, 0.3)

    • sd: exponential(2)

  • rcrit:

    • main: normal(0, 1)

    • effects: normal(0, 0.5)

    • sd: exponential(2)

  • sigmar:

    • main: normal(0, 0.5)

    • effects: normal(0, 0.3)

    • sd: exponential(2)

  • rho:

    • main: normal(0, 0.5)

    • effects: normal(0, 0.3)

    • sd: exponential(2)

  • kcrit:

    • main: normal(0, 1)

    • effects: normal(0, 0.5)

    • sd: exponential(2)