tests/testthat/test-model_sdt_cdp.R

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

sim_cdp_data <- function(n_subjects = 4, n_trials = 100, n_new = 3, n_old = 3,
                         kcrit = NULL, thresholds = NULL, ...) {
  dat <- expand.grid(id = seq_len(n_subjects), stimulus = c(0L, 1L))
  thresholds <- thresholds %||%
    .cdp_make_thresholds(0, -0.3, n_new, n_old, "parsimonious")
  cbind(dat, rsdt_cdp(nrow(dat), n_trials, dat$stimulus,
                      dfam = 0.8, drec = 1.0, thresholds = thresholds,
                      rcrit = 0.5, kcrit = kcrit, n_new = n_new, ...))
}

test_that("sdt_cdp gives every parameter an sd default prior", {
  rate1 <- c("dfam", "drec")
  for (tt in c("parsimonious", "log_distance")) {
    m <- sdt_cdp(stimulus = "s", n_new = 3, n_old = 3, threshold_type = tt)
    sds <- vapply(m$default_priors, function(p) p$sd %||% NA_character_,
                  character(1))
    expect_false(anyNA(sds), info = tt)
    expect_true(all(sds[rate1] == "exponential(1)"), info = tt)
    expect_true(all(sds[setdiff(names(sds), rate1)] == "exponential(2)"), info = tt)
  }
})

test_that("sdt_cdp only supports dist = 'normal'", {
  expect_error(
    sdt_cdp(stimulus = "s", n_new = 3, n_old = 3, dist = "logistic"),
    "supports only dist"
  )
})

test_that("sdt_cdp requires stimulus, n_new, and n_old", {
  expect_error(sdt_cdp(n_new = 3, n_old = 3))
  expect_error(sdt_cdp(stimulus = "s", n_old = 3))
  expect_error(sdt_cdp(stimulus = "s", n_new = 3))
  expect_error(sdt_cdp(stimulus = "s", n_new = 1, n_old = 1), "at least 3")
})

test_that("sdt_cdp validates n_new, n_old and the response prefix", {
  for (bad in list(NA, Inf, 0, c(2, 3), "3")) {
    expect_error(sdt_cdp(stimulus = "s", n_new = bad, n_old = 3),
                 "n_new must be a single integer >= 1")
    expect_error(sdt_cdp(stimulus = "s", n_new = 3, n_old = bad),
                 "n_old must be a single integer >= 1")
  }
  expect_warning(m <- sdt_cdp(stimulus = "s", n_new = 2.6, n_old = 3),
                 "2.6 was truncated to 2")
  expect_identical(m$other_vars$n_new, 2L)
  expect_error(
    sdt_cdp(response = c("new1", "know2", "remember2"), stimulus = "s",
            n_new = 1, n_old = 1),
    "single prefix"
  )
})

test_that(".cdp_make_thresholds anchors the old/new boundary at criterion", {
  # symmetric: boundary is the middle threshold
  thr <- .cdp_make_thresholds(0.3, -0.2, n_new = 3, n_old = 3, "parsimonious")
  expect_length(thr, 5)
  expect_equal(thr[3], 0.3)             # boundary between bin 3 and 4
  expect_true(!is.unsorted(thr))
  # asymmetric (1 new + 5 old): boundary is thr[1]
  thr_a <- .cdp_make_thresholds(0.5, -0.1, n_new = 1, n_old = 5, "equidistant")
  expect_length(thr_a, 5)
  expect_equal(thr_a[1], 0.5)
  expect_true(!is.unsorted(thr_a))
})

test_that(".cdp_make_thresholds log_distance builds cumulative log-distances", {
  # symmetric 3/3: anchor at thr[3]; delta_i is the log width of the interval
  # between thresholds i and i + 1, as in sdt_rating
  d <- log(c(0.8, 0.6, 0.6, 0.8))   # delta1 ... delta4
  thr <- .cdp_make_thresholds(0.5, 0, 3, 3, "log_distance", deltas = d)
  expect_length(thr, 5)
  expect_equal(thr[3], 0.5)
  expect_equal(thr[4] - thr[3], 0.6)   # exp(delta3)
  expect_equal(thr[5] - thr[4], 0.8)   # exp(delta4)
  expect_equal(thr[3] - thr[2], 0.6)   # exp(delta2)
  expect_equal(thr[2] - thr[1], 0.8)   # exp(delta1)
  expect_true(!is.unsorted(thr))
  # asymmetric 1/5: anchor thr[1], all four distances above
  thr_a <- .cdp_make_thresholds(0.2, 0, 1, 5, "log_distance",
                                deltas = log(c(0.6, 0.7, 0.8, 0.9)))
  expect_equal(thr_a[1], 0.2)
  expect_equal(thr_a[2] - thr_a[1], 0.6)
  expect_true(!is.unsorted(thr_a))
  # asymmetric 4/2: anchor thr[4], three distances below + one above
  thr_b <- .cdp_make_thresholds(0.1, 0, 4, 2, "log_distance",
                                deltas = log(c(0.5, 0.6, 0.7, 0.8)))
  expect_equal(thr_b[4], 0.1)
  expect_equal(thr_b[5] - thr_b[4], 0.8)
  expect_true(!is.unsorted(thr_b))
})

test_that(".cdp_make_thresholds vectorizes over draws", {
  crit <- c(0.1, -0.2, 0.4)
  sp <- c(-0.3, 0, 0.2)
  thr <- .cdp_make_thresholds(crit, sp, 3, 3, "parsimonious")
  expect_equal(dim(thr), c(3, 5))
  for (i in 1:3) {
    expect_equal(thr[i, ],
                 .cdp_make_thresholds(crit[i], sp[i], 3, 3, "parsimonious"))
  }
  deltas <- rbind(log(c(0.8, 0.6, 0.6, 0.8)), log(c(0.5, 0.7, 0.9, 1.1)))
  thr_ld <- .cdp_make_thresholds(c(0.5, -0.1), 0, 3, 3, "log_distance",
                                 deltas = deltas)
  expect_equal(dim(thr_ld), c(2, 5))
  for (i in 1:2) {
    expect_equal(thr_ld[i, ],
                 .cdp_make_thresholds(deltas = deltas[i, ],
                                      c(0.5, -0.1)[i], 0, 3, 3,
                                      "log_distance"))
  }
})

test_that("sdt_cdp log_distance names its deltas in interval order and has no spacing", {
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3,
               threshold_type = "log_distance")
  deltas <- paste0("delta", 1:4)
  expect_equal(.sdt_threshold_delta_names(m), deltas)
  expect_false("spacing" %in% names(m$parameters))
  expect_true(all(deltas %in% names(m$init_ranges)))
})

test_that("sdt_cdp log_distance integrates with the pipeline via mock backend", {
  skip_on_cran() # over 0.4 s in CRAN mode (ADR 0004, R10)
  thr <- .cdp_make_thresholds(0.3, 0, 3, 3, "log_distance",
                              deltas = log(c(0.8, 0.6, 0.6, 0.8)))
  dat <- sim_cdp_data(thresholds = thr)
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3,
               threshold_type = "log_distance")
  expect_silent(
    bmm(bmf(dfam ~ 1, drec ~ 1, criterion ~ 1, rcrit ~ 1,
            delta1 ~ 1, delta2 ~ 1, delta3 ~ 1, delta4 ~ 1),
        dat, m, backend = "mock", mock_fit = 1, rename = FALSE)
  )
})

test_that(".sdt_cdp_logmu_stan generates a log_distance wrapper with delta args", {
  m <- sdt_cdp(stimulus = "s", n_new = 3, n_old = 3,
               threshold_type = "log_distance")
  sc <- .sdt_cdp_logmu_stan(m)
  expect_true(grepl("real sdt_cdp_logmu", sc))
  expect_true(grepl("real delta1.*real delta2.*real delta3.*real delta4", sc))
  expect_true(grepl("array\\[4\\] real deltas", sc))
  expect_true(grepl("cdp_make_thresholds\\(criterion, spacing, deltas", sc))
})

test_that(".sdt_cdp_category_probs sums to 1 across variants", {
  thr <- .cdp_make_thresholds(0, -0.3, 3, 3, "parsimonious")
  # R/K, target and lure
  pt <- .sdt_cdp_category_probs(thr, 0.8, 1.0, 0, 0.5, NULL, 1, 3, 3, "normal", 0)
  pl <- .sdt_cdp_category_probs(thr, 0.8, 1.0, 0, 0.5, NULL, 0, 3, 3, "normal", 0)
  expect_length(pt, 9)
  expect_equal(sum(pt), 1, tolerance = 1e-10)
  expect_equal(sum(pl), 1, tolerance = 1e-10)
  # unequal variance + correlation
  puv <- .sdt_cdp_category_probs(thr, 0.8, 1.0, log(1.5), 0.5, NULL, 1, 3, 3,
                                 "normal", 0.4)
  expect_equal(sum(puv), 1, tolerance = 1e-10)
  # 3-way R/K/G
  p3 <- .sdt_cdp_category_probs(thr, 0.8, 1.0, 0, 0.6, 0.1, 1, 3, 3, "normal", 0)
  expect_length(p3, 12)
  expect_equal(sum(p3), 1, tolerance = 1e-10)
})

test_that(".sdt_cdp_category_probs guess + know recover the 2-way know total", {
  thr <- .cdp_make_thresholds(0, -0.3, 3, 3, "parsimonious")
  p2 <- .sdt_cdp_category_probs(thr, 0.8, 1.0, 0, 0.6, NULL, 1, 3, 3, "normal", 0)
  p3 <- .sdt_cdp_category_probs(thr, 0.8, 1.0, 0, 0.6, 0.1, 1, 3, 3, "normal", 0)
  # know(2-way) == guess + know (3-way) per old bin, after renormalising the
  # 2-way vector is identical because the bin totals are unchanged
  know2 <- p2[4:6]
  guess3 <- p3[4:6]
  know3 <- p3[7:9]
  expect_equal(know2, guess3 + know3, tolerance = 1e-8)
})

test_that(".sdt_cdp_category_probs vectorizes over observations", {
  n <- 10
  thr <- .cdp_make_thresholds(rnorm(n, 0, 0.2), rnorm(n, -0.3, 0.1), 3, 3,
                              "parsimonious")
  df <- rnorm(n, 0.8, 0.2)
  dr <- rnorm(n, 1.0, 0.2)
  sg <- rnorm(n, 0.1, 0.1)
  rh <- rnorm(n, 0.2, 0.2)
  rc <- rnorm(n, 0.5, 0.1)
  kc <- rnorm(n, -0.3, 0.1)
  st <- rep_len(c(0L, 1L), n)
  pv <- .sdt_cdp_category_probs(thr, df, dr, sg, rc, kc, st, 3, 3, "normal", rh)
  expect_equal(dim(pv), c(n, 12))
  for (j in seq_len(n)) {
    expect_equal(pv[j, ],
                 .sdt_cdp_category_probs(thr[j, ], df[j], dr[j], sg[j], rc[j],
                                         kc[j], st[j], 3, 3, "normal", rh[j]))
  }
})

test_that("sdt_cdp_logmu matches .sdt_cdp_category_probs and preserves shape", {
  thr_name <- "parsimonious"
  args <- list(n_new = 3, n_old = 3, thresh = 1L, has_guess = 0L,
               dfam = 0.8, drec = 1.0, criterion = 0.1, spacing = -0.3,
               rcrit = 0.5, sigmar = 0, rho = 0.2, kcrit = -100, stimulus = 1)
  thr <- .cdp_make_thresholds(args$criterion, args$spacing, 3, 3, thr_name)
  probs <- .sdt_cdp_category_probs(thr, args$dfam, args$drec, args$sigmar,
                                   args$rcrit, NULL, 1, 3, 3, "normal", args$rho)
  for (cat in seq_len(9)) {
    val <- sdt_cdp_logmu(cat, 3, 3, 1L, 0L, args$dfam, args$drec,
                         args$criterion, args$spacing, args$rcrit, args$sigmar,
                         args$rho, args$kcrit, args$stimulus)
    expect_equal(val, log(probs[cat]), tolerance = 1e-10)
  }
  # draws-by-observation shape is preserved
  mat <- matrix(0.8, nrow = 4, ncol = 2)
  out <- sdt_cdp_logmu(1, 3, 3, 1L, 0L, mat, mat, mat, matrix(-0.3, 4, 2),
                       mat, matrix(0, 4, 2), matrix(0, 4, 2),
                       matrix(-100, 4, 2), matrix(1, 4, 2))
  expect_equal(dim(out), c(4, 2))
})

test_that("sdt_cdp_logmu vectorized path matches per-element evaluation", {
  nd <- 6
  nobs <- 4
  withr::local_seed(11)
  mk <- function(center, spread) matrix(rnorm(nd * nobs, center, spread), nd, nobs)
  a <- list(dfam = mk(0.8, 0.2), drec = mk(1.0, 0.2),
            criterion = mk(0.1, 0.1), spacing = mk(-0.3, 0.1),
            rcrit = mk(0.5, 0.1), sigmar = mk(0.1, 0.1), rho = mk(0.2, 0.1),
            kcrit = mk(-0.3, 0.1),
            stimulus = matrix(rep(c(0, 0, 1, 1), each = nd), nd, nobs))
  for (cat in c(1L, 5L, 9L, 12L)) {
    vec <- sdt_cdp_logmu(cat, 3, 3, 1L, 1L, a$dfam, a$drec, a$criterion,
                         a$spacing, a$rcrit, a$sigmar, a$rho, a$kcrit,
                         a$stimulus)
    ref <- matrix(0, nd, nobs)
    for (i in seq_len(nd)) {
      for (j in seq_len(nobs)) {
        ref[i, j] <- sdt_cdp_logmu(cat, 3, 3, 1L, 1L, a$dfam[i, j],
                                   a$drec[i, j], a$criterion[i, j],
                                   a$spacing[i, j], a$rcrit[i, j],
                                   a$sigmar[i, j], a$rho[i, j], a$kcrit[i, j],
                                   a$stimulus[i, j])
      }
    }
    expect_equal(vec, ref, tolerance = 1e-12)
  }
})

test_that("rsdt_cdp adds guess columns and honours per-observation parameters", {
  withr::local_seed(2203)
  dat <- expand.grid(id = 1:5, stimulus = c(0L, 1L))
  thr <- .cdp_make_thresholds(0, -0.3, 3, 3, "parsimonious")
  # guess columns appear with a finite kcrit
  cnts_g <- rsdt_cdp(nrow(dat), 100, dat$stimulus, 0.8, 1.0, thresholds = thr,
                     rcrit = 0.5, kcrit = -0.3, n_new = 3)
  expect_true(all(c("guess4", "guess5", "guess6") %in% colnames(cnts_g)))
  # per-observation parameters are honoured
  df_v <- seq(0.2, 1.4, length.out = nrow(dat))
  cnts_v <- rsdt_cdp(nrow(dat), 10000, dat$stimulus, df_v, 1.0, thresholds = thr,
                     rcrit = 0.5, n_new = 3)
  expect_equal(dim(cnts_v), c(nrow(dat), 9))
})

test_that("dsdt_cdp matches dmultinom against the category probabilities", {
  thr <- .cdp_make_thresholds(0, -0.3, 3, 3, "parsimonious")
  cnt <- c(40, 20, 10, 15, 25, 35, 5, 20, 80)
  dl <- dsdt_cdp(counts = cnt, stimulus = 1, dfam = 0.8, drec = 1.0,
                 thresholds = thr, rcrit = 0.5, n_new = 3, log = TRUE)
  probs <- .sdt_cdp_category_probs(thr, 0.8, 1.0, 0, 0.5, NULL, 1, 3, 3,
                                   "normal", 0)
  expect_equal(dl, dmultinom(cnt, prob = probs, log = TRUE), tolerance = 1e-10)
})

test_that("aggregate_sdt_cdp_data pivots long data to the wide count columns", {
  long <- expand.grid(id = 1:4, stimulus = c(0L, 1L),
                      judgment = c("new", "know", "remember"),
                      confidence = 1:6)
  long <- long[(long$judgment == "new") == (long$confidence <= 3), ]
  long$count <- sample(5:20, nrow(long), replace = TRUE)
  wide <- aggregate_sdt_cdp_data(long, "judgment", "confidence", "count")
  expect_equal(nrow(wide), 8)                # 4 subjects x 2 stimulus classes
  expect_true(all(c("new1", "know4", "remember6") %in% colnames(wide)))
  expect_equal(sum(wide[grep("new|know|remember", colnames(wide))]),
               sum(long$count))
  # a response prefix is applied to the count columns
  wide_p <- aggregate_sdt_cdp_data(long, "judgment", "confidence", "count",
                                   response = "cdp")
  expect_true(all(c("cdpnew1", "cdpknow4") %in% colnames(wide_p)))
})

test_that("aggregate_sdt_cdp_data aggregates trial-level data without counts", {
  long <- expand.grid(id = 1:3, stimulus = c(0L, 1L),
                      judgment = c("new", "know", "remember"),
                      confidence = 1:6)
  long <- long[(long$judgment == "new") == (long$confidence <= 3), ]
  long$count <- sample(2:8, nrow(long), replace = TRUE)
  trials <- long[rep(seq_len(nrow(long)), long$count),
                 c("id", "stimulus", "judgment", "confidence")]
  wide <- aggregate_sdt_cdp_data(trials, "judgment", "confidence")
  expect_equal(nrow(wide), 6)
  expect_equal(sum(wide[grep("new|know|remember", colnames(wide))]),
               sum(long$count))
})

test_that("aggregate_sdt_cdp_data emits guess columns and validates input", {
  base <- data.frame(
    id = 1, stimulus = 1,
    judgment = c("new", "know", "remember", "guess"),
    confidence = c(1, 2, 2, 2), count = c(10, 5, 7, 3)
  )
  wide <- aggregate_sdt_cdp_data(base, "judgment", "confidence", "count")
  expect_true("guess2" %in% colnames(wide))
  # bad judgment label
  bad_j <- base
  bad_j$judgment[1] <- "maybe"
  expect_error(aggregate_sdt_cdp_data(bad_j, "judgment", "confidence", "count"),
               "not allowed")
  # missing a required judgment type
  miss <- base[base$judgment != "remember", ]
  expect_error(aggregate_sdt_cdp_data(miss, "judgment", "confidence", "count"),
               "must contain")
  # confidence partition violated: an old judgment on a new level
  bad_c <- base
  bad_c$confidence[2] <- 1
  expect_error(aggregate_sdt_cdp_data(bad_c, "judgment", "confidence", "count"),
               "must occupy levels")
})

test_that("check_data.sdt_cdp binds the count columns into a multinomial matrix", {
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  dat <- sim_cdp_data(n_trials = 200)
  cd <- check_data(m, dat, bmf(dfam ~ 1))
  expect_false(attr(cd, "has_guess"))
  expect_equal(ncol(cd$Y), 9)
  expect_equal(nrow(cd), 8)                  # 4 subjects x 2 stimulus classes
  expect_true(all(cd$nTrials == 200))
})

test_that("check_data.sdt_cdp handles asymmetric scales and guess columns", {
  m_a <- sdt_cdp(stimulus = "stimulus", n_new = 1, n_old = 5)
  thr_a <- .cdp_make_thresholds(0.2, -0.2, 1, 5, "parsimonious")
  dat_a <- sim_cdp_data(n_subjects = 3, n_trials = 300, n_new = 1, n_old = 5,
                        thresholds = thr_a, sigmar = log(1.5))
  cd_a <- check_data(m_a, dat_a, bmf(dfam ~ 1))
  expect_equal(ncol(cd_a$Y), 11)

  m_g <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  dat_g <- sim_cdp_data(n_subjects = 3, n_trials = 200, kcrit = 0.0)
  cd_g <- check_data(m_g, dat_g, bmf(dfam ~ 1))
  expect_true(attr(cd_g, "has_guess"))
  expect_equal(ncol(cd_g$Y), 12)
})

test_that("check_data.sdt_cdp checks the declared scale and the guess columns", {
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  dat <- sim_cdp_data()
  # a declared scale that mismatches the columns is a missing-column error
  m_big <- sdt_cdp(stimulus = "stimulus", n_new = 2, n_old = 4)
  expect_error(check_data(m_big, dat, bmf(dfam ~ 1)), "missing in the data")
  # partial guess columns
  dat_pg <- dat
  dat_pg$guess4 <- 1L
  expect_error(check_data(m, dat_pg, bmf(dfam ~ 1)), "all of them or none")
})

test_that("sdt_cdp accepts a response prefix for the count columns", {
  skip_on_cran() # over 0.4 s in CRAN mode (ADR 0004, R10)
  dat <- sim_cdp_data()
  resp_cols <- grep("new|know|remember", colnames(dat))
  colnames(dat)[resp_cols] <- paste0("cdp", colnames(dat)[resp_cols])
  m <- sdt_cdp(response = "cdp", stimulus = "stimulus", n_new = 3, n_old = 3)
  cd <- check_data(m, dat, bmf(dfam ~ 1))
  expect_equal(ncol(cd$Y), 9)
  expect_silent(
    bmm(bmf(dfam ~ 1, drec ~ 1, criterion ~ 1, spacing ~ 1, rcrit ~ 1),
        dat, m, backend = "mock", mock_fit = 1, rename = FALSE)
  )
})

test_that("bmf2bf builds a multinomial non-linear formula with one mu per category", {
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  dat <- sim_cdp_data(n_subjects = 3)
  cd <- check_data(m, dat, bmf(dfam ~ 1))
  m$other_vars$has_guess <- attr(cd, "has_guess")
  bf <- bmf2bf(m, bmf(dfam ~ 1))
  expect_s3_class(bf, "brmsformula")
  # 9 categories => 1 main + 8 mu non-linear formulas referencing sdt_cdp_logmu
  pforms <- vapply(bf$pforms, function(f) paste(deparse(f), collapse = " "),
                   character(1))
  expect_true(sum(grepl("sdt_cdp_logmu", pforms)) >= 8)
  # each category calls the kernel with its own index
  expect_match(deparse(bf$formula, width.cutoff = 500), "sdt_cdp_logmu(1, 3,",
               fixed = TRUE)
  for (k in 2:9) {
    pform <- deparse(bf$pforms[[paste0("mucdp", k)]], width.cutoff = 500)
    expect_match(pform, paste0("sdt_cdp_logmu(", k, ", 3,"), fixed = TRUE)
  }
})

test_that("a link on dfam, drec, criterion or rcrit is inverted inside the cdp formula", {
  # the multinomial family has no link of its own for these parameters, so a
  # link the user sets reaches the kernel only through the formula
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3,
               links = list(dfam = "log", drec = "softplus",
                            criterion = "probit", rcrit = "logit"))
  m$other_vars$has_guess <- FALSE
  bf <- bmf2bf(m, bmf(dfam ~ 1))
  for (f in c(list(bf$formula), bf$pforms[paste0("mucdp", 2:9)])) {
    txt <- paste(deparse(f, width.cutoff = 500), collapse = "")
    expect_match(txt, "(\\d+, ){4}exp\\(dfam\\), log1p_exp\\(drec\\), Phi\\(criterion\\), spacing, inv_logit\\(rcrit\\), sigmar,")
  }
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  m$other_vars$has_guess <- FALSE
  expect_match(deparse(bmf2bf(m, bmf(dfam ~ 1))$formula, width.cutoff = 500),
               "(\\d+, ){4}dfam, drec, criterion, spacing, rcrit, sigmar,")
})

test_that("sdt_cdp produces multinomial stancode with the CDP functions", {
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  dat <- sim_cdp_data()
  sc <- stancode(bmf(dfam ~ 1, drec ~ 1, criterion ~ 1, spacing ~ 1,
                     rcrit ~ 1), dat, m)
  expect_true(grepl("sdt_cdp_logmu", sc))
  expect_true(grepl("cdp_Phi2", sc))
  expect_true(grepl("owens_t", sc))
  expect_true(grepl("multinomial", sc))
})

test_that("freeing rho and sigmar adds them to the stancode", {
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  dat <- sim_cdp_data()
  sc <- stancode(bmf(dfam ~ 1, drec ~ 1, criterion ~ 1, spacing ~ 1,
                     rcrit ~ 1, sigmar ~ 1, rho ~ 1), dat, m)
  expect_true(grepl("rho", sc))
  expect_true(grepl("sigmar", sc))
})

############################################################################# !
# KERNEL ACCURACY: independent reference and Stan parity                  ####
############################################################################# !

# A grid of category cells out to the tails for every threshold type and the
# three scales 3/3, 1/5 and 4/2 (K = 6 throughout, so log_distance always has
# four deltas): strength bins far above and below both means, a Remember
# criterion inside and far beyond the recollection distribution, strong
# positive and negative F-R correlation, R/K and R/K/G. The thresholds use a
# non-zero spacing and unequal deltas so the builders' anchors matter, and
# sigmar is crossed with criterion rather than cycled alongside it.
cdp_kernel_configs <- function() {
  cf <- expand.grid(n_new = c(3L, 1L, 4L), thresh = 1:3)
  cf$n_old <- 6L - cf$n_new
  cf
}

cdp_kernel_grid <- function() {
  pars <- expand.grid(criterion = c(-5, 0, 5), stimulus = c(0, 1),
                      rho = c(-1.5, 1.5), rcrit = c(0.9, 5), has_guess = 0:1)
  idx <- as.integer(factor(pars$criterion)) + as.integer(factor(pars$rho)) +
    as.integer(factor(pars$rcrit))
  pars$sigmar <- c(0, 0.7, -0.5)[idx %% 3 + 1]
  pars$kcrit <- ifelse(pars$has_guess == 1, 0.5, -100)
  pars$spacing <- -0.3
  cf <- cdp_kernel_configs()
  cells <- do.call(rbind, lapply(seq_len(nrow(cf)), function(i) {
    g <- pars
    g$n_new <- cf$n_new[i]
    g$n_old <- cf$n_old[i]
    g$thresh <- cf$thresh[i]
    n_cat <- ifelse(g$has_guess == 1, 3L, 2L) * g$n_old + g$n_new
    g <- g[rep(seq_len(nrow(g)), n_cat), ]
    g$cat <- unlist(lapply(n_cat, seq_len))
    g
  }))
  cells$dfam <- 0.8
  cells$drec <- 1.0
  rownames(cells) <- NULL
  cells
}

cdp_kernel_deltas <- log(c(0.4, 1.6, 0.4, 1.6))

# the documented threshold rules, written out independently of
# .cdp_make_thresholds() and cdp_make_thresholds()
cdp_kernel_thresholds <- function(g) {
  K <- g$n_new + g$n_old
  k <- seq_len(K - 1L)
  if (g$thresh == 1L) {
    g$criterion + exp(g$spacing) *
      (log(k / (K - k)) - log(g$n_new / (K - g$n_new)))
  } else if (g$thresh == 2L) {
    g$criterion + (k - g$n_new) * exp(g$spacing)
  } else {
    w <- exp(cdp_kernel_deltas)
    thr <- numeric(K - 1L)
    thr[g$n_new] <- g$criterion
    for (j in seq_len(K - 1L - g$n_new)) thr[g$n_new + j] <- thr[g$n_new + j - 1L] + w[g$n_new + j - 1L]
    for (j in rev(seq_len(g$n_new - 1L))) thr[j] <- thr[j + 1L] - w[j]
    thr
  }
}

# the exported R companion of the Stan wrapper, with the model's own
# threshold builder and delta order
cdp_kernel_r <- function(cells) {
  vapply(seq_len(nrow(cells)), function(i) {
    g <- cells[i, ]
    args <- c(list(g$cat, g$n_new, g$n_old, g$thresh, g$has_guess, g$dfam,
                   g$drec, g$criterion, g$spacing, g$rcrit, g$sigmar, g$rho,
                   g$kcrit, g$stimulus),
              if (g$thresh == 3L) as.list(cdp_kernel_deltas))
    do.call(sdt_cdp_logmu, args)
  }, numeric(1))
}

# log P(S in the cell's bin, R and F in its region), integrated over R with
# stats::integrate() and F | R normal: shares no code with Owen's T, TVPACK or
# the kernels' Gauss-Legendre band quadrature
cdp_kernel_reference <- function(g) {
  thr <- cdp_kernel_thresholds(g)
  K <- g$n_new + g$n_old
  if (g$cat <= g$n_new) {
    type <- 1L
    k <- g$cat
  } else {
    r <- g$cat - g$n_new
    block <- (r - 1L) %/% g$n_old
    type <- if (g$has_guess == 1) block + 2L else block + 3L
    k <- g$n_new + r - block * g$n_old
  }
  c_lo <- c(-Inf, thr)[k]
  c_hi <- c(thr, Inf)[k]
  r_lo <- if (type == 4L) g$rcrit else -Inf
  r_hi <- if (type %in% 2:3) g$rcrit else Inf
  f_lo <- if (type == 3L && g$has_guess == 1) g$kcrit else -Inf
  f_hi <- if (type == 2L) g$kcrit else Inf
  old <- g$stimulus == 1
  mu_F <- if (old) g$dfam else 0
  mu_R <- if (old) g$drec else 0
  sd_R <- if (old) exp(g$sigmar) else 1
  corr <- tanh(g$rho)
  log_f <- function(r) {
    m <- mu_F + corr * (r - mu_R) / sd_R
    s <- sqrt(1 - corr^2)
    lo <- (pmax(c_lo - r, f_lo) - m) / s
    hi <- (pmin(c_hi - r, f_hi) - m) / s
    band <- ifelse(lo > 0,
                   stats::pnorm(-lo) - stats::pnorm(-hi),
                   stats::pnorm(hi) - stats::pnorm(lo))
    stats::dnorm(r, mu_R, sd_R, log = TRUE) + log(pmax(band, 0))
  }
  grid <- seq(max(r_lo, mu_R - 40 * sd_R), min(r_hi, mu_R + 40 * sd_R),
              length.out = 4001)
  lg <- log_f(grid)
  if (!any(is.finite(lg))) return(-Inf)
  top <- max(lg)
  keep <- range(which(lg > top - 60))
  a <- grid[max(1, keep[1] - 1)]
  b <- grid[min(4001, keep[2] + 1)]
  # the band's bounds switch expression where a strength threshold meets the
  # kcrit bound; those points are piece boundaries so QUADPACK never straddles
  # a kink
  kinks <- c(c_lo - f_lo, c_hi - f_hi, c_lo - f_hi, c_hi - f_lo)
  br <- sort(unique(c(seq(a, b, length.out = 21), kinks[kinks > a & kinks < b])))
  top + log(sum(vapply(seq_len(length(br) - 1L), function(j) {
    stats::integrate(function(r) exp(log_f(r) - top), br[j], br[j + 1],
                     rel.tol = 1e-10, abs.tol = 0, subdivisions = 2000L,
                     stop.on.error = FALSE)$value
  }, numeric(1))))
}

test_that("the R cdp kernel matches an independent reference in both tails", {
  skip_on_cran() # over 0.4 s in CRAN mode (ADR 0004, R10)
  cells <- cdp_kernel_grid()
  ref <- vapply(seq_len(nrow(cells)), function(i) {
    cdp_kernel_reference(cells[i, ])
  }, numeric(1))
  r <- cdp_kernel_r(cells)
  # down to log p = -30 the kernel resolves the category mass; the old Guess
  # quadrature was off by up to 43 nats at log p near -3 on this grid
  ok <- ref > -30
  expect_gt(sum(ok & ref < -15), 100)
  for (cf in split(seq_len(nrow(cells)), cells[c("n_new", "thresh")])) {
    expect_lt(max(abs(r - ref)[cf][ok[cf]]), 1e-5)
  }
})

# the chunk plus the generated wrappers, as bmm emits them: sdt_cdp_logmu for
# parsimonious/equidistant and, renamed, the four-delta log_distance one
cdp_kernel_stan <- function(cells) {
  m_ld <- sdt_cdp(stimulus = "s", n_new = 3, n_old = 3,
                  threshold_type = "log_distance")
  wrapper_ld <- sub("real sdt_cdp_logmu(", "real sdt_cdp_logmu_ld(",
                    .sdt_cdp_logmu_stan(m_ld), fixed = TRUE)
  program <- paste0(
    "functions {\n",
    read_lines2(system.file("stan_chunks", "sdt_cdp_funs.stan", package = "bmm")),
    "\n", .sdt_cdp_logmu_stan(sdt_cdp(stimulus = "s", n_new = 3, n_old = 3)),
    wrapper_ld, "}\n",
    "data { int N; array[N] int cat; array[N] int n_new; array[N] int n_old;\n",
    "  array[N] int thresh; array[N] int has_guess; vector[N] criterion;\n",
    "  vector[N] spacing; vector[N] dfam; vector[N] drec; vector[N] sigmar;\n",
    "  vector[N] rho; vector[N] rcrit; vector[N] kcrit; vector[N] stimulus;\n",
    "  vector[4] deltas; }\n",
    "generated quantities { vector[N] lp;\n",
    "  for (i in 1:N) {\n",
    "    if (thresh[i] == 3) {\n",
    "      lp[i] = sdt_cdp_logmu_ld(cat[i], n_new[i], n_old[i], 3, has_guess[i],\n",
    "        dfam[i], drec[i], criterion[i], 0, rcrit[i], sigmar[i], rho[i],\n",
    "        kcrit[i], stimulus[i], deltas[1], deltas[2], deltas[3], deltas[4]);\n",
    "    } else {\n",
    "      lp[i] = sdt_cdp_logmu(cat[i], n_new[i], n_old[i], thresh[i], has_guess[i],\n",
    "        dfam[i], drec[i], criterion[i], spacing[i], rcrit[i], sigmar[i], rho[i],\n",
    "        kcrit[i], stimulus[i]);\n",
    "    }\n",
    "  }\n}\n"
  )
  data <- c(list(N = nrow(cells), deltas = cdp_kernel_deltas),
            as.list(cells[c("cat", "n_new", "n_old", "thresh", "has_guess",
                            "criterion", "spacing", "dfam", "drec", "sigmar",
                            "rho", "rcrit", "kcrit", "stimulus")]))
  fit <- cmdstanr::cmdstan_model(cmdstanr::write_stan_file(program))$sample(
    data = data, fixed_param = TRUE, chains = 1, iter_sampling = 1,
    iter_warmup = 0, refresh = 0, show_messages = FALSE, sig_figs = 18
  )
  csv <- utils::read.csv(fit$output_files()[1], comment.char = "#",
                         check.names = FALSE)
  as.numeric(csv[1, paste0("lp.", seq_len(nrow(cells)))])
}

test_that("the Stan cdp kernel matches the R kernel in both tails", {
  skip_unless_stan_tests()
  skip_on_cran()

  cells <- cdp_kernel_grid()
  stan <- cdp_kernel_stan(cells)
  r <- cdp_kernel_r(cells)
  # Stan's Owen's T and R's bivariate normal share no code; on the #372
  # review's 109,080-cell grid they agree to 2.6e-4 in log p for true log p
  # in (-30, -15] and to 1e-8 above -15, where the Owen's T constant term has
  # not started to cancel
  body <- r > -15
  tail <- r > -30 & !body
  expect_gt(sum(tail), 100)
  for (cf in split(seq_len(nrow(cells)), cells[c("n_new", "thresh")])) {
    expect_lt(max(abs(stan - r)[cf][body[cf]]), 1e-8)
    expect_lt(max(abs(stan - r)[cf][tail[cf]]), 1e-3)
  }
})

test_that("the Stan cdp likelihood has a finite gradient with an unbounded top bin", {
  # the parity test compares values only; an infinite bin bound that reaches
  # Phi() as a function of the parameters gives a NaN gradient (zero density
  # times an infinite adjoint) while every value stays right
  skip_unless_stan_tests()
  skip_on_cran()

  program <- paste0(
    "functions {\n",
    read_lines2(system.file("stan_chunks", "sdt_cdp_funs.stan", package = "bmm")),
    "\n}\n",
    "data { array[2, 11] int y; }\n",
    "parameters { real dfam; real drec; real criterion; real spacing;\n",
    "  real rcrit; real sigmar; }\n",
    "model {\n",
    "  array[0] real deltas;\n",
    "  vector[5] thr = cdp_make_thresholds(criterion, spacing, deltas, 1, 5, 1);\n",
    "  for (s in 1:2) {\n",
    "    vector[11] lp;\n",
    "    for (c in 1:11) lp[c] = log(cdp_category_prob(c, thr, dfam, drec,\n",
    "      sigmar, 0, rcrit, -100, s - 1, 1, 5, 0));\n",
    "    y[s] ~ multinomial_logit(lp);\n",
    "  }\n}\n"
  )
  # Rotello et al. (2005) neutral condition, as in the dual-process article
  y <- rbind(c(745, 199, 124, 84, 45, 13, 46, 29, 47, 51, 55),
             c(290, 122, 106, 101, 92, 59, 26, 38, 61, 109, 433))
  model <- cmdstanr::cmdstan_model(cmdstanr::write_stan_file(program))
  inits <- list(dfam = 0.8, drec = 0.8, criterion = 0.1, spacing = -0.4,
                rcrit = 0.9, sigmar = 0.4)
  diag <- suppressMessages(model$diagnose(data = list(y = y), init = list(inits)))
  grad <- diag$gradients()
  expect_equal(nrow(grad), length(inits))
  expect_true(all(is.finite(grad$model)))
  expect_lt(max(abs(grad$error)), 1e-5)
})

############################################################################# !
# KNOW/GUESS SWITCH: guess columns and kcrit must agree                   ####
############################################################################# !

test_that("guess columns without kcrit in the formula are refused, kcrit without guess columns warns", {
  skip_on_cran() # over 0.4 s in CRAN mode (ADR 0004, R10)
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  dat_rk <- sim_cdp_data()
  dat_rkg <- sim_cdp_data(kcrit = 0.2)
  f_rk <- bmf(dfam ~ 1, drec ~ 1, criterion ~ 1, spacing ~ 1, rcrit ~ 1)
  f_rkg <- bmf(dfam ~ 1, drec ~ 1, criterion ~ 1, spacing ~ 1, rcrit ~ 1,
               kcrit ~ 1)
  mock <- function(f, d) {
    bmm(f, d, m, backend = "mock", mock_fit = 1, rename = FALSE)
  }
  # the documented R/K formula on R/K/G data: with kcrit fixed at -100 every
  # Guess category sits at the floor and the fit converges to wrong values
  expect_error(mock(f_rk, dat_rkg), "kcrit must be estimated")
  # the mirror: nothing in R/K data identifies kcrit
  expect_warning(mock(f_rkg, dat_rk), "posterior will equal its prior")
  # the two intended pairings run
  expect_silent(mock(f_rkg, dat_rkg))
  expect_silent(mock(f_rk, dat_rk))
  # a constant kcrit in the formula is a fixed kcrit, not a free one
  expect_error(
    mock(bmf(dfam ~ 1, drec ~ 1, criterion ~ 1, spacing ~ 1, rcrit ~ 1,
             kcrit = 0.2), dat_rkg),
    "kcrit must be estimated"
  )
})

test_that("the Know/Guess switch is checked on update()'s path too", {
  # update.bmmfit() runs check_data() and then check_formula() on the new data
  # with the fit's formula, exactly this call sequence; the guard therefore
  # fires for update(rk_fit, newdata = <R/K/G data>) as it does for bmm()
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  f_rk <- bmf(dfam ~ 1, drec ~ 1, criterion ~ 1, spacing ~ 1, rcrit ~ 1)
  m <- check_model(m, sim_cdp_data(), f_rk)
  cd <- check_data(m, sim_cdp_data(kcrit = 0.2), f_rk)
  expect_true(attr(cd, "has_guess"))
  expect_error(check_formula(m, cd, f_rk), "kcrit must be estimated")
})

test_that("check_data.sdt_cdp refuses count columns beyond the declared scale", {
  # 1 new + 3 old levels in the data, 1 + 2 declared: the know4/remember4
  # columns used to be dropped without a word, and their trials with them
  thr <- .cdp_make_thresholds(0.2, -0.2, 1, 3, "parsimonious")
  dat <- sim_cdp_data(n_subjects = 3, n_trials = 300, n_new = 1, n_old = 3,
                      thresholds = thr)
  m <- sdt_cdp(stimulus = "stimulus", n_new = 1, n_old = 2)
  expect_error(check_data(m, dat, bmf(dfam ~ 1)),
               "'know4', 'remember4' look like response counts beyond")
  # the same with a column prefix
  resp <- grep("new|know|remember", colnames(dat))
  colnames(dat)[resp] <- paste0("cdp", colnames(dat)[resp])
  m_p <- sdt_cdp(response = "cdp", stimulus = "stimulus", n_new = 1, n_old = 2)
  expect_error(check_data(m_p, dat, bmf(dfam ~ 1)), "'cdpknow4', 'cdpremember4'")
  # a column that only looks similar is left alone
  dat$new_trials <- 300L
  dat$knowledge2 <- 1L
  m_ok <- sdt_cdp(response = "cdp", stimulus = "stimulus", n_new = 1, n_old = 3)
  expect_equal(ncol(check_data(m_ok, dat, bmf(dfam ~ 1))$Y), 7)
})

test_that("check_data and aggregate_sdt_cdp_data refuse malformed input", {
  # an empty row
  m <- sdt_cdp(stimulus = "stimulus", n_new = 3, n_old = 3)
  dat <- sim_cdp_data()
  dat[1, grep("new|know|remember", colnames(dat))] <- 0L
  expect_error(check_data(m, dat, bmf(dfam ~ 1)), "no empty rows")
  # the aggregator's count and confidence checks
  long <- data.frame(
    id = 1, stimulus = 1,
    judgment = c("new", "know", "remember", "guess"),
    confidence = c(1, 2, 2, 2), count = c(10, 5, 7, 3)
  )
  expect_error(aggregate_sdt_cdp_data(replace(long, "count", c(10, -5, 7, 3)),
                                      "judgment", "confidence", "count"),
               "non-negative")
  expect_error(aggregate_sdt_cdp_data(replace(long, "count", c(10, NA, 7, 3)),
                                      "judgment", "confidence", "count"),
               "must not contain NA")
  expect_warning(
    wide <- aggregate_sdt_cdp_data(replace(long, "count", c(10, 2.5, 7, 3)),
                                   "judgment", "confidence", "count"),
    "rounded"
  )
  expect_equal(wide$know2, 2L)
  expect_error(aggregate_sdt_cdp_data(replace(long, "confidence", c(0, 2, 2, 2)),
                                      "judgment", "confidence", "count"),
               ">= 1")
  expect_error(aggregate_sdt_cdp_data(replace(long, "confidence", c(1, NA, 2, 2)),
                                      "judgment", "confidence", "count"),
               "must not contain NA")
})

test_that("the Stan cdp gradient is right where a threshold meets the strength mean and where rcrit equals kcrit", {
  # two measure-zero sets that init = 0 starts on: a lure criterion of 0 puts
  # a threshold exactly on mu_S, so the bivariate-normal CDF sees h = 0; and
  # rcrit == kcrit on lures makes all three Guess breakpoints coincide. The
  # values were always right there; the gradients were not (0 for the first,
  # -3.65 against a finite difference of -1.51 for the second)
  skip_unless_stan_tests()
  skip_on_cran()

  program <- paste0(
    "functions {\n",
    read_lines2(system.file("stan_chunks", "sdt_cdp_funs.stan", package = "bmm")),
    "\n}\n",
    "data { array[2, 12] int y; }\n",
    "parameters { real dfam; real drec; real criterion; real spacing;\n",
    "  real rcrit; real kcrit; real rho; }\n",
    "model {\n",
    "  array[0] real deltas;\n",
    "  vector[5] thr = cdp_make_thresholds(criterion, spacing, deltas, 3, 3, 1);\n",
    "  for (s in 1:2) {\n",
    "    vector[12] lp;\n",
    "    for (c in 1:12) lp[c] = log(cdp_category_prob(c, thr, dfam, drec,\n",
    "      0, rho, rcrit, kcrit, s - 1, 3, 3, 1));\n",
    "    y[s] ~ multinomial_logit(lp);\n",
    "  }\n}\n"
  )
  y <- rbind(c(60, 50, 40, 12, 10, 3, 8, 15, 20, 4, 10, 30),
             c(10, 15, 20, 10, 12, 5, 8, 20, 40, 10, 30, 90))
  model <- cmdstanr::cmdstan_model(cmdstanr::write_stan_file(program))
  ties <- list(
    threshold_on_mu_S = list(dfam = 0.8, drec = 1.0, criterion = 0, spacing = -0.3,
                             rcrit = 0.9, kcrit = 0.2, rho = 0),
    # lures have beta = 1/2 for every rho, so rcrit == kcrit makes all three
    # breakpoints coincide; the review measured the defect at rho = 1.5
    rcrit_equals_kcrit = list(dfam = 0.8, drec = 1.0, criterion = 0.2, spacing = -0.3,
                              rcrit = 0.9, kcrit = 0.9, rho = 1.5)
  )
  for (tie in names(ties)) {
    # cmdstan exits non-zero when a gradient misses its own error threshold;
    # the threshold is lifted so the comparison below is the expectation
    diag <- suppressMessages(model$diagnose(data = list(y = y), init = list(ties[[tie]]),
                                            error = 100))
    grad <- diag$gradients()
    expect_true(all(is.finite(grad$model)), info = tie)
    expect_lt(max(abs(grad$error)), 1e-5, label = paste(tie, "max |autodiff - finite difference|"))
  }
})