R/model_sdt_ranking.R

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

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

.model_sdt_ranking <- function(response = NULL, m = NULL,
                               dist = "gumbel_min",
                               links = NULL, call = NULL, ...) {
  parameters <- list(
    d = paste0(
      "Sensitivity: d' under equal variance (g' for gumbel_min). When sdratio ",
      "is estimated, d is d_a, the distance between the target and lure ",
      "distributions in units of their root-mean-square SD"
    )
  )
  # The sd rates match sdt_yn, so d and sdratio shrink the same way whichever
  # SDT model the same subjects are fitted with. Rate 1 for d: its
  # between-subject SD is 0.30 [0.24, 0.37] on meyer_grant_jakob_2025 and ~0.6
  # for sdt_yn on broeder_schuetz_2009_e3, and only rate 1's median (0.69)
  # covers the larger of the two. Rate 2 for sdratio, whose between-subject SD
  # is 0.06 [0.00, 0.12] here and ~0.2 there.
  default_priors <- list(
    d = list(main = "normal(1, 1)", effects = "normal(0, 0.5)", sd = "exponential(1)")
  )
  param_links <- list(d = "identity")
  fixed_pars <- list()
  init_ranges <- list(d = c(0.5, 1.5))

  # Gaussian ranking carries sdratio as an overridable fixed parameter (fixed to
  # 0 = equal variance, sampled when the user adds sdratio ~ ...), as sdt_yn
  # does.
  if (dist == "normal") {
    parameters$sdratio <- paste0(
      "Log SD ratio: log ratio of signal to noise standard deviations ",
      "(exp(sdratio) is the ratio, 0 = equal variance)"
    )
    # Identical to sdt_yn's, whose comment gives the empirical grounds: it
    # reaches past the widest group-level ratio on record (1.71) and keeps its
    # mass off ratios above 2, which recognition does not produce. The
    # quadrature no longer bounds it; it holds 1e-6 over ratios in [0.2, 8]
    # (see .ranking_gh_n and .ranking_gl_rule), 5.4 prior SDs below 0 and 6.9
    # above.
    default_priors$sdratio <- list(
      main = "normal(0, 0.3)", effects = "normal(0, 0.3)", sd = "exponential(2)"
    )
    param_links$sdratio <- "identity"
    fixed_pars$sdratio <- 0
    init_ranges$sdratio <- c(-0.3, 0.3)
  }

  requirements <- glue(
    "Provide pre-aggregated ranking counts in wide format:", "\n\n",
    "  - Rank-count columns (response): one column per rank position, each ",
    "giving the number of trials in which the target received that rank ",
    "(column 1 = most likely target, column m = least)", "\n",
    "  - Set size (m): a constant or a column giving the number of ranked ",
    "items per row; rows with fewer ranks leave the surplus columns at 0", "\n",
    "  No stimulus column needed (all trials include exactly one target)"
  )

  out <- structure(
    list(
      resp_vars = nlist(response),
      other_vars = nlist(m, dist),
      domain = "Perception & Recognition Memory",
      task = "Ranking Task",
      name = "Signal Detection Theory (Ranking)",
      citation = glue(
        "Meyer-Grant, C. G., Kellen, D., Harding, S. M., & Singmann, H. ",
        "(2026). Extreme-value signal detection theory for recognition memory: ",
        "The parametric road not taken. Psychological Review. Advance online ",
        "publication. ",
        "https://doi.org/10.1037/rev0000615"
      ),
      version = "NA",
      requirements = requirements,
      parameters = parameters,
      links = param_links,
      fixed_parameters = fixed_pars,
      default_priors = default_priors,
      init_ranges = init_ranges
    ),
    class = c("bmmodel", "sdt", "sdt_ranking"),
    call = call
  )
  set_links(out, links)
}

# `sdratio` is the log SD ratio itself -- the Stan code reads exp(sdratio) --
# and it is fixed at 0, which means equal variance only while the link is
# identity. Any other link both rescales the fixed value and applies a second
# transformation on top of the exp(), so the model would sample a ratio the
# user never asked for.
#
# `d` is settable only for gumbel_min, whose closed-form kernel is exact far
# beyond any plausible d. The Gaussian quadrature's nodes stop short of the
# last rank's integrand peak once d is large: at equal variance and set sizes
# up to 8 the Stan kernel is off by several nats at d = 15 and returns -Inf
# from d = 30. Under the identity link the normal(1, 1) prior never gets
# there, but under a log link it puts 4.4% of its mass above d = 15 and 0.8%
# above d = 30.
#' @exportS3Method
settable_links.sdt_ranking <- function(model) {
  if (model$other_vars$dist == "gumbel_min") "d" else character(0)
}

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


#' @title Ranking Signal Detection Theory Model
#' @name sdt_ranking
#' @details `r model_docs(.model_sdt_ranking())`
#'
#' Models the rank ordering of `m` items by perceived strength. Only `d`
#' is estimated (no criterion). Supports `dist = "gumbel_min"` (closed-form via
#' lgamma ratios) and `dist = "normal"` (numerical quadrature, accurate to
#' 1e-6 in each rank probability for SD ratios from 0.2 to 8, `d` from -1 to 4
#' and set sizes up to 20).
#'
#' The model uses the native brms multinomial family: each rank position is a
#' multinomial category whose logit is set to `log p(rank)`, so `softmax`
#' recovers the rank distribution exactly. This means `log_lik`,
#' `posterior_predict`, `posterior_epred`, and `pp_check` come from brms as
#' proper joint multinomial draws.
#'
#' @section Sensitivity is on the same scale as [sdt_yn()]:
#' `d` is \eqn{d'} whenever the target and lure distributions share an SD:
#' always for `dist = "gumbel_min"`, where it is the \eqn{g'} of Meyer-Grant et
#' al. (2026), and for `dist = "normal"` unless you give `sdratio` a formula.
#' With `sdratio` estimated, `d` is the balanced index \eqn{d_a} that
#' [sdt_yn()] reports, the separation divided by the root-mean-square of the two
#' SDs; the noise-standardized separation is then
#' `d * sqrt((1 + exp(sdratio)^2) / 2)`.
#'
#' Ranking has no criterion, so \eqn{d_a} is not read off an ROC here. It is
#' adopted to keep `d` on the same scale across the SDT family, and it has a
#' direct meaning for rankings: it fixes the two-item accuracy at
#' \eqn{\Phi(d/\sqrt{2})} whatever `sdratio` is (see below).
#'
#' Ranking is the one SDT design that identifies the variance ratio from a
#' single condition. The rank distribution supplies `m - 1` free probabilities
#' per set size, so with `m >= 3` there is enough information to separate `d`
#' from `sdratio` without the criterion sweep that [sdt_yn()] needs -- the
#' *shape* of the rank distribution, not just its mean, carries the ratio.
#'
#' At `m = 2` the model reduces to 2AFC and the two parameters are no longer
#' separable: the probability of ranking the target first is the area under the
#' yes/no ROC, which for Gaussian noise is \eqn{\Phi(d/\sqrt{2})} whatever
#' `sdratio` is. Keep `sdratio` fixed for two-item designs.
#'
#' For Gaussian ranking (`dist = "normal"`), `sdratio` is fixed to 0 by default.
#' Add `sdratio ~ 1` to the formula for unequal-variance ranking.
#'
#' The set size `m` may be a constant or the name of a data column. Supply a
#' column to fit trials with different set sizes in a single model: the response
#' has `max(m)` columns, and rows with a smaller set size switch off the surplus
#' rank categories (the multinomial then renormalizes over the valid ranks).
#' @param response A character vector of column names with the target rank-count
#'   columns, ordered from rank 1 (most likely target) to rank `m` (least).
#'   With a constant `m`, supply exactly `m` columns; with a varying set size,
#'   supply `max(m)` columns.
#' @param m Either a single integer >= 2 giving the number of ranked items
#'   (constant across all rows), or a single string naming a data column that
#'   gives the number of ranked items per row.
#' @param dist Character. The distribution assumed for the latent evidence,
#'   given here by its cumulative distribution function:
#'   \itemize{
#'     \item "gumbel_min" (default): smallest extreme value,
#'       \eqn{1 - \exp(-\exp(x))} (complementary log-log). Closed form via
#'       gamma-function ratios.
#'     \item "normal": Gaussian, \eqn{\Phi(x)} (supports unequal variance via
#'       `sdratio`)
#'   }
#' @param links A named list of link functions for the parameters. Only the
#'   link of `d` can be changed, and only for `dist = "gumbel_min"`: with
#'   `dist = "normal"` the quadrature loses accuracy at the large `d` a
#'   log link reaches, so `d` keeps the identity link. `sdratio` always
#'   keeps the identity link.
#' @param ... used internally for testing, ignore it
#' @return An object of class `bmmodel`
#' @references
#' Meyer-Grant, C. G., Kellen, D., Harding, S. M., & Singmann, H. (2026).
#'   Extreme-value signal detection theory for recognition memory: The
#'   parametric road not taken. \emph{Psychological Review}. Advance online
#'   publication. \doi{10.1037/rev0000615}
#' @keywords bmmodel
#' @export
#' @examples
#' \dontrun{
#' dat <- data.frame(id = 1:20)
#' dat <- cbind(dat, rsdt_ranking(20, 200, m = 4, d = 1.5))
#'
#' model <- sdt_ranking(
#'   response = c("rank1", "rank2", "rank3", "rank4"),
#'   m = 4
#' )
#'
#' fit <- bmm(
#'   formula = bmf(d ~ 1),
#'   data = dat,
#'   model = model,
#'   cores = 4,
#'   backend = "cmdstanr"
#' )
#'
#' # Mixed set sizes in one model: response has max(m) columns, m is a column
#' model_mixed <- sdt_ranking(
#'   response = c("rank1", "rank2", "rank3", "rank4", "rank5"),
#'   m = "set_size"
#' )
#' }
sdt_ranking <- function(response, m,
                        dist = c("gumbel_min", "normal"),
                        links = NULL, ...) {
  call <- match.call()
  stop_missing_args()
  dist <- match.arg(dist)

  stopif(!((is.numeric(m) && length(m) == 1 && is.finite(m) && m >= 2) ||
           (is.character(m) && length(m) == 1)),
         "m must be a single integer >= 2, or the name of a set-size column in the data")
  if (is.numeric(m)) {
    warnif(m != trunc(m), "m should be an integer value; {m} was truncated to {as.integer(m)}")
    m <- as.integer(m)
  }

  stopif(length(response) < 2,
         "response must name at least 2 rank-count columns")
  stopif(is.numeric(m) && length(response) != m,
         "With a constant m ({m}), response must name exactly m columns ({length(response)} supplied)")

  .model_sdt_ranking(response = response, m = m, dist = dist,
                     links = links, call = call, ...)
}


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

#' @export
check_data.sdt_ranking <- function(model, data, formula) {
  resp_cols <- model$resp_vars$response
  .validate_sdt_count_cols(data, resp_cols)

  max_rank <- .sdt_resolve_set_size(model$other_vars$m, data)
  stopif(anyNA(max_rank),
         "Set-size column '{model$other_vars$m}' must not contain NA")
  n_ranks <- length(resp_cols)
  stopif(any(max_rank > n_ranks, na.rm = TRUE),
         "Set size must not exceed the number of rank columns ({n_ranks})")

  Y <- as.matrix(data[resp_cols])
  # NAs are only allowed as structural blanks beyond a row's set size; an NA
  # within the set size would silently shrink nTrials if zeroed.
  stopif(anyNA(Y[col(Y) <= max_rank]),
         "Response columns must not contain NA counts within the row's set \\
         size (m); only rank columns beyond the set size may be NA")
  Y[is.na(Y)] <- 0
  stopif(any(rowSums(Y) <= 0, na.rm = TRUE),
         "Row sums of response columns must be positive (no empty rows)")

  # Counts in rank columns beyond a row's set size must be structural zeros: a
  # nonzero count there means the set size is mislabelled for that row.
  surplus <- col(Y) > max_rank
  stopif(any(Y[surplus] != 0),
         "Rank columns beyond the row's set size (m) must be 0")

  reserved <- intersect(c("Y", "nTrials", "max_rank"), 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)
  data$max_rank <- as.numeric(max_rank)

  free_sdratio <- "sdratio" %in% names(model$parameters) &&
    !"sdratio" %in% names(model$fixed_parameters)
  warnif(free_sdratio && max(max_rank, na.rm = TRUE) > 20,
         "Estimating `sdratio` with set sizes above 20 goes past the range bmm \\
         has verified its Gaussian quadrature over, so the rank probabilities \\
         may carry errors above 1e-6. Consider fixing `sdratio`, or \\
         `dist = \"gumbel_min\"`, which is closed form")
  warnif(free_sdratio && max(max_rank, na.rm = TRUE) < 3,
         "`sdratio` cannot be identified from two-item rankings: the model \\
         reduces to 2AFC, where the rank probabilities depend only on `d`. \\
         Remove `sdratio` from the formula so it stays fixed at equal variance")

  NextMethod("check_data")
}


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

# sdt_ranking uses brms' native multinomial family: each rank position's logit
# is set to log(p(rank)) so softmax recovers the ranking probabilities. The
# ranking math is computed by the static Stan function sdt_ranking_logmu for
# fitting and by the exported R companion sdt_ranking_logmu() for
# posterior_predict/epred (brms evaluates the non-linear formula in R for
# prediction). log_lik and posterior_predict therefore come from brms — proper
# joint multinomial draws — while both noise distributions stay supported.

#' @export
bmf2bf.sdt_ranking <- function(model, formula) {
  resp_cats <- model$resp_vars$response

  # max_rank (per-row set size) travels as a covariate; sdratio is the
  # parameter for the normal distribution and the literal 0 for gumbel_min;
  # the distribution id selects the kernel in both the Stan function and the
  # R companion.
  sdratio_arg <- if ("sdratio" %in% names(model$parameters)) "sdratio" else "0"
  # d reaches the kernel on the natural scale, so a non-identity link is
  # inverted here, as apply_links() does for m3
  args <- paste("max_rank", deparse(inv_link("d", model$links$d)), sdratio_arg,
                .sdt_dist_id(model$other_vars$dist), sep = ", ")

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


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

#' @export
configure_model.sdt_ranking <- 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)

  # sdt_dist_funs.stan supplies sdt_rms_scale(), which the Gaussian branch needs
  # to turn d_a into the noise-standardized separation
  sc_path <- system.file("stan_chunks", package = "bmm")
  stan_funs <- paste(
    read_lines2(paste0(sc_path, "/sdt_dist_funs.stan")),
    .ranking_fill_quadrature(
      read_lines2(paste0(sc_path, "/sdt_ranking_funs.stan")),
      max_m = max(data$max_rank),
      free_sdratio = "sdratio" %in% names(model$parameters) &&
        !"sdratio" %in% names(model$fixed_parameters)
    ),
    sep = "\n"
  )
  stanvars <- brms::stanvar(scode = stan_funs, block = "functions")

  nlist(formula, data, stanvars)
}


# Substitute the quadrature tokens in sdt_ranking_funs.stan with the rules that
# .ranking_gh_n() and .ranking_gl_rule() select for the data's largest set size;
# .ranking_prob_r() selects per row, so rows with a smaller set size can use a
# different rule (they agree to ~1e-8). 17 significant digits round-trips a
# double exactly.
.ranking_fill_quadrature <- function(scode, max_m, free_sdratio) {
  gh <- .gh_rule(.ranking_gh_n(max_m, free_sdratio))
  gl <- .ranking_gl_rule(max_m)
  brace <- function(x) paste0("({", paste(sprintf("%.17g", x), collapse = ", "), "})")
  tokens <- list(
    "{{N_GH}}" = length(gh$nodes),
    "({{GH_NODES}})" = brace(gh$nodes),
    "({{GH_WEIGHTS}})" = brace(gh$weights),
    "{{N_GL}}" = length(gl$nodes),
    "{{GL_HALF_WIDTH}}" = sprintf("%.17g", gl$half_width),
    "({{GL_NODES}})" = brace(gl$nodes),
    "({{GL_WEIGHTS}})" = brace(gl$weights),
    "({{GL_CDF}})" = brace(gl$cdf)
  )
  for (token in names(tokens)) {
    scode <- gsub(token, tokens[[token]], scode, fixed = TRUE)
  }
  scode
}


#' @title Ranking SDT rank log-probability (multinomial logit)
#' @description R companion to the Stan `sdt_ranking_logmu` function. It returns
#'   `log(p(rank))` for rank position `cat`, which the [sdt_ranking()]
#'   multinomial formula uses as the category logit (so `softmax` recovers the
#'   ranking 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 rank position index.
#' @param max_rank Set size (number of ranked items) for the observation; ranks
#'   above it return the finite `-100` sentinel so the category switches off.
#' @param d Sensitivity (draws-by-observation matrix supplied by brms).
#' @param sdratio Log SD ratio; `0` for the gumbel_min distribution.
#' @param dist Integer noise-distribution id (2 = gumbel_min, 1 = normal).
#' @return `log(p(cat))`, matching the shape of `d`.
#' @keywords internal
#' @export
sdt_ranking_logmu <- function(cat, max_rank, d, sdratio = 0, dist = 2L) {
  dist_name <- .sdt_dist_names[dist]

  shape <- dim(d)
  d <- as.vector(d)
  n <- length(d)
  max_rank <- rep_len(as.vector(max_rank), n)
  sdratio <- rep_len(as.vector(sdratio), n)

  # grouping by the (few) unique set sizes vectorizes the quadrature over all
  # draws-by-observations elements sharing a set size
  out <- rep(-100, n)
  for (m_j in unique(max_rank[max_rank >= cat])) {
    idx <- which(max_rank == m_j)
    probs <- rbind(.ranking_all_probs_r(d[idx], m_j, dist_name,
                                        sdratio[idx]))
    out[idx] <- log(probs[, cat])
  }

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