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") |>
  set_parameters(separate = list(m = TRUE, s = TRUE)) |>
  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(
  trials = sim$trials,
  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,
    n_trials = n_trials
  )
  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)
)

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))