Skip to contents

Pooling Accumulators

Source: vignettes/pool_model.Rmd
pool_model.Rmd

This vignette shows how to define a response from a pool of accumulators. Here response A is produced by a 2-of-3 pool: the response is observed once any two of A1, A2, and A3 have finished. Response B is generated by a single competing accumulator.

## 
## Attaching package: 'AccumulatR'
## The following object is masked from 'package:stats':
## 
##     simulate

Define the model We use three lognormal accumulators feeding a pooled response A, plus one lognormal accumulator feeding response B.

model <- race_spec() |>
  add_accumulator("A1", "lognormal") |>
  add_accumulator("A2", "lognormal") |>
  add_accumulator("A3", "lognormal") |>
  add_accumulator("B", "lognormal") |>
  add_pool("A_pool", c("A1", "A2", "A3"), k = 2L) |>
  add_outcome("A", "A_pool") |>
  add_outcome("B", "B") |>
  finalize_model()

true_params <- c(
  A1.m = log(0.28),
  A1.s = 0.16,
  A2.m = log(0.28),
  A2.s = 0.16,
  A3.m = log(0.28),
  A3.s = 0.16,
  B.m = log(0.28),
  B.s = 0.18
)

Simulate data We simulate a response and response time for each trial.

set.seed(123456)

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

sim <- simulate(model, params_df)

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

table(data_df$R)
## 
##   A   B 
## 751 749

When response A is observed, its response time is the finishing time of the second finisher of the pool.

Estimate parameters with optim() We estimate a shared location and spread for the three pool members, plus the location and spread of accumulator B. The spread parameters are optimized on the log scale.

prepared <- prepare_data(model, data_df)
ctx <- make_context(model)
neg_loglik <- function(theta) {
  est <- true_params
  est[c("A1.m", "A2.m", "A3.m")] <- theta[["A_pool.m"]]
  est[c("A1.s", "A2.s", "A3.s")] <- exp(theta[["log_A_pool.s"]])
  est["B.m"] <- theta[["B.m"]]
  est["B.s"] <- exp(theta[["log_B.s"]])
  params_df <- build_param_matrix(
    model,
    est,
    trial_df = prepared
  )
  ll <- log_likelihood(ctx, prepared, params_df)
  -as.numeric(ll)
}

start <- c(
  A_pool.m = log(0.24),
  log_A_pool.s = log(0.12),
  B.m = log(0.34),
  log_B.s = log(0.12)
)

set.seed(123456)
fit <- optim(
  start,
  neg_loglik,
  method = "Nelder-Mead",
  control = list(maxit = 4000, reltol = 1e-9)
)
fit_params <- c(
  A_pool.m = fit$par[["A_pool.m"]],
  A_pool.s = exp(fit$par[["log_A_pool.s"]]),
  B.m = fit$par[["B.m"]],
  B.s = exp(fit$par[["log_B.s"]])
)
target <- c(
  A_pool.m = true_params[["A1.m"]],
  A_pool.s = true_params[["A1.s"]],
  B.m = true_params[["B.m"]],
  B.s = true_params[["B.s"]]
)

data.frame(true = target, recovered = fit_params, miss = abs(target - fit_params))
##               true  recovered         miss
## A_pool.m -1.272966 -1.2749227 0.0019570346
## A_pool.s  0.160000  0.1605954 0.0005954359
## B.m      -1.272966 -1.2731653 0.0001996564
## B.s       0.180000  0.1768839 0.0031160588

The code path is the same as in a basic race model. The difference is that the observed response A now depends on a pooled completion rule rather than a single accumulator.