Two-Accumulator Linear Ballistic Accumulator (LBA) Model

The Linear Ballistic Accumulator (LBA) treats a choice as a race between two accumulators that rise linearly - no within-trial noise at all - each at a rate drawn afresh on every trial. The first to reach the threshold determines both the reaction time and the choice. The observed RT is that decision time shifted by a non-decision time ndt, and a fixed proportion poutlier of responses is generated by an outlier process instead of by the race.

Functions:

rcogmod_lba2(): Simulates random draws from the LBA.
dcogmod_lba2(): Computes the density (likelihood).
cogmod_lba2(): Creates a brms::custom_family() for use in brms models.
cogmod_lba2_stanvars(): Generates the stanvars to pass to brm().
p_outlier(): Per-trial posterior probability of being an outlier.

Usage

rcogmod_lba2(
  n,
  driftzero = 3,
  driftone = 3,
  sigmazero = 1,
  sigmaone = 1,
  sigmabias = 0.5,
  boundary = 0.5,
  ndt = 0.2,
  poutlier = 0
)

dcogmod_lba2(
  x,
  driftzero = 3,
  driftone = 3,
  sigmazero = 1,
  sigmaone = 1,
  sigmabias = 0.5,
  boundary = 0.5,
  ndt = 0.2,
  response,
  poutlier = 0,
  log = FALSE
)

cogmod_lba2(
  link_mu = "identity",
  link_driftone = "identity",
  link_sigmazero = "softplus",
  link_sigmaone = "softplus",
  link_sigmabias = "softplus",
  link_boundary = "softplus",
  link_ndt = "log",
  link_poutlier = "logit",
  predict_outliers = FALSE
)

cogmod_lba2_lpdf_expose()

cogmod_lba2_stanvars()

log_lik_cogmod_lba2(i, prep)

posterior_predict_cogmod_lba2(i, prep, predict_outliers = NULL, ...)

posterior_epred_cogmod_lba2(prep)

Arguments

n: Number of simulated trials. If length(n) > 1, the length is taken to be the number required.
driftzero, driftone: Mean drift rate of each accumulator (choice 0 and 1). Any real value; larger means faster. See Details on negative drifts.
sigmazero, sigmaone: Between-trial SD of each accumulator's drift rate. Must be positive.
sigmabias: Maximum starting point. The starting point of each accumulator on each trial is drawn from Uniform(0, sigmabias). Must be non-negative; 0 means both accumulators start at zero on every trial, so only the drift rates vary - the choice counterpart of the recinormal special case described in rcogmod_lba1().
boundary: Threshold offset, so the threshold is b = boundary + sigmabias. Must be positive.
ndt: Non-decision time (shift parameter), in seconds. Represents the time taken for processes unrelated to the decision (e.g., encoding, motor response). Must be non-negative. Range: [0, Inf).
poutlier: Proportion of responses generated by the outlier process rather than by the race. Range: [0, 1]. At poutlier = 0 the distribution reduces to the plain shifted LBA.
x: The observed reaction time (RT).
response: The winning accumulator (0 or 1). This gives the defective density f_response(x) * S_other(x), mixed with the outlier component, which is what a race likelihood needs.
log: Logical; if TRUE, returns the log-density. Default: FALSE.
link_mu, link_driftone: Link functions for the two mean drift rates. mu is driftzero: brms requires the first distributional parameter of a custom family to be called mu.
link_sigmazero, link_sigmaone: Link functions for the between-trial drift SDs.
link_sigmabias, link_boundary: Link functions for the start-point range and the threshold offset.
link_ndt, link_poutlier: Link functions for the non-decision time and the outlier rate.
predict_outliers: Logical; whether posterior_predict() should include the outlier component. FALSE (the default) fixes poutlier to zero for prediction, so predictions describe the race alone; the likelihood is always the full mixture either way. On the prediction method itself the default is NULL, which defers to the flag carried on the model - see with_outliers() to change it after fitting. See Details.
i, prep: For brms' functions to run: index of the observation and a brms preparation object.
...: Additional arguments.

Value

rcogmod_lba2() returns a data frame with n rows and two columns:

rt: The simulated reaction time.
response: The winning accumulator, coded 0 or 1, matching the dec() coding used by the brms families.

dcogmod_lba2() returns the defective density at each element of x - the log density if log = TRUE - for the response given in response, recycled to the length of the longest argument. cogmod_lba2() returns a brms::custom_family object, to put on a brms::bf() formula. cogmod_lba2_stanvars() returns a brms::stanvars object holding the family's Stan functions block, to pass to brms::brm(), and cogmod_lba2_lpdf_expose() compiles that Stan code and returns it as an R function, for checking the density outside of a model. The remaining functions are brms post-processing methods, called by brms rather than directly: log_lik_cogmod_lba2() returns a numeric vector holding one log-likelihood value per posterior draw for observation i, and posterior_predict_cogmod_lba2() a draws x 2 matrix of reaction times and choices simulated for observation i. posterior_epred_cogmod_lba2() returns nothing: the expected reaction time of a race has no closed form, so it errors rather than report one - summarise posterior_predict() draws instead.

Parameterization

Accumulator k starts at a point z ~ Uniform(0, sigmabias), drawn afresh on every trial, and rises at a constant rate v ~ Normal(drift_k, sigma_k), also drawn afresh on every trial, until it reaches the threshold b = boundary + sigmabias. Its finishing time is therefore (b - z) / v, and boundary is the threshold offset: the distance from the highest possible starting point to the threshold. All the randomness is between trials; within a trial the path is a straight line, which is what makes the density closed-form.

This is the B parameterization of DMC and EMC2, where b = B + A with A the start-point range. It is used for a reason rather than for taste: the threshold has to sit above the highest possible starting point, and writing the offset makes b > A hold automatically for any positive value. The alternative - estimating the absolute threshold, as rtdists does - needs an order constraint between two estimated parameters, which has to hold in every cell of the design once either of them carries a predictor. The cost is that boundary alone is not the quantity to read off a fitted model; boundary + sigmabias is.

ndt is expressed directly, in seconds (through a log link in the brms family). Nothing about it is taken from the data: it is not bounded by the fastest observed response, so a non-decision time that varies by condition or by participant can exceed the sample minimum wherever the data support it.

Tying ndt to the fastest observed response would cap it at an order statistic of the sample, so any condition or participant whose true ndt exceeded that response would be inexpressible, and the misfit would surface as spurious effects on the race parameters. Expressing it directly is what avoids that.

Negative drift rates

A Normal drift rate can come out negative, and such an accumulator rises away from the threshold and never responds. Brown and Heathcote (2008) noted the problem and left it; every implementation since has had to decide what to do about it, and this one follows the convention of rtdists (posdrift = TRUE, its default), DMC, EMC2 and ggdmc: each drift rate is a Normal truncated at zero. Every accumulator is then guaranteed a positive rate on every trial, every trial produces a response, and the losing accumulator's survival is that of a truncated Normal. rcogmod_lba2() draws each drift from exactly that truncated distribution.

The density is normalised to match: the winner's defective density is divided by its own pnorm(drift / sigma), and the loser's survival is taken conditional on its own drift being positive. Without the truncation the density integrates to the probability that at least one accumulator finishes rather than to one, which at low drift rates is a long way short - 0.83 at drifts of 0.5 and 0.2 with SDs of 1.5, so the likelihood would be wrong by 17% and wrong by different amounts at different parameter values, which is what biases estimates rather than merely offsetting them.

Two consequences are worth knowing. First, because the convention is the field's, parameter estimates are directly comparable with those packages' and with the published LBA literature built on them, and dcogmod_lba2() reproduces rtdists::dLBA() at the same parameter values. Second, drift and sigma are the location and scale of the untruncated Normal, not the mean and SD of the drifts actually realised: where a drift is small relative to its SD, the realised mean is higher and the realised SD lower than the parameters say. That is the price every implementation pays for a race that always finishes; with both drifts large relative to their SDs it is negligible.

The truncation also creates a flat direction. A Normal truncated at zero whose location runs off to minus infinity while its scale grows, with |drift| / sigma^2 held fixed, converges to an Exponential, so once an accumulator rarely wins its drift and sigma are identified only through that ratio, and the likelihood is nearly flat along the ray. The error accumulator in a task with a 5% error rate is exactly such a case: left flat, its drift wandered to -12 with an interval of -23 to -6.5, on an evidence scale where the correct accumulator's drift is 3. cogmod_priors() therefore puts normal(1, 2) on driftone and normal(0, 1.5) on its slopes, as it does for cogmod_lnr()'s nuone. If the rarely chosen option is the one on mu, mirror that prior onto mu by hand. Fixing both SDs (sigmazero = 1, sigmaone = 1, a single sv) removes the ray altogether and is common practice in the LBA literature.

Versions before 0.3.1 used a different convention - the pair of drifts was conditioned on at least one being positive, with a negative loser allowed to run away - so LBA estimates from those versions are not comparable with these where a drift is small relative to its SD, and fits made with them cannot be post-processed with this version.

The evidence scale is arbitrary

Multiply driftzero, driftone, sigmazero, sigmaone, sigmabias and boundary all by any c > 0 and every finishing time (b - z) / v is unchanged. The likelihood is therefore exactly constant along that ray, which runs to infinity in both directions: the six parameters are identified only up to a common scale factor.

cogmod_priors() puts priors on all four positive parameters, which makes the posterior proper and the sampler well behaved, but a prior does not identify a direction the likelihood cannot see. If the individual parameters are to be interpreted - rather than the RT distribution they jointly generate, which is perfectly well identified - fix one SD in the formula, the usual convention being sigmazero = 1:

f <- brms::bf(RT | dec(Error) ~ Condition, driftone ~ Condition,
              sigmazero = 1, sigmaone ~ 1, sigmabias ~ 1, boundary ~ 1,
              ndt ~ 1, poutlier ~ 1, family = cogmod_lba2())

A second, milder flat direction remains either way: sigmabias and boundary enter only through the sum b = boundary + sigmabias, so they trade off along a ridge. The sum is the trustworthy quantity to interpret and to compare across conditions. cogmod_lba1() and cogmod_rdm() share it.

The outlier component

A shifted distribution assigns exactly zero density to any response faster than ndt, which puts a hard boundary in the likelihood at the fastest observed RT. Mixing in a component with support over the whole positive line removes it: every response keeps positive density whatever ndt is, so the boundary becomes a finite cost rather than a wall and the log-density stays smooth and differentiable.

Because this model produces a choice as well as a time, the contaminant has to produce both. It is a guess: the choice is uniform over the two options, and the RT is a half Normal with scale 0.2 seconds.

$$f(t, k) = p \frac{1}{K} g(t) + (1 - p) f_k(t - ndt)$$

The 1 / K is what keeps the total summing to one over the response options; without it it would come to 1 + poutlier.

poutlier is a rate, not a classification: the model never labels individual trials, it estimates what share of them came from elsewhere. Use p_outlier() for per-trial posterior probabilities.

Reaction times must be in seconds

The outlier component's scale is a constant in seconds, and so are the priors cogmod_priors() supplies. There is no argument for changing the unit. Millisecond data fails silently rather than loudly - the outlier component contributes nothing and the min-RT boundary comes back. See the corresponding section of cogmod_lognormal() for the full account, which applies unchanged here.

Fitting

f <- brms::bf(RT | dec(Error) ~ Condition, driftone ~ Condition,
              sigmazero = 1, sigmaone ~ 1, sigmabias ~ 1, boundary ~ 1,
              ndt ~ 1, poutlier ~ 1, family = cogmod_lba2())
brms::brm(f, data = df,
          prior    = cogmod_priors(f, df),
          init     = cogmod_inits(f, df),
          stanvars = cogmod_stanvars(f))

The brms family names the drift of the first accumulator mu (as brms requires) and that of the second driftone. Use cogmod_inits() rather than init = 0: brms initialises on the unconstrained scale, so init = 0 puts ndt at exp(0) = 1 second - above nearly every sub-second RT, which leaves every response attributed to the outlier component and the race parameters with no gradient at all.

Predictions exclude the outlier component

posterior_predict() describes the race alone by default, as if poutlier were zero, because the outlier component is a fixed regularizer rather than a claim about how guesses are distributed. Use with_outliers() for the fitted mixture - chiefly for brms::pp_check() - and without_outliers() to go back. log_lik() is always the full mixture.

posterior_epred() is not provided: for a race model the expectation needs numerical integration per draw and per observation, and users are better off summarising posterior_predict() draws.

References

Brown, S. D., & Heathcote, A. (2008). The simplest complete model of choice response time: Linear ballistic accumulation. Cognitive Psychology, 57(3), 153-178. doi:10.1016/j.cogpsych.2007.12.002

Examples

# Simulate data, with 2% of trials from the outlier process
data <- rcogmod_lba2(1000,
  driftzero = 3, driftone = 2, sigmazero = 1, sigmaone = 1,
  sigmabias = 0.5, boundary = 0.5, ndt = 0.2, poutlier = 0.02
)
head(data)
#>          rt response
#> 1 0.4344558        0
#> 2 0.6071895        1
#> 3 0.4257365        1
#> 4 0.6650784        0
#> 5 0.3747992        1
#> 6 0.4242575        0

# Responses faster than ndt keep positive density, unlike the unmixed model
dcogmod_lba2(0.1, ndt = 0.2, response = 0, poutlier = 0.02)
#> [1] 0.03520653
dcogmod_lba2(0.1, ndt = 0.2, response = 0, poutlier = 0)
#> [1] 0

if (FALSE) { # \dontrun{
# Needs cmdstanr and a CmdStan toolchain, which live outside CRAN - see the
# package website to install them. Not run under R CMD check, which executes
# every example in one R session: once brms has fitted a model there (the
# cogmod_inits() and p_outlier() examples do), rstan is live in the process
# and loading an exposed Stan function next to it segfaults on Linux.
lpdf <- cogmod_lba2_lpdf_expose()
lpdf(
  Y = 0.5, mu = 3, driftone = 2, sigmazero = 1, sigmaone = 1,
  sigmabias = 0.5, boundary = 0.5, ndt = 0.2, poutlier = 0.02, dec = 0
)
} # }