Two-Accumulator Racing Diffusion Model (RDM)

The Racing Diffusion Model (RDM) treats a choice as a race between two diffusion processes, one per response option, each accumulating evidence at its own rate until it reaches a common threshold. The winner determines both the observed 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_rdm(): Simulates random draws from the RDM.
dcogmod_rdm(): Computes the density (likelihood).
pcogmod_rdm(): Computes the CDF of the reaction time, marginally over the choice or defectively for one response.
qcogmod_rdm(): Computes the corresponding quantiles.
cogmod_rdm(): Creates a brms::custom_family() for use in brms models.
cogmod_rdm_stanvars(): Generates the stanvars to pass to brm().
p_outlier(): Per-trial posterior probability of being an outlier.

Usage

rcogmod_rdm(
  n,
  vzero = 3,
  vone = 2,
  boundary = 0.5,
  bias = 0.2,
  ndt = 0.2,
  poutlier = 0
)

dcogmod_rdm(
  x,
  vzero = 3,
  vone = 2,
  boundary = 0.5,
  bias = 0.2,
  ndt = 0.2,
  response = NULL,
  poutlier = 0,
  log = FALSE
)

pcogmod_rdm(
  q,
  vzero = 3,
  vone = 2,
  boundary = 0.5,
  bias = 0.2,
  ndt = 0.2,
  poutlier = 0,
  response = NULL,
  lower.tail = TRUE,
  log.p = FALSE
)

qcogmod_rdm(
  p,
  vzero = 3,
  vone = 2,
  boundary = 0.5,
  bias = 0.2,
  ndt = 0.2,
  poutlier = 0,
  response = NULL,
  scale_p = FALSE,
  lower.tail = TRUE,
  log.p = FALSE,
  interval = c(0, 10)
)

cogmod_rdm(
  link_mu = "softplus",
  link_driftone = "softplus",
  link_sigmabias = "softplus",
  link_boundary = "softplus",
  link_ndt = "log",
  link_poutlier = "logit",
  predict_outliers = FALSE
)

cogmod_rdm_lpdf_expose()

cogmod_rdm_stanvars()

log_lik_cogmod_rdm(i, prep)

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

posterior_epred_cogmod_rdm(prep)

Arguments

n: Number of simulated trials. If length(n) > 1, the length is taken to be the number required.
vzero, vone: Drift rates of the two accumulators (choice 0 and 1). Must be non-negative; larger means faster. Zero is allowed - such an accumulator is slow, but it still finishes, and can still win. Range: [0, Inf).
boundary: Threshold offset, boundary = b - bias, where b is the decision threshold and bias the maximum starting point. Must be positive.
bias: Maximum starting point. The starting point of each accumulator on each trial is drawn from Uniform(0, bias). Must be non-negative; zero is allowed and gives the plain Wald race, in which both accumulators start at 0 on every trial. Range: [0, Inf). Called sigmabias in the brms family, to match cogmod_lba2().
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 RDM.
x: The observed reaction time (RT).
response: Accumulator whose finishing time is being scored: 0 for the vzero accumulator, 1 for the vone accumulator. This gives the defective density f_response(x) * S_other(x), mixed with the outlier component, which is what a race likelihood needs. The default NULL instead returns the marginal density of the RT, ignoring which accumulator won - the sum of the two.
log: Logical; if TRUE, returns the log-density. Default: FALSE.
q: Vector of quantiles (reaction times).
lower.tail: If TRUE (default) return P(RT <= q), otherwise the survival P(RT > q). With a response, both are defective - see Details.
log.p: If TRUE, probabilities are returned on the log scale.
p: Vector of probabilities. With response = NULL these are ordinary probabilities of the marginal RT distribution. With a response they are read off the defective CDF unless scale_p = TRUE, so they must be below the probability of that response; anything above it has no quantile and comes back NA with a warning.
scale_p: Logical. If TRUE, p is taken as a fraction of the chosen response's own probability rather than of the whole distribution, so that p = 0.5 is that response's median. This is what a quantile-probability plot wants. Ignored when response is NULL. Default FALSE.
interval: Length-2 numeric giving the initial bracket, in seconds, for the root search. The upper end is doubled until it covers the requested probability, so this only affects speed.
link_mu, link_driftone: Link functions for the two drift rates. mu is vzero: brms requires the first distributional parameter of a custom family to be called mu, so that is the name the formula and this argument use, and the drift of accumulator 0 is what it means.
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_rdm() returns a data frame with n rows and two columns:

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

dcogmod_rdm() returns the density at each element of x - the log density if log = TRUE - pcogmod_rdm() the cumulative probability at each element of q, and qcogmod_rdm() the quantile at each element of p, in seconds. With a response the latter two are defective, i.e. scaled to that response's own probability rather than to one. All are numeric vectors, recycled to the length of the longest argument. cogmod_rdm() returns a brms::custom_family object, to put on a brms::bf() formula. cogmod_rdm_stanvars() returns a brms::stanvars object holding the family's Stan functions block, to pass to brms::brm(), and cogmod_rdm_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_rdm() returns a numeric vector holding one log-likelihood value per posterior draw for observation i, and posterior_predict_cogmod_rdm() a draws x 2 matrix of reaction times and choices simulated for observation i. posterior_epred_cogmod_rdm() 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.

Details

pcogmod_rdm() with response = NULL (the default) describes the RT of the trial as a whole - whichever accumulator wins, and whether or not the trial came from the outlier component - since P(min(T0, T1) > q) = S0(q) * S1(q). That is a closed form and is exact.

With a response, it returns the defective CDF P(RT <= q, choice = response), which is what a defective-CDF or quantile-probability plot needs. It does not reach one: its limit is the probability of that response, which pcogmod_rdm(Inf, response = k) gives. There is no closed form for it, so it is obtained by quadrature over the defective density: accurate to about 1e-8 rather than to machine precision, and about ten times slower per element (roughly 3 ms against 0.3 ms), since the marginal is a vectorised closed form and this is a loop. lower.tail = FALSE integrates the upper side directly rather than subtracting, so the defective survival stays accurate into the tail.

qcogmod_rdm() inverts pcogmod_rdm() by root-finding, and so inherits its quadrature error where a response is given. It is the natural way to get the RT quantiles of each response for a quantile-probability plot: ask for p = c(0.1, 0.3, 0.5, 0.7, 0.9) with scale_p = TRUE, once per response.

Parameterization

Each accumulator is a diffusion with drift rate v and unit diffusion coefficient, starting from a point z ~ Uniform(0, bias) drawn afresh on every trial and finishing when it reaches the threshold b = boundary + bias. The distance it has to cover is therefore b - z = boundary + bias - z, which is what makes boundary the threshold offset: the distance from the highest possible starting point to the threshold. Its first passage time is Wald (inverse Gaussian) with the start point integrated out. The observed reaction time is ndt + min(T_0, T_1) and the observed choice is whichever accumulator got there first.

This is the B parameterization of DMC and EMC2, where b = B + A with A the start-point range (bias here). 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 + bias is.

A drift rate of exactly zero is allowed, and is not the same as an accumulator that never responds: driftless Brownian motion still reaches any positive level with probability one, so a zero-drift accumulator is slow but still finishes, and can still win the race.

A start-point range of exactly zero is allowed too, and is a model rather than a degenerate parameter: both accumulators then start at 0 on every trial and the race is between two plain Walds - equation 2 of Tillman et al. (2020), which is the limit the density already takes. cogmod_lba1() and cogmod_lba2() have always allowed it. Note that this does not make the sigmabias direction any better identified - see Fitting below.

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.

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. That is what makes the direct parameterization of ndt workable without taking a bound from the data.

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. The half Normal is used for the timing because it is flat at the origin (zero derivative), so the very fastest responses - the ones least plausibly decisions - are not starved of density, and because it dies fast enough above that range to leave the slow tail to the race itself rather than claiming it.

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,
              sigmabias ~ 1, boundary ~ 1, ndt ~ 1, poutlier ~ 1,
              family = cogmod_rdm())
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, and calls the start-point range sigmabias to match cogmod_lba2(), where it denotes the same quantity. Note that this is not the same thing as bias in cogmod_ddm(), which is a relative starting point in [0, 1]. Both drifts use a softplus link with a lower bound of zero, following cogmod_invgaussian(): a Wald drift must be non-negative for the accumulator to be a proper first passage time.

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.

cogmod_priors() is not a convenience here either. Beyond ndt and poutlier, sigmabias and boundary are only weakly identified from each other, because they enter the threshold only through the sum b = boundary + sigmabias and trade off almost freely: on simulated data with 4000 trials the profile log-likelihood varies by only about 3 units as sigmabias ranges from 0 to half the threshold, while boundary slides to compensate. With flat priors the sampler tends to wander down the sigmabias -> 0 ridge (the plain Wald race) and produce divergent transitions, and a softplus link reaches zero only at minus infinity - a flat prior over an unbounded flat region, which is an improper posterior. That the endpoint is now a legal parameter value does not help: the link never reaches it, so the flat direction is as long as it ever was. cogmod_priors() fences both off, exactly as it does for cogmod_lba1(), which shares this parameterisation.

The sum boundary + sigmabias is well identified either way, so it is the more trustworthy quantity to interpret and to compare across conditions. The same caveat applies to cogmod_lba2().

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

Michael, J. R., Schucany, W. R., & Haas, R. W. (1976). Generating Random Variates Using Transformations with Multiple Roots. The American Statistician, 30(2), 88-90. doi:10.2307/2683801
Tillman, G., Van Zandt, T., & Logan, G. D. (2020). Sequential sampling models without random between-trial variability: The racing diffusion model of speeded decision making. Psychonomic Bulletin & Review, 27, 911-936. doi:10.3758/s13423-020-01719-6
Folks, J. L., & Chhikara, R. S. (1978). The inverse Gaussian distribution and its statistical application-a review. Journal of the Royal Statistical Society Series B: Statistical Methodology, 40(3), 263-275.

Examples

# Simulate data, with 2% of trials from the outlier process
data <- rcogmod_rdm(1000,
  vzero = 2.5, vone = 1.6, boundary = 0.5, bias = 0.2,
  ndt = 0.2, poutlier = 0.02
)
head(data)
#>          rt response
#> 1 0.3342623        0
#> 2 0.3085687        0
#> 3 0.2719745        0
#> 4 0.4592891        0
#> 5 0.5737855        0
#> 6 0.3282243        1

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

# Defective CDF of one response: at q = Inf it is that response's probability
pcogmod_rdm(c(0.4, 0.6, Inf), vzero = 2.5, vone = 1.6, response = 0)
#> [1] 0.4449708 0.5822743 0.6112906

# The RT quantiles of each response, for a quantile-probability plot
sapply(0:1, function(k) {
  qcogmod_rdm(c(0.1, 0.3, 0.5, 0.7, 0.9),
    vzero = 2.5, vone = 1.6, response = k, scale_p = TRUE
  )
})
#>           [,1]      [,2]
#> [1,] 0.2604085 0.2601144
#> [2,] 0.2961675 0.2951968
#> [3,] 0.3345118 0.3325338
#> [4,] 0.3894861 0.3857909
#> [5,] 0.5113970 0.5034841

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_rdm_lpdf_expose()
lpdf(
  Y = 0.5, mu = 2, driftone = 1.5, sigmabias = 0.2, boundary = 0.5,
  ndt = 0.2, poutlier = 0.02, dec = 0
)
} # }