R/model_sdt_rating.R

In bmm: Easy and Accessible Bayesian Measurement Models Using 'brms'

############################################################################# !
# CONFIDENCE-THRESHOLD PARAMETERIZATIONS                                  ####
############################################################################# !

# Threshold-parameterization registry: the position defines the integer
# thresh_type code passed to Stan -- reordering entries changes the R <-> Stan
# contract (see sdt_make_thresholds_rating in sdt_rating_funs.stan)
.sdt_threshold_types <- c("parsimonious", "equidistant", "log_distance",
                          "log_ratio", "softmax")

# Parameter spec for one threshold parameterization. parsimonious/equidistant
# need a single spacing; log_distance/log_ratio need K-2 interval deltas, named
# in interval order (delta_i lies between thresholds i and i + 1); softmax needs
# spacing plus K-3 allocation deltas.
# Every threshold parameter is a log-scale quantity -- the Stan builders read
# exp(spacing) and exp(delta) -- so its sd default follows sdratio's rate 2
# rather than d's rate 1: at rate 2 an individual's spacing stays within a
# factor of about 2.8 of the group value at the one-sided 95% quantile (90%
# central), at rate 1 within a factor of 8.
.sdt_threshold_parameter_parts <- function(n_ratings, threshold_type) {
  parameters <- list()
  default_priors <- list()
  param_links <- list()
  delta_prior <- list(
    main = "normal(0, 1)", effects = "normal(0, 0.5)", sd = "exponential(2)"
  )

  if (threshold_type %in% c("equidistant", "parsimonious")) {
    parameters$spacing <- paste0(
      "Threshold spacing: controls distance between adjacent thresholds ",
      "(exp(spacing) ensures positive spacing)"
    )
    default_priors$spacing <- list(
      main = "normal(0, 0.5)", effects = "normal(0, 0.3)", sd = "exponential(2)"
    )
    param_links$spacing <- "identity"
  } else if (threshold_type == "log_distance") {
    for (i in seq_len(n_ratings - 2L)) {
      pname <- paste0("delta", i)
      parameters[[pname]] <- glue(
        "Log width of the interval between thresholds {i} and {i + 1}"
      )
      default_priors[[pname]] <- delta_prior
      param_links[[pname]] <- "identity"
    }
  } else if (threshold_type == "log_ratio") {
    labels <- .sdt_log_ratio_delta_labels(n_ratings)
    for (i in seq_along(labels)) {
      pname <- paste0("delta", i)
      parameters[[pname]] <- labels[i]
      default_priors[[pname]] <- delta_prior
      param_links[[pname]] <- "identity"
    }
  } else if (threshold_type == "softmax") {
    parameters$spacing <- paste0(
      "Average log spacing across threshold intervals ",
      "(exp(spacing) is the mean interval size)"
    )
    default_priors$spacing <- list(
      main = "normal(0, 0.5)", effects = "normal(0, 0.3)", sd = "exponential(2)"
    )
    param_links$spacing <- "identity"

    n_deltas <- max(0L, n_ratings - 3L)
    for (i in seq_len(n_deltas)) {
      pname <- paste0("delta", i)
      parameters[[pname]] <- glue(
        "Softmax threshold allocation parameter for interval {i}, relative ",
        "to the last interval"
      )
      default_priors[[pname]] <- delta_prior
      param_links[[pname]] <- "identity"
    }
  }

  nlist(parameters, default_priors, param_links)
}

# The log_ratio deltas are not exchangeable (see .sdt_log_ratio_widths()), so
# each label names the role its interval plays.
.sdt_log_ratio_delta_labels <- function(n_ratings) {
  n_gaps <- n_ratings - 2L
  labels <- character(n_gaps)
  if (n_ratings %% 2L == 0L) {
    m <- n_ratings %/% 2L
    labels[m] <- "Log spread: width of the interval just above the criterion"
    labels[m - 1L] <- paste0(
      "Log ratio of the interval just below the criterion to the spread above it"
    )
    first <- c(rep(m - 1L, max(0L, m - 2L)), rep(m, max(0L, n_gaps - m)))
    rest <- setdiff(seq_len(n_gaps), c(m - 1L, m))
  } else {
    g <- (n_ratings - 1L) %/% 2L
    labels[g] <- "Log spread: width of the middle category, centred on the criterion"
    labels[c(g - 1L, g + 1L)] <- glue(
      "Log ratio of interval {c(g - 1L, g + 1L)} to the middle category"
    )
    rest <- setdiff(seq_len(n_gaps), c(g - 1L, g, g + 1L))
    first <- ifelse(rest < g, g - 1L, g + 1L)
  }
  labels[rest] <- glue("Log ratio of interval {rest} to interval {first}")
  labels
}

.sdt_threshold_delta_names <- function(object) {
  grep("^delta", names(object$parameters), value = TRUE)
}


############################################################################# !
# MODELS                                                                 ####
############################################################################# !

# Per-version specification for the rating SDT lattice. `standard` is a single
# familiarity process; `dpsdt` adds a recollection threshold process (Yonelinas,
# 1994; recall-to-reject of new items after Yonelinas, 2024) whose Ro/Rn are
# fixed near zero by default (recovering standard SDT) and freed through the
# formula; `metad` adds a type-2 metacognitive process
# (Maniscalco & Lau, 2012) parameterized by the log M-ratio. Each entry carries
# the extra parameters plus the Stan
# logmu function name and call the shared multinomial pipeline assembles, so the
# formula, family, and prediction code stay version-agnostic.
.sdt_rating_variants <- list(
  standard = list(
    parameters = list(),
    default_priors = list(),
    links = list(),
    init_ranges = list(),
    fixed_parameters = list(),
    extra_params = character(0),
    logmu_fun = "sdt_rating_logmu",
    logmu_cat_call = "sdt_rating_logmu_cat(cat, thr, d, sdratio, stimulus, dist_type)",
    stan_chunk = "sdt_rating_funs.stan",
    citation = glue(
      "Green, D. M., & Swets, J. A. (1966). Signal detection theory ",
      "and psychophysics. Wiley."
    )
  ),
  dpsdt = list(
    parameters = list(
      Ro = glue("Recollection of old items: probability inv_logit(Ro) that an ",
                "old item is recollected as old, loading the most-confident ",
                "'signal' category"),
      Rn = glue("Recollection of new items: probability inv_logit(Rn) that a ",
                "new item is recall-rejected, loading the most-confident ",
                "'noise' category")
    ),
    # Ro and Rn are probabilities on a logit scale, which is the mixing-weight
    # case the package's sd rule assigns rate 1: an individual then sits
    # between p = .12 and p = .89 around a group value of .5 at the 95%
    # quantile, where rate 2 would hold them inside [.26, .74].
    default_priors = list(
      Ro = list(main = "normal(0, 1)", effects = "normal(0, 0.5)", sd = "exponential(1)"),
      Rn = list(main = "normal(0, 1)", effects = "normal(0, 0.5)", sd = "exponential(1)")
    ),
    links = list(Ro = "identity", Rn = "identity"),
    init_ranges = list(Ro = c(-0.5, 0.5), Rn = c(-0.5, 0.5)),
    # inv_logit(-100) is numerically 0: recollection off unless freed by formula
    fixed_parameters = list(Ro = -100, Rn = -100),
    extra_params = c("Ro", "Rn"),
    logmu_fun = "sdt_dpsdt_logmu",
    logmu_cat_call = "sdt_dpsdt_logmu_cat(cat, thr, d, sdratio, stimulus, dist_type, Ro, Rn)",
    stan_chunk = "sdt_dpsdt_funs.stan",
    citation = c(
      glue(
        "Yonelinas, A. P. (1994). Receiver-operating characteristics in ",
        "recognition memory: Evidence for a dual-process model. Journal of ",
        "Experimental Psychology: Learning, Memory, and Cognition, 20(6), ",
        "1341-1354. https://doi.org/10.1037/0278-7393.20.6.1341"
      ),
      glue(
        "Yonelinas, A. P. (2024). The role of recollection and familiarity in ",
        "visual working memory: A mixture of threshold and signal detection ",
        "processes. Psychological Review, 131(2), 321-348. ",
        "https://doi.org/10.1037/rev0000432"
      )
    )
  ),
  metad = list(
    parameters = list(
      logmratio = glue("Log M-ratio, log(meta-d/d): metacognitive efficiency. ",
                       "0 is ideal metacognition (meta-d = d), negative is ",
                       "inefficiency, positive is hyper-efficiency. meta-d is ",
                       "recovered as exp(logmratio) * d, on the same scale as d")
    ),
    # Estimating log(meta-d/d) rather than meta-d directly keeps meta-d
    # positive, regularizes it toward d, and anchors the field-standard
    # metacognitive-efficiency measure (M-ratio) at the ideal point of 0
    # (Maniscalco & Lau, 2012; Fleming, 2017). Both sensitivities share one
    # scale (d', or d_a when sdratio is estimated), so the ratio is invariant to
    # sdratio.
    # logmratio is a log ratio, so its sd rate follows sdratio's 2 rather than
    # d's 1: an individual's M-ratio then stays within a factor of about 2.8
    # of the group value at the 95% quantile.
    default_priors = list(
      logmratio = list(main = "normal(0, 0.5)", effects = "normal(0, 0.3)", sd = "exponential(2)")
    ),
    links = list(logmratio = "identity"),
    init_ranges = list(logmratio = c(-0.3, 0.3)),
    fixed_parameters = list(),
    extra_params = "logmratio",
    logmu_fun = "sdt_metad_logmu",
    logmu_cat_call = "sdt_metad_logmu_cat(cat, thr, d, exp(logmratio) * d, sdratio, stimulus, dist_type)",
    stan_chunk = "sdt_metad_funs.stan",
    citation = glue(
      "Maniscalco, B., & Lau, H. (2012). A signal detection theoretic ",
      "approach for estimating metacognitive sensitivity from confidence ",
      "ratings. Consciousness and Cognition, 21(1), 422-430. ",
      "https://doi.org/10.1016/j.concog.2011.09.021"
    )
  )
)

.model_sdt_rating <- function(response = NULL, stimulus = NULL,
                              dist = "normal",
                              threshold_type = "parsimonious",
                              version = "standard",
                              links = NULL, call = NULL, ...) {
  # one count column per category, so the columns fix the number of categories
  n_ratings <- length(response)

  parameters <- list(
    d = paste0(
      "Sensitivity: d' under equal variance (the default). When sdratio is ",
      "estimated, d is d_a, the distance between the signal and noise ",
      "distributions in units of their root-mean-square SD"
    ),
    criterion = paste0(
      "Decision criterion, measured from the mean of the noise distribution in ",
      "noise-SD units: the middle threshold (the old/new boundary) for an even ",
      "number of categories, the centre of the middle category for an odd number"
    ),
    bias = paste0(
      "Alternative to criterion: the same middle threshold measured from the ",
      "midpoint between the noise and signal means (0 = neutral, positive = ",
      "conservative). Use one of the two"
    )
  )
  # sd rates as in sdt_yn, so the same subjects shrink the same way whichever
  # SDT model they are fitted with: rate 1 for the sensitivity d, whose
  # between-subject SD is ~0.6 on broeder_schuetz_2009_e3 and 0.30 on
  # meyer_grant_jakob_2025, and rate 2 for bias (~0.15 there) and for the
  # log-scale sdratio; criterion carries half of d's spread on top of bias's
  # and takes rate 1, as in sdt_yn, whose comment gives the measurement.
  default_priors <- list(
    d = list(main = "normal(1, 1)", effects = "normal(0, 0.5)", sd = "exponential(1)"),
    criterion = list(main = "normal(0.5, 1.5)", effects = "normal(0, 0.5)", sd = "exponential(1)"),
    bias = list(main = "normal(0, 1.5)", effects = "normal(0, 0.5)", sd = "exponential(2)")
  )
  param_links <- list(d = "identity", criterion = "identity", bias = "identity")

  threshold_parts <- .sdt_threshold_parameter_parts(n_ratings, threshold_type)
  parameters <- c(parameters, threshold_parts$parameters)
  default_priors <- c(default_priors, threshold_parts$default_priors)
  param_links <- c(param_links, threshold_parts$param_links)

  parameters$sdratio <- paste0(
    "Log SD ratio: the log of the signal-to-noise standard deviation ratio, ",
    "so exp(sdratio) is the ratio itself and 0 means equal variance"
  )
  # Identical to sdt_yn's, whose comment gives the empirical grounds.
  default_priors$sdratio <- list(
    main = "normal(0, 0.3)", effects = "normal(0, 0.3)", sd = "exponential(2)"
  )
  param_links$sdratio <- "identity"

  # Bounded starting values for every estimated parameter. Without them the
  # flexible threshold types start from brms' wide random init, where adjacent
  # thresholds can land in the same distribution tail and the category mass
  # underflows to log_diff_exp(-Inf, -Inf) = NaN, so every chain rejects its
  # initial value.
  init_ranges <- list(
    d = c(0.5, 1.5), criterion = c(0, 1), bias = c(-0.5, 0.5),
    sdratio = c(-0.3, 0.3)
  )
  for (p in names(threshold_parts$parameters)) {
    init_ranges[[p]] <- if (p == "spacing") c(-0.7, -0.2) else c(-0.5, 0.2)
  }

  variant <- .sdt_rating_variants[[version]]
  parameters <- c(parameters, variant$parameters)
  default_priors <- c(default_priors, variant$default_priors)
  param_links <- c(param_links, variant$links)
  init_ranges <- c(init_ranges, variant$init_ranges)
  fixed_parameters <- c(list(sdratio = 0), variant$fixed_parameters)

  requirements <- glue(
    "Provide pre-aggregated data with the following columns:", "\n\n",
    "  - Response counts: one column per rating category (K columns)", "\n",
    "  - Stimulus type (stimulus): 0 = noise, 1 = signal", "\n",
    "  Categories should be ordered: 1 = 'definitely noise' to ",
    "K = 'definitely signal'"
  )

  out <- structure(
    list(
      resp_vars = nlist(response),
      other_vars = nlist(stimulus, dist, n_ratings, threshold_type),
      domain = "Perception & Recognition Memory",
      task = "Signal/Noise or Old/New Recognition",
      name = "Signal Detection Theory (Confidence Rating)",
      citation = variant$citation,
      version = version,
      requirements = requirements,
      parameters = parameters,
      links = param_links,
      fixed_parameters = fixed_parameters,
      alternative_parameters = list(c("criterion", "bias")),
      default_priors = default_priors,
      init_ranges = init_ranges
    ),
    class = c("bmmodel", "sdt", "sdt_rating", paste0("sdt_rating_", version)),
    call = call
  )
  set_links(out, links)
}

# `sdratio` is fixed at 0 and read as a log SD ratio (the Stan code takes its
# exp()), and every threshold parameter is read the same way, so neither
# survives a change of link: it would rescale the fixed value and stack a
# second transformation on the exp(). `d`, `criterion` and `bias` fix nothing,
# so their links stay settable, as in sdt_yn.
#' @exportS3Method
settable_links.sdt_rating <- function(model) {
  c("d", "criterion", "bias")
}

# the links are applied by substituting the inverse link into the multinomial
# formula (.sdt_rating_logmu_args), so the model can honour inv_link()'s links,
# not every link a brms family can emit
#' @exportS3Method
settable_link_functions.sdt_rating <- function(model) {
  eval(formals(inv_link)$link)
}


#' @title Confidence Rating Signal Detection Theory Model
#' @name sdt_rating
#' @details Three versions share the confidence-rating response interface,
#'   selected with `version`:
#'
#' #### Version: `standard` (default)
#' `r model_docs(.model_sdt_rating(version = "standard"))`
#'
#' #### Version: `dpsdt`
#' Dual-process SDT (Yonelinas, 1994): a familiarity SDT process plus an
#' all-or-none recollection process that loads the most-confident category.
#' `Ro` is recollection of old items (loads the most-confident "signal"
#' category), the only recollection term of the classic model; `Rn` is
#' recall-to-reject of new items (loads the most-confident "noise" category),
#' the two-sided extension that Yonelinas (2024) applies to visual working
#' memory. `inv_logit(Ro)`/`inv_logit(Rn)` are the recollection
#' probabilities. Both are fixed off by default
#' (recovering `standard`); add `Ro ~ 1` for the one-sided model and
#' `Ro ~ 1, Rn ~ 1` for the two-sided model. `d` describes the familiarity
#' distributions only -- the observed ROC is a mixture of familiarity and
#' recollection, so it is not the discriminability of that mixture.
#' `summary()` reports `Ro` and `Rn` on the logit scale; [latent_sdt()]
#' returns them as probabilities (attribute `extra`).
#' Two cautions. Freeing `sdratio` alongside `Ro` is weakly identified from
#' a single ROC: in simulation the two are correlated at about -0.8 in the
#' posterior and `Ro` is pulled down while `exp(sdratio)` is pulled above 1;
#' keep the familiarity process equal-variance unless the design separates
#' them. And a recollection probability near zero is reported as the tail of
#' its prior: on the probability scale the interval cannot include 0, so the
#' test of "no recollection" is a comparison with the fit that leaves the
#' parameter fixed off, not the interval.
#' `r model_docs(.model_sdt_rating(version = "dpsdt"))`
#'
#' #### Version: `metad`
#' Meta-d' (Maniscalco & Lau, 2012): a type-2 metacognitive sensitivity governs
#' confidence-threshold placement, with the total "old"/"new" response rates
#' held to what type-1 `d` implies. Rather than estimating meta-d'
#' directly, the model estimates `logmratio`, the log M-ratio
#' \eqn{\log(\mathrm{meta\text{-}d'}/d')}, and recovers meta-d' as
#' `exp(logmratio) * d`. The M-ratio is the field-standard measure of
#' metacognitive efficiency (Maniscalco & Lau, 2012; Fleming, 2017): estimating
#' it on the log scale keeps meta-d' positive, regularizes it toward `d`,
#' and anchors the ideal point (meta-d' = `d`, perfect metacognition) at
#' `logmratio = 0`, which recovers `standard`. Type-1 and type-2 sensitivity
#' share one scale (\eqn{d'}, or \eqn{d_a} when `sdratio` is estimated), so the
#' M-ratio is unaffected by `sdratio`.
#' The type-1 boundary is the middle threshold (`criterion` or `bias`), so this
#' version needs an even number of rating categories: with an odd number the
#' middle category straddles the boundary (see "Where `criterion` and `bias`
#' sit").
#' Every threshold keeps its distance from the midpoint between the means when
#' `d` is replaced by meta-d', as the HMeta-d toolbox's model equations do
#' (Fleming, 2017, Appendix). Whether the formula estimates `criterion` or
#' `bias` changes how the type-1 boundary is parameterized, not this
#' constraint.
#' Maximum-likelihood meta-d' instead constrains meta-c' = c' (Maniscalco &
#' Lau, 2014), so its estimates differ from bmm's when the criterion is far
#' from the midpoint.
#' Extract the M-ratio posterior with [mratio()].
#' `r model_docs(.model_sdt_rating(version = "metad"))`
#'
#' By default, the model assumes equal variance (sdratio fixed to 0). To
#' estimate unequal variance, add `sdratio ~ 1` (or `sdratio ~ predictors`)
#' to the formula.
#'
#' @section Sensitivity is on the same scale as [sdt_yn()]:
#' `d` is \eqn{d'} whenever `sdratio` stays fixed at 0, which is every fit
#' that does not give `sdratio` a formula. With `sdratio` estimated, `d` is the
#' balanced index \eqn{d_a} that [sdt_yn()] reports: the separation between the
#' signal and noise distributions divided by the root-mean-square of their SDs.
#' Unlike the noise-standardized \eqn{d'}, it remains comparable across
#' conditions that differ in `sdratio`; see the sensitivity section of
#' [sdt_yn()] for the reasoning, for how far the two indices lie apart, and for
#' the caveat that under the Gumbel distributions \eqn{d_a} is not the
#' AUC-equivalent index once `sdratio` is estimated.
#'
#' The criterion and the confidence thresholds are **not** rescaled. They stay
#' on the noise-standardized axis, so under unequal variance they and `d` are
#' in different units.
#'
#' The `dpsdt` and `metad` versions inherit the same convention: `d` there
#' describes the familiarity (type-1) distributions, and meta-d' is scaled the
#' same way, which leaves the M-ratio invariant to `sdratio`.
#'
#' @section Where `criterion` and `bias` sit:
#' Every `threshold_type` places the criterion the same way. With an even
#' number of categories K it is the middle threshold, the boundary between the
#' K/2 "noise" categories and the K/2 "signal" categories. With an odd number
#' there is no such boundary -- the middle category straddles it -- so the
#' criterion is the centre of that category, and the two thresholds around it
#' sit half an interval below and above. A shift in the criterion therefore
#' moves the whole threshold set, and the threshold parameters (`spacing`,
#' `delta`) describe its shape around that point.
#'
#' As in [sdt_yn()], the formula estimates that point either as `criterion`,
#' measured from the mean of the noise distribution, or as `bias`, measured
#' from the midpoint between the noise and signal means (`bias = criterion -
#' d * sqrt((1 + exp(sdratio)^2) / 2) / 2`); a formula that names neither
#' estimates `criterion ~ 1`. The two differ in what a linear formula holds
#' constant when `d` or `sdratio` varies; the section "Choosing between
#' `criterion` and `bias`" of [sdt_yn()] says which to choose. Threshold locations reported by
#' [sdt_thresholds()] and [latent_sdt()], and the `thresholds` taken by
#' [rsdt_rating()] and its siblings, are positions measured from the noise
#' mean.
#'
#' @section Identifying `sdratio`:
#' Unlike [sdt_yn()], a rating design identifies `sdratio` from a single
#' condition: the K - 1 thresholds give K - 1 operating points per condition,
#' which trace the ROC; under `dist = "normal"` its z-transform is a line with
#' slope `1 / exp(sdratio)`. `sdratio ~ 1` on a one-condition dataset is
#' identified from K = 3 categories, where it uses every degree of freedom
#' (the fit is saturated), and leaves the zROC testable from K = 4; the
#' identification caveats on the [sdt_yn()] page do not carry over.
#' @param response A character vector of K column names containing response
#'   counts per rating category, ordered from "definitely noise" to
#'   "definitely signal".
#' @param stimulus The name of the variable coding the stimulus type.
#'   Must be coded as 0 (noise/new) and 1 (signal/old).
#' @param dist The distribution assumed for the latent evidence, given here by
#'   its cumulative distribution function. One of:
#'   \itemize{
#'     \item "normal" (default): Gaussian SDT, \eqn{\Phi(x)}
#'     \item "gumbel_min": smallest-extreme-value SDT,
#'       \eqn{1 - \exp(-\exp(x))} (complementary log-log)
#'     \item "gumbel_max": largest-extreme-value SDT, \eqn{\exp(-\exp(-x))}
#'       (log-log, as in \code{evd::pgumbel})
#'     \item "logistic": logistic SDT, \eqn{1 / (1 + \exp(-x))}
#'   }
#' @param threshold_type Character. Threshold parameterization:
#'   \itemize{
#'     \item "parsimonious" (default): 2 parameters (criterion + spacing).
#'       Thresholds follow logit-spaced canonical positions (Selker et al.,
#'       2019).
#'     \item "equidistant": 2 parameters (criterion + spacing). Thresholds
#'       are equally spaced, exp(spacing) apart.
#'     \item "log_distance": K-2 parameters. `delta<i>` is the log width of
#'       the interval between thresholds i and i + 1, so the intervals are
#'       free and the ordering is guaranteed (Paulewicz & Blaut, 2022).
#'     \item "log_ratio": K-2 parameters (Paulewicz & Blaut, 2022; the odd-K
#'       form is bmm's). One interval is the spread, `exp(delta)`: for even K
#'       the interval just above `criterion`, for odd K the middle category.
#'       The first interval on the other side (even K) or on each side (odd K)
#'       is a ratio times the spread, and every further interval a ratio
#'       times the first interval on its side. The deltas are therefore not
#'       exchangeable; `model$parameters` names the role of each.
#'     \item "softmax": K-2 parameters. A shared spacing parameter sets the
#'       mean interval width, exp(spacing), while K-3 delta parameters share
#'       the total width out over the intervals through a softmax (each delta
#'       is the log ratio of its interval to the last one).
#'   }
#' @param version Character. The latent-process version of the rating model.
#'   One of `"standard"` (default, a single familiarity SDT process),
#'   `"dpsdt"` (dual-process SDT with recollection parameters `Ro`/`Rn`), or
#'   `"metad"` (meta-d' with a metacognitive efficiency parameter `logmratio`,
#'   the log M-ratio). All three use the same confidence-rating response
#'   interface. See Details.
#' @param links A named list of link functions for the parameters, e.g.
#'   `links = list(d = "log")`. Only `d`, `criterion` and `bias` can be set, to
#'   `"identity"`, `"log"`, `"softplus"`, `"logit"` or `"probit"`. `sdratio`
#'   and the threshold parameters keep their identity links, because the model
#'   reads each of them through `exp()` and fixes `sdratio` at 0 for equal
#'   variance. The same holds for `Ro` and `Rn`, read through `inv_logit()`
#'   and fixed off by default, and for `logmratio`, read through `exp()`.
#' @param ... used internally for testing, ignore it
#' @return An object of class `bmmodel`
#' @references
#' Yonelinas, A. P. (1994). Receiver-operating characteristics in recognition
#'   memory: Evidence for a dual-process model. \emph{Journal of Experimental
#'   Psychology: Learning, Memory, and Cognition}, \emph{20}(6), 1341--1354.
#'
#' Yonelinas, A. P. (2024). The role of recollection and familiarity in visual
#'   working memory: A mixture of threshold and signal detection processes.
#'   \emph{Psychological Review}, \emph{131}(2), 321--348.
#'   \doi{10.1037/rev0000432}
#'
#' Maniscalco, B., & Lau, H. (2012). A signal detection theoretic approach for
#'   estimating metacognitive sensitivity from confidence ratings.
#'   \emph{Consciousness and Cognition}, \emph{21}(1), 422--430.
#'
#' 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.), \emph{The
#'   cognitive neuroscience of metacognition} (pp. 25--66). Springer.
#'
#' Fleming, S. M. (2017). HMeta-d: hierarchical Bayesian estimation of
#'   metacognitive efficiency from confidence ratings. \emph{Neuroscience of
#'   Consciousness}, \emph{2017}(1), nix007. \doi{10.1093/nc/nix007}
#'
#' Green, D. M., & Swets, J. A. (1966). \emph{Signal detection theory and
#'   psychophysics}. Wiley.
#'
#' Selker, R., van den Bergh, D., Criss, A. H., & Wagenmakers, E.-J. (2019).
#'   Parsimonious estimation of signal detection models from confidence ratings.
#'   \emph{Behavior Research Methods}, \emph{51}(5), 1953--1967.
#'   \doi{10.3758/s13428-019-01231-3}
#'
#' Paulewicz, B., & Blaut, A. (2022). The general causal cumulative model of
#'   ordinal response. \emph{PsyArXiv}. \doi{10.31234/osf.io/e7a3x}
#' @keywords bmmodel
#' @export
#' @examples
#' \dontrun{
#' # EV-SDT rating model
#' dat <- expand.grid(id = 1:20, stimulus = c(0L, 1L))
#' dat <- cbind(dat, rsdt_rating(nrow(dat), 200, dat$stimulus,
#'                               d = 1.5, thresholds = c(0.25, 0.75, 1.25)))
#'
#' model <- sdt_rating(
#'   response = c("r1", "r2", "r3", "r4"),
#'   stimulus = "stimulus"
#' )
#'
#' fit <- bmm(
#'   formula = bmf(d ~ 1, criterion ~ 1, spacing ~ 1),
#'   data = dat,
#'   model = model,
#'   cores = 4,
#'   backend = "cmdstanr"
#' )
#'
#' # UV-SDT: add sdratio to the formula
#' fit_uv <- bmm(
#'   formula = bmf(d ~ 1, criterion ~ 1, spacing ~ 1, sdratio ~ 1),
#'   data = dat,
#'   model = model,
#'   cores = 4,
#'   backend = "cmdstanr"
#' )
#'
#' # Dual-process SDT: free recollection of old items (one-sided) or both
#' # old and new items (two-sided) via the formula
#' model_dp <- sdt_rating(c("r1", "r2", "r3", "r4"), "stimulus", version = "dpsdt")
#' fit_dp <- bmm(bmf(d ~ 1, criterion ~ 1, spacing ~ 1, Ro ~ 1, Rn ~ 1),
#'               data = dat, model = model_dp, backend = "cmdstanr")
#'
#' # Meta-d': estimate metacognitive efficiency (log M-ratio)
#' model_md <- sdt_rating(c("r1", "r2", "r3", "r4"), "stimulus", version = "metad")
#' fit_md <- bmm(bmf(d ~ 1, criterion ~ 1, spacing ~ 1, logmratio ~ 1),
#'               data = dat, model = model_md, backend = "cmdstanr")
#' mratio(fit_md) # posterior M-ratio (meta-d'/d')
#' }
sdt_rating <- function(response, stimulus,
                       dist = c("normal", "gumbel_min", "gumbel_max", "logistic"),
                       threshold_type = c("parsimonious", "equidistant",
                                          "log_distance", "log_ratio",
                                          "softmax"),
                       version = c("standard", "dpsdt", "metad"),
                       links = NULL, ...) {
  call <- match.call()
  stop_missing_args()
  dist <- match.arg(dist)
  threshold_type <- match.arg(threshold_type)
  version <- match.arg(version)

  stopif(length(response) <= 2,
         "response must name more than 2 rating-count columns, one per category")
  # each column name becomes part of a brms parameter name (mu<column>)
  unusable <- response[grepl("[._]", response)]
  stopif(length(unusable) > 0,
         "Response column names must not contain '.' or '_', because brms \\
         builds a parameter name from each of them. Rename {collapse_comma(unusable)}")
  stopif(threshold_type == "log_ratio" && length(response) < 4,
         "log_ratio thresholds require at least 4 rating categories (the anchor ratio needs an interval above the criterion)")
  stopif(version == "metad" && length(response) %% 2L != 0L,
         "version = 'metad' needs an even number of rating categories: \\
         meta-d' splits the scale at the old/new boundary, which an odd \\
         number of categories does not have (the middle category straddles \\
         it). response names {length(response)} columns")

  .model_sdt_rating(response = response, stimulus = stimulus,
                    dist = dist, threshold_type = threshold_type,
                    version = version,
                    links = links, call = call, ...)
}


############################################################################# !
# CHECK_DATA S3 METHODS                                                  ####
############################################################################# !

#' @export
check_data.sdt_rating <- function(model, data, formula) {
  stim_var <- model$other_vars$stimulus
  data[[stim_var]] <- .validate_sdt_stimulus(data, stim_var)

  resp_cols <- model$resp_vars$response
  .validate_sdt_count_cols(data, resp_cols)

  Y <- as.matrix(data[resp_cols])
  stopif(anyNA(Y), "Response columns must not contain NA counts")
  stopif(any(rowSums(Y) <= 0),
         "Row sums of response columns must be positive (no empty rows)")

  reserved <- intersect(c("Y", "nTrials"), colnames(data))
  warnif(length(reserved) > 0,
         "Column(s) {collapse_comma(reserved)} in your data are reserved by \\
         {model$name} and will be overwritten")
  data <- data[!colnames(data) %in% resp_cols]
  data$Y <- Y
  data$nTrials <- rowSums(Y)

  NextMethod("check_data")
}


############################################################################# !
# MULTINOMIAL FORMULA & FAMILY CONSTRUCTION                              ####
############################################################################# !

# sdt_rating uses brms' native multinomial family: each category's logit is set
# to log(p_k^SDT) so softmax recovers the SDT probabilities exactly. The SDT math
# (thresholds + noise CDF) is computed by a per-model Stan function for fitting
# and by the exported R companion sdt_rating_logmu() for posterior_predict/epred
# (brms evaluates the non-linear formula in R for prediction, looking the
# function up on the search path). log_lik and posterior_predict therefore come
# from brms — proper joint multinomial draws — while the four distributions and
# five threshold parameterizations stay supported.

# The dispatch pieces (kernel names, extra parameters, Stan chunk) are pure
# derived state: look them up in the variant registry instead of storing them
# in the model object. Legacy model objects from before the version machinery
# carry version = "rating" or "NA" and resolve to the standard variant.
.sdt_rating_variant <- function(model) {
  .sdt_rating_variants[[model$version]] %||% .sdt_rating_variants$standard
}

# cat, K, dist, threshold type and anchor travel as integer literals so the
# same non-linear call resolves against both the generated Stan function and
# the R companion. spacing is the parameter when present, otherwise the literal
# 0 (threshold types without spacing ignore it). `anchor` is the parameter the
# formula estimates for the middle threshold, `criterion` (1) or `bias` (2).
.sdt_rating_logmu_args <- function(model, anchor) {
  has_spacing <- "spacing" %in% names(model$parameters)
  # d and criterion reach the kernel on the natural scale, so a non-identity
  # link is inverted here, as apply_links() does for m3
  c(model$other_vars$n_ratings, .sdt_dist_id(model$other_vars$dist),
    match(model$other_vars$threshold_type, .sdt_threshold_types),
    match(anchor, c("criterion", "bias")),
    deparse(inv_link("d", model$links$d)),
    deparse(inv_link(anchor, model$links[[anchor]])),
    if (has_spacing) "spacing" else "0", "sdratio",
    .sdt_rating_variant(model)$extra_params,
    model$other_vars$stimulus, .sdt_threshold_delta_names(model))
}

# Stan function generated per model: the delta arity is fixed here, while cat,
# K, dist, and threshold type arrive as the integer arguments described above.
.sdt_rating_logmu_stan <- function(model) {
  variant <- .sdt_rating_variant(model)
  delta_names <- .sdt_threshold_delta_names(model)
  nd <- length(delta_names)
  extra <- variant$extra_params

  signature <- paste(
    c("int cat", "int K", "int dist_type", "int thresh_type", "int anchor",
      "real d", "real criterion", "real spacing", "real sdratio",
      if (length(extra)) paste("real", extra),
      "real stimulus", if (nd) paste("real", delta_names)),
    collapse = ", "
  )
  delta_decl <- if (nd == 0L) {
    "  array[0] real deltas;\n"
  } else {
    paste0("  array[", nd, "] real deltas;\n",
           paste(sprintf("  deltas[%d] = %s;", seq_len(nd), delta_names),
                 collapse = "\n"), "\n")
  }

  paste0(
    "real ", variant$logmu_fun, "(", signature, ") {\n",
    delta_decl,
    "  // a bias is measured from the midpoint, the thresholds from the noise mean\n",
    "  real crit = anchor == 2 ? criterion + d * sdt_rms_scale(exp(sdratio)) / 2.0 : criterion;\n",
    "  vector[K - 1] thr = sdt_make_thresholds_rating(crit, spacing, deltas, K, thresh_type);\n",
    "  return ", variant$logmu_cat_call, ";\n",
    "}\n"
  )
}


############################################################################# !
# Convert bmmformula to brmsformula methods                              ####
############################################################################# !

# Base multinomial brmsformula: Y | trials(nTrials) carries the counts and each
# category gets a non-linear mu calling the version's logmu function. The user's
# parameter formulas (d, criterion, spacing/deltas, sdratio) are added afterwards
# by bmf2bf.bmmodel, so they are deliberately not added here.
#' @export
bmf2bf.sdt_rating <- function(model, formula) {
  resp_cats <- model$resp_vars$response
  fun <- .sdt_rating_variant(model)$logmu_fun
  args <- paste(.sdt_rating_logmu_args(model, .sdt_anchor_parameter(formula)),
                collapse = ", ")

  bform <- brms::bf(
    glue("Y | trials(nTrials) ~ {fun}(1, {args})"),
    nl = TRUE
  )
  for (k in seq_along(resp_cats)[-1]) {
    bform <- bform + brms::nlf(stats::as.formula(
      glue("mu{resp_cats[k]} ~ {fun}({k}, {args})")
    ))
  }
  bform
}


############################################################################# !
# CONFIGURE_MODEL S3 METHODS                                             ####
############################################################################# !

#' @export
configure_model.sdt_rating <- function(model, data, formula) {
  resp_cats <- model$resp_vars$response

  formula <- bmf2bf(model, formula)
  formula$family <- brms::multinomial(refcat = NA)
  formula$family$cats <- resp_cats
  formula$family$dpars <- paste0("mu", resp_cats)

  sc_path <- system.file("stan_chunks", package = "bmm")
  chunks <- unique(c("sdt_dist_funs.stan", "sdt_rating_funs.stan",
                     .sdt_rating_variant(model)$stan_chunk))
  stan_funs <- paste(
    paste(vapply(chunks, function(f) read_lines2(paste0(sc_path, "/", f)),
                 character(1)), collapse = "\n"),
    .sdt_rating_logmu_stan(model),
    sep = "\n"
  )
  stanvars <- brms::stanvar(scode = stan_funs, block = "functions")

  nlist(formula, data, stanvars)
}


#' @title Rating SDT category log-probability (multinomial logit)
#' @description R companion to the Stan `sdt_rating_logmu` function. It returns
#'   `log(p_k)` for rating category `cat`, which the [sdt_rating()] multinomial
#'   formula uses as the category logit (so `softmax` recovers the SDT category
#'   probabilities). `brms` evaluates the non-linear formula in R for
#'   `posterior_predict()` and `posterior_epred()`, so this function must be on
#'   the search path; it is exported for that reason and is not called directly.
#' @param cat Integer rating category index.
#' @param K Integer number of rating categories.
#' @param dist Integer noise-distribution id (see the `.sdt_dists` registry).
#' @param thresh Integer threshold-parameterization id.
#' @param anchor Integer: 1 when `criterion` is measured from the noise mean, 2
#'   when it is the `bias` parameter, measured from the midpoint between the
#'   means.
#' @param d,criterion,spacing,sdratio Model parameters (draws-by-observation
#'   matrices supplied by brms). `d` is d', or d_a when sdratio is not 0;
#'   `spacing` is `0` for threshold types without it.
#' @param Ro,Rn Recollection parameters of the `dpsdt` version on the logit
#'   (linear-predictor) scale; `inv_logit(Ro)`/`inv_logit(Rn)` are the
#'   recollection and recall-to-reject probabilities.
#' @param logmratio Log M-ratio for the `metad` version; meta-d' is recovered
#'   as `exp(logmratio) * d`, on the same scale as `d`.
#' @param stimulus Stimulus covariate (0 = noise, 1 = signal).
#' @param ... Threshold `delta` parameters, when the threshold type uses them.
#' @return `log(p_cat)`, matching the shape of `d`.
#' @keywords internal
#' @export
sdt_rating_logmu <- function(cat, K, dist, thresh, anchor, d, criterion, spacing,
                             sdratio, stimulus, ...) {
  dist_name <- .sdt_dist_names[dist]
  thresh_name <- .sdt_threshold_types[thresh]

  shape <- dim(d)
  d <- as.vector(d)
  n <- length(d)
  criterion <- rep_len(as.vector(criterion), n)
  spacing <- rep_len(as.vector(spacing), n)
  sdratio <- rep_len(as.vector(sdratio), n)
  stimulus <- rep_len(as.vector(stimulus), n)
  deltas <- if (...length() > 0L) {
    do.call(cbind, lapply(list(...), function(x) rep_len(as.vector(x), n)))
  }

  if (anchor == 2L) criterion <- .sdt_bias_to_criterion(criterion, d, exp(sdratio))
  thr <- .sdt_make_thresholds(criterion, K, thresh_name, spacing, deltas)
  out <- rbind(.sdt_category_log_probs(rbind(thr), d, exp(sdratio),
                                       stimulus, dist_name))[, cat]

  if (!is.null(shape)) dim(out) <- shape
  out
}


#' @rdname sdt_rating_logmu
#' @export
sdt_dpsdt_logmu <- function(cat, K, dist, thresh, anchor, d, criterion, spacing,
                            sdratio, Ro, Rn, stimulus, ...) {
  dist_name <- .sdt_dist_names[dist]
  thresh_name <- .sdt_threshold_types[thresh]

  shape <- dim(d)
  d <- as.vector(d)
  n <- length(d)
  criterion <- rep_len(as.vector(criterion), n)
  spacing <- rep_len(as.vector(spacing), n)
  sdratio <- rep_len(as.vector(sdratio), n)
  Ro <- rep_len(as.vector(Ro), n)
  Rn <- rep_len(as.vector(Rn), n)
  stimulus <- rep_len(as.vector(stimulus), n)
  deltas <- if (...length() > 0L) {
    do.call(cbind, lapply(list(...), function(x) rep_len(as.vector(x), n)))
  }

  if (anchor == 2L) criterion <- .sdt_bias_to_criterion(criterion, d, exp(sdratio))
  thr <- .sdt_make_thresholds(criterion, K, thresh_name, spacing, deltas)
  out <- rbind(.sdt_dpsdt_category_log_probs(rbind(thr), d, exp(sdratio),
                                             stimulus, dist_name, Ro, Rn))[, cat]

  if (!is.null(shape)) dim(out) <- shape
  out
}


#' @rdname sdt_rating_logmu
#' @export
sdt_metad_logmu <- function(cat, K, dist, thresh, anchor, d, criterion, spacing,
                            sdratio, logmratio, stimulus, ...) {
  dist_name <- .sdt_dist_names[dist]
  thresh_name <- .sdt_threshold_types[thresh]

  shape <- dim(d)
  d <- as.vector(d)
  n <- length(d)
  criterion <- rep_len(as.vector(criterion), n)
  spacing <- rep_len(as.vector(spacing), n)
  sdratio <- rep_len(as.vector(sdratio), n)
  # meta-d is derived from the estimated log M-ratio, mirroring the Stan call
  metad <- exp(rep_len(as.vector(logmratio), n)) * d
  stimulus <- rep_len(as.vector(stimulus), n)
  deltas <- if (...length() > 0L) {
    do.call(cbind, lapply(list(...), function(x) rep_len(as.vector(x), n)))
  }

  if (anchor == 2L) criterion <- .sdt_bias_to_criterion(criterion, d, exp(sdratio))
  thr <- .sdt_make_thresholds(criterion, K, thresh_name, spacing, deltas)
  out <- rbind(.sdt_metad_category_log_probs(rbind(thr), d, exp(sdratio),
                                             stimulus, dist_name, metad))[, cat]

  if (!is.null(shape)) dim(out) <- shape
  out
}