Drift Diffusion Model (DDM)

The Drift Diffusion Model (DDM) describes a two-choice decision as noisy evidence accumulating between two boundaries until one of them is reached. The boundary reached is the choice and the time taken is the decision time. 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 diffusion.

Functions:

• rcogmod_ddm(): Simulates random draws from the DDM.

• dcogmod_ddm(): Computes the density (likelihood).

• pcogmod_ddm(): Computes the cumulative distribution function.

• cogmod_ddm(): Creates a brms::custom_family() for use in brms models.

• cogmod_ddm_stanvars(): Generates the stanvars to pass to brm().

• p_outlier(): Per-trial posterior probability of being an outlier.

pcogmod_ddm() is the cumulative distribution function: the probability that a response has been made by time q. With response = NULL it is the RT distribution marginally over the choice; with response set it is the defective CDF that boundary carries, rising to the probability of that boundary rather than to one.

Usage

rcogmod_ddm(
  n,
  drift = 0,
  boundary = 1,
  bias = 0.5,
  ndt = 0.2,
  sigmadrift = 0,
  sigmabias = 0,
  sigmandt = 0,
  poutlier = 0
)

dcogmod_ddm(
  x,
  drift = 0,
  boundary = 1,
  bias = 0.5,
  ndt = 0.2,
  response,
  sigmadrift = 0,
  sigmabias = 0,
  sigmandt = 0,
  poutlier = 0,
  log = FALSE
)

pcogmod_ddm(
  q,
  drift = 0,
  boundary = 1,
  bias = 0.5,
  ndt = 0.2,
  response = NULL,
  poutlier = 0,
  lower.tail = TRUE,
  log.p = FALSE
)

cogmod_ddm(
  link_mu = "identity",
  link_boundary = "softplus",
  link_bias = "logit",
  link_sigmadrift = "softplus",
  link_sigmabias = "logit",
  link_sigmandt = "log",
  link_ndt = "log",
  link_poutlier = "logit",
  predict_outliers = FALSE
)

cogmod_ddm_lpdf_expose()

cogmod_ddm_stanvars()

log_lik_cogmod_ddm(i, prep)

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

posterior_epred_cogmod_ddm(prep, predict_outliers = NULL)

Arguments

n

Number of simulated trials. If length(n) > 1, the length is taken to be the number required.

drift

Drift rate. Any real value; positive pushes the accumulator towards the boundary coded 1.

boundary

Boundary separation. Must be positive.

bias

Starting point, as a proportion of the boundary separation measured from the boundary coded 0. Must be in (0, 1).

ndt

Non-decision time (shift parameter), in seconds. Must be non-negative. With sigmandt > 0 it is the lower bound of the between-trial distribution rather than its midpoint.

sigmadrift

Between-trial SD of the drift rate (sv). Must be non-negative. Default 0.

sigmabias

Between-trial start-point range, as a fraction in [0, 1) of the widest range that keeps the start point inside the boundaries: sw = sigmabias * min(2 * bias, 2 * (1 - bias)). Default 0.

sigmandt

Between-trial range of the non-decision time (st0), in the same unit as the data, with ndt its lower bound. Default 0.

poutlier

Proportion of responses generated by the outlier process rather than by the diffusion. Range: [0, 1].

x

The observed reaction time (RT).

response

The boundary reached: 1 for the upper boundary, 0 for the lower one. This gives the defective density that boundary carries, mixed with the outlier component.

log

Logical; if TRUE, returns the log-density. Default: FALSE.

q

Vector of quantiles (reaction times, in seconds).

lower.tail

Logical; if TRUE (default) the probability is P(RT <= q), otherwise P(RT > q). With a response, both are defective

• see Details.

log.p

Logical; if TRUE, probabilities are returned on the log scale.

link_mu, link_boundary, link_bias

Link functions for the drift rate, the boundary separation and the starting point. mu is the drift: brms requires the first distributional parameter of a custom family to be called mu.

link_sigmadrift, link_sigmabias, link_sigmandt

Link functions for the between-trial variability parameters. Fix them in the brms::bf() formula (e.g. sigmadrift = 0) to recover the classic 4-parameter DDM.

link_ndt, link_poutlier

Link functions for the non-decision time and the outlier rate.

predict_outliers

Logical; whether posterior_predict() and posterior_epred() should include the outlier component. FALSE (the default) fixes poutlier to zero for prediction, so predictions describe the diffusion alone; the likelihood is always the full mixture either way. On the prediction methods themselves the default is NULL, which defers to the flag carried on the model - see with_outliers(). See Details.

i, prep

For brms' functions to run: index of the observation and a brms preparation object.

...

Additional arguments.

Value

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

rt

The simulated reaction time.

response

The boundary reached, 1 for upper and 0 for lower, matching the dec() coding used by the brms families.

dcogmod_ddm() returns the defective density at each element of x - the log density if log = TRUE - and pcogmod_ddm() the defective cumulative probability at each element of q, both for the response given in response and recycled to the length of the longest argument. cogmod_ddm() returns a brms::custom_family object, to put on a brms::bf() formula. cogmod_ddm_stanvars() returns a brms::stanvars object holding the family's Stan functions block, to pass to brms::brm(), and cogmod_ddm_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_ddm() returns a numeric vector holding one log-likelihood value per posterior draw for observation i, posterior_predict_cogmod_ddm() a draws x 2 matrix of reaction times and choices simulated for observation i, and posterior_epred_cogmod_ddm() a draws x observations matrix of expected reaction times (marginal over the two responses, and only approximate once the between-trial variability parameters are non-zero).

Response coding

The response coded 1 (response here, dec() in a brms formula) is the upper boundary and the response coded 0 is the lower one, following brms's own wiener() family. Two consequences are worth keeping in mind when reading a fitted model:

• bias is measured from the lower boundary, so bias > 0.5 places the starting point closer to the response coded 1, and bias < 0.5 closer to the response coded 0. Since first-passage times are shorter for the nearer boundary, bias > 0.5 makes responses coded 1 faster than responses coded 0, and bias < 0.5 makes them slower. At exactly bias = 0.5 the two conditional RT distributions are identical, whatever the drift rate.

• A positive drift pushes the accumulator towards the response coded 1. When 0 codes correct responses and 1 codes errors (a common choice when the task has no natural stimulus-to-boundary mapping), good performance therefore corresponds to a negative drift rate.

Parameterization

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 other parameters. Expressing it directly is what avoids that.

sigmandt is the between-trial range of the non-decision time (st0 in the usual notation), expressed directly in the same unit as the data, with ndt the lower bound of the resulting Uniform.

Between-trial variability

sigmadrift, sigmabias and sigmandt extend the classic 4-parameter process to the full 7-parameter one. Each is legitimately zero, and setting all three to zero in the formula recovers the classic model:

brms::bf(RT | dec(Error) ~ Condition, sigmadrift = 0, sigmabias = 0,
         sigmandt = 0, ndt ~ 1, poutlier ~ 1, family = cogmod_ddm())

Writing sigmadrift = 0 fixes the parameter; leaving it out of bf() altogether estimates it, which is not the same thing. All three are hard to recover even from a lot of data, and each has a flat direction at its own floor - the link only reaches zero at minus infinity, and the likelihood stops changing well before then - so cogmod_priors() gives all three deliberately tight priors. Fixing the ones a design cannot identify is usually better than estimating them behind a prior.

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.

Implementation

The 4-parameter density is the Navarro and Fuss (2009) series, evaluated in log space and vectorised over parameter sets, so a response in the far tail has a finite log-density rather than log(0). It agrees with brms::dwiener() and with rtdists to about 1e-10 wherever those return a number, and is about eight times cheaper per element than the former. Draws are taken by inverting the first-passage CDF, which - unlike brms::rwiener(), and unlike rtdists - vectorises over parameter sets, so the per-call setup is paid once rather than once per posterior draw. That matters for rstantools::posterior_predict(), where every draw carries its own parameters; it is several times faster there and agrees with both packages' samplers to within sampling error.

The full 7-parameter model is built on top of these rather than delegated to another package: between-trial variability is simulated by drawing the per-trial parameters, and evaluated by combining a closed-form drift correction with Gauss-Legendre quadrature over the starting point and non-decision time, also in log space.

In Stan the decision component is wiener_lpdf(), called with its own non-decision time set to zero because the shift is applied by the mixture around it, and short-circuited to the cheaper 4- and 5-parameter forms whenever the start-point and non-decision-time ranges both vanish. The cheapest of the three, Stan's classic 4-parameter density, is used only where it is sound: it returns -inf with NaN derivatives once the rescaled decision time t / boundary^2 falls below about 6.6e-4, or once the density underflows, and a NaN derivative on a trial the mixture gives no weight to still turns the gradient of the whole model to NaN. Those calls go to the sv-capable form instead, which the two agree with to 1e-13 where they meet.

Fitting

f <- brms::bf(RT | dec(Error) ~ Condition, boundary ~ Condition, bias ~ 1,
              sigmadrift = 0, sigmabias = 0, sigmandt = 0,
              ndt ~ 1, poutlier ~ 1, family = cogmod_ddm())
brms::brm(f, data = df,
          prior    = cogmod_priors(f, df),
          init     = cogmod_inits(f, df),
          stanvars = cogmod_stanvars(f))

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 diffusion parameters with no gradient at all.

Predictions exclude the outlier component

posterior_predict() and posterior_epred() describe the diffusion 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.

Accuracy of pcogmod_ddm()

This CDF is more accurate than the one in rtdists, which is the usual reference implementation. Deviation from numerical integration of the density, at drift = -4, boundary = 0.6, bias = 0.25, by decision time:

decision time	this function	rtdists::pdiffusion()
0.02	+2.8e-16	+5.5e-04
0.05	-1.1e-16	+4.0e-04
0.30	-1.1e-16	+1.5e-04

Over a grid of 360 (drift, boundary, bias, boundary-reached, time) cells the worst deviation from integrating dcogmod_ddm() is 1e-15, i.e. rounding. The choice probabilities are exact to the same order, where rtdists::pdiffusion(Inf, ...) is out by up to 1.4e-04. The difference is not academic: this function exists because rcogmod_ddm() inverts it to draw from, so any error in it would land directly in the draws.

pcogmod_ddm() covers the classic 4-parameter DDM plus the outlier component. The between-trial variability parameters would each need their own quadrature layer on top - and sigmadrift, which the density handles with a closed-form correction, has no such form here - so they are not arguments at all: passing one is an error rather than a silently wrong number. Integrate dcogmod_ddm() over q if you need them.

With a response, pcogmod_ddm() returns the defective CDF P(RT <= q, choice = response), which does not reach one: its limit is the probability of that boundary, given by pcogmod_ddm(Inf, response = k). The upper tail is then the matching defective survival P(RT > q, choice = response), so the two add to that response's own probability rather than to one. Marginally (response = NULL) they add to one as usual. pcogmod_rdm() follows the same convention.

References

• Ratcliff, R., & McKoon, G. (2008). The diffusion decision model: Theory and data for two-choice decision tasks. Neural Computation, 20(4), 873-922. doi:10.1162/neco.2008.12-06-420

• Navarro, D. J., & Fuss, I. G. (2009). Fast and accurate calculations for first-passage times in Wiener diffusion models. Journal of Mathematical Psychology, 53(4), 222-230. doi:10.1016/j.jmp.2009.02.003

See also

rcogmod_rdm(), rcogmod_lba2(), rcogmod_lnr()

Examples

# Simulate data, with 2% of trials from the outlier process
data <- rcogmod_ddm(1000,
  drift = 0.5, boundary = 1, bias = 0.5, ndt = 0.2, poutlier = 0.02
)
head(data)
#>          rt response
#> 1 0.6881871        0
#> 2 0.4426854        0
#> 3 0.3992289        1
#> 4 0.6083311        0
#> 5 0.3406405        0
#> 6 0.3247461        0

# Responses faster than ndt keep positive density, unlike the unmixed model
dcogmod_ddm(0.1, ndt = 0.2, response = 1, poutlier = 0.02)
#> [1] 0.03520653
dcogmod_ddm(0.1, ndt = 0.2, response = 1, 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_ddm_lpdf_expose()
lpdf(
  Y = 0.5, mu = 0.5, boundary = 1, bias = 0.5, sigmadrift = 0,
  sigmabias = 0, sigmandt = 0, ndt = 0.2, poutlier = 0.02, dec = 1
)
} # }