Multiple Outcomes

This vignette shows how to work with data in which more than one ordered response is observed on each trial. Here each trial records the first response (R, rt) and the second response (R2, rt2), and so on for any number of responses (R3, rt3 etc.).

library(AccumulatR)

Define the model

We use a simple two-accumulator race with direct responses A and B. Setting n_outcomes = 2 tells the model to retain the first and second finishing responses.

model <- race_spec(n_outcomes = 2L) |>
  add_accumulator("A", "lognormal") |>
  add_accumulator("B", "lognormal") |>
  add_outcome("A", "A") |>
  add_outcome("B", "B") |>
  set_parameters(separate = list(m = TRUE, s = TRUE)) |>
  finalize_model()

true_params <- c(
  A.m = log(0.30),
  A.s = 0.18,
  B.m = log(0.38),
  B.s = 0.22
)

Simulate data

We generate data and keep both ordered responses for each trial.

set.seed(123456)

n_trials <- 500
params_df <- build_param_matrix(model, true_params, n_trials = n_trials)

sim <- simulate(model, params_df)

data_df <- data.frame(
  trials = sim$trials,
  R = factor(sim$R),
  rt = sim$rt,
  R2 = factor(sim$R2),
  rt2 = sim$rt2,
  stringsAsFactors = FALSE
)

head(data_df)
##   trials R        rt R2       rt2
## 1      1 A 0.3394130  B 0.3514512
## 2      2 A 0.2373401  B 0.4565748
## 3      3 A 0.3709420  B 0.4913625
## 4      4 B 0.2974072  A 0.3441526
## 5      5 A 0.3297834  B 0.4790856
## 6      6 A 0.3671582  B 0.4663272

Estimate parameters with optim()

We estimate A.m, A.s, B.m, and B.s. The lognormal standard deviations are optimized on the log scale and transformed back inside the function.

prepared <- prepare_data(model, data_df)
ctx <- make_context(model)
neg_loglik <- function(theta) {
  est <- true_params
  est[c("A.m", "A.s", "B.m", "B.s")] <- theta[c("A.m", "A.s", "B.m", "B.s")]
  est[c("A.s", "B.s")] <- exp(est[c("A.s", "B.s")])
  params_df <- build_param_matrix(
    model,
    est,
    n_trials = n_trials
  )
  ll <- log_likelihood(ctx, prepared, params_df)
  -as.numeric(ll)
}

start <- c(
  A.m = log(0.24),
  A.s = log(0.12),
  B.m = log(0.48),
  B.s = log(0.12)
)

fit <- optim(start, neg_loglik, method = "Nelder-Mead")

fit_params <- fit$par
fit_params[c("A.s", "B.s")] <- exp(fit_params[c("A.s", "B.s")])
target <- true_params[c("A.m", "A.s", "B.m", "B.s")]

data.frame(
  true = target,
  recovered = fit_params,
  miss = abs(target - fit_params)
)

Ranked observations require direct accumulator or pool outcomes. Chained onsets are supported, but logical outcome expressions, guessing, remapping, censoring, and truncation are not. Observed response times must increase strictly and response labels must be distinct. Later unobserved ranks use NA in both columns; the first response must be observed.