How to Convert DMC's parameters"
In dRiftDM: Estimating (Time-Dependent) Drift Diffusion Models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(dRiftDM)
set.seed(1014)

When working with Drift-Diffusion Models (DDMs), like the Diffusion Model for Conflict Tasks [DMC, @Ulrichetal.2015], a recurring nuisance is that the size of the model's parameters depends on the time scale (seconds vs. milliseconds) and the diffusion constant. For example, in the original article by @Ulrichetal.2015, DMC was introduced on the time scale of milliseconds and with a diffusion constant of $\sigma = 4$. However, it is not uncommon for DDMs to be realized in the time scale of seconds and/or with a diffusion constant of $\sigma = 1$.

The motivation for this vignette is to show how to convert DMC's parameters for different time scales and diffusion constants, thereby answering a frequently asked question [see also Appendix B in @Koobetal.2023]. While the following functions and equations only hold for DMC, the problem addressed here is a general one and partially applies to other DDMs as well. Therefore, we hope that this vignette may even be interesting to researchers that are are not specifically concerned with DMC.

The Mathematical Equations

DMC has 8 parameters in its most general form:

muc the drift rate of the controlled process
b the (constant) boundary
non_dec the mean of the (normally-distributed) non-decision time
sd_non_dec the standard deviation of the (normally-distributed) non-decision time
tau the scale of the automatic process
a the shape of the automatic process
A the amplitude of the automatic process
alpha the shape parameter of the starting point distribution (a symmetric beta-distribution)

Scaling for the Diffusion Constant

We have to consider the diffusion constant whenever a parameter depends on the scale of the evidence space (i.e., if it is somehow related to the scale of the y-axis of the evidence accumulation process). This is the case for muc, b, and A of DMC. We can transform these parameters with:

$$ \theta_{new} = \theta_{\sigma_{old}} \cdot \frac{\sigma_{new}}{\sigma_{old}}\,. $$ Here, $\theta_{\sigma_{old}}$ is a parameter (e.g., muc, b, and A) in the scale of the previous diffusion constant, $\sigma_{old}$, and $\theta_{\sigma_{new}}$ is the transformed parameter after scaling it to the target diffusion constant, $\sigma_{new}$.

Example: To transform muc $= 0.5$ from a diffusion constant of $\sigma_{old} = 4$ to match a DDM with a diffusion constant of $\sigma_{new} = 1$, we compute $$0.5\cdot \frac{1}{4} = 0.125\,.$$

Scaling the Time Space

Time scaling is a bit more difficult, because we have to consider if and how a parameter is affected by the time and/or evidence-space scale.

non_dec, sd_non_dec, and tau depend primarily on the time space. The transformation from seconds to milliseconds, or vice versa, is straightforward:

$$\theta_{s} = \theta_{ms} / 1000 \quad \text{and} \quad \theta_{ms} = \theta_{s} \cdot 1000\,.$$

b and A depend primarily on the evidence space scale. Yet, we still transform them slightly.^[They depend indirectly on the time scale through the variance of the diffusion process.] For these parameters we can use

$$\theta_{s} = \theta_{ms} / \sqrt{1000} \quad \text{and} \quad \theta_{ms} = \theta_{s} \cdot \sqrt{1000}\,.$$

Finally, the drift rate muc describes the rate of evidence increase per time step. For it we have to use

$$\theta_s = \theta_{ms} \cdot \frac{1000}{\sqrt{1000}} \quad \text{and} \quad \theta_{ms} = \theta_{s} \cdot \frac{\sqrt{1000}}{1000}\,.$$

Example: To transform muc $= 0.5$ from milliseconds to seconds, we calculate $$0.5\cdot \frac{1000}{\sqrt{1000}} = 15.8\,.$$

An R Helper Function

It would be tedious to do all the transformations "by hand," so we present here a small helper function that implements these transformations. The function takes a named numeric vector of parameter values and returns the transformed values.^[We did not export it with dRiftDM because it does not apply to all DDMs, only to DMC]

# Input documentation:
# named_values: a named numeric vector
# sigma_old, sigma_new: the previous and target diffusion constants
# t_from_to: scaling of time (options: ms->s, s->ms, or none)
convert_prms <- function(named_values,
                         sigma_old = 4,
                         sigma_new = 1,
                         t_from_to = "ms->s") {
  # Some rough input checks
  stopifnot(is.numeric(named_values), is.character(names(named_values)))
  stopifnot(is.numeric(sigma_old), is.numeric(sigma_new))
  t_from_to <- match.arg(t_from_to, choices = c("ms->s", "s->ms", "none"))

  # Internal conversion function (takes a name and value pair, and transforms it)
  internal <- function(name, value) {
    name <- match.arg(
      name,
      choices = c("muc", "b", "non_dec", "sd_non_dec", "tau", "a", "A", "alpha")
    )

    # 1. scale for the diffusion constant
    if (name %in% c("muc", "b", "A")) {
      value <- value * (sigma_new / sigma_old)
    }


    # 2. scale for the time
    # determine the scaling per parameter (assuming s->ms)
    scale <- 1
    if (name %in% c("non_dec", "sd_non_dec", "tau")) scale <- 1000
    if (name %in% c("b", "A")) scale <- sqrt(1000)
    if (name %in% c("muc")) scale <- sqrt(1000) / 1000

    # adapt, depending on the t_from_to argument
    if (t_from_to == "ms->s") scale <- 1 / scale
    if (t_from_to == "none") scale <- 1

    value <- value * scale
  }

  # Apply the internal function to each element
  converted_values <- mapply(FUN = internal, names(named_values), named_values)

  return(converted_values)
}

Let's see if the conversion works by checking if model predictions are the same before and after scaling:

dmc_s <- dmc_dm()
prms_solve(dmc_s) # current parameter settings for sigma = 1 and seconds
quants_s <- calc_stats(dmc_s, "quantiles") # calculate predicted quantiles
head(quants_s) # show quantiles

# now the same with new diffusion constant of 4 and time scale in milliseconds
dmc_ms <- dmc_dm()
prms_solve(dmc_ms)["sigma"] <- 4 # new diffusion constant
prms_solve(dmc_ms)["t_max"] <- 3000 # 3000 ms is new max time
prms_solve(dmc_ms)["dt"] <- 1 # 1 ms steps
coef(dmc_ms) <- convert_prms(
  named_values = coef(dmc_ms), # the previous parameters
  sigma_old = 1, # diffusion constants
  sigma_new = 4,
  t_from_to = "s->ms" # how shall the time be scaled?
)
quants_ms <- calc_stats(dmc_ms, "quantiles") # calculate predicted quantiles
head(quants_ms) # show quantiles

Note: The transformed values in the unit of milliseconds and with a diffusion constant of 4 are similar in size to those obtained by @Ulrichetal.2015.

coef(dmc_s)
coef(dmc_ms)

References

# "Unit test" the function

# TEST 1 -> converting twice should lead to the same result as previously
a_model <- dmc_dm(instr = "a ~!")
convert_1 <- convert_prms(
  coef(a_model),
  sigma_old = 1,
  sigma_new = 4,
  t_from_to = "s->ms"
)

convert_2 <- convert_prms(
  convert_1,
  sigma_old = 4,
  sigma_new = 1,
  t_from_to = "ms->s"
)

stopifnot(convert_2 == coef(a_model))


# TEST 2 -> quantiles from above should be very similar
stopifnot(
  abs(quants_ms$Quant_corr - quants_s$Quant_corr * 1000) <= 0.001
)


# TEST 3 -> expectation based on "hand" formula
DMCfun_def <- c(A = 20, tau = 30, muc = 0.5, b = 75, non_dec = 300, sd_non_dec = 30, a = 2, alpha = 3)
exp <- convert_prms(DMCfun_def, sigma_new = 0.1, sigma_old = 4, t_from_to = "ms->s")
stopifnot(abs(exp["A"] - 0.01581139) < 0.0001)
stopifnot(abs(exp["tau"] - 0.030) < 0.0001)
stopifnot(abs(exp["muc"] - 0.3952847) < 0.0001)
stopifnot(abs(exp["b"] - 0.05929271) < 0.0001)
stopifnot(abs(exp["non_dec"] - 0.300) < 0.0001)
stopifnot(abs(exp["sd_non_dec"] - 0.030) < 0.0001)
stopifnot(abs(exp["a"] - 2) < 0.0001)
stopifnot(abs(exp["alpha"] - 3) < 0.0001)


# TEST 4 -> expectation based on "hand" formula (no time scaling)
exp <- convert_prms(DMCfun_def, sigma_new = 0.1, sigma_old = 4, t_from_to = "none")
stopifnot(abs(exp["A"] - 0.5) < 0.0001)
stopifnot(abs(exp["tau"] - 30) < 0.0001)
stopifnot(abs(exp["muc"] - 0.0125) < 0.0001)
stopifnot(abs(exp["b"] - 1.875) < 0.0001)
stopifnot(abs(exp["non_dec"] - 300) < 0.0001)
stopifnot(abs(exp["sd_non_dec"] - 30) < 0.0001)
stopifnot(abs(exp["a"] - 2) < 0.0001)
stopifnot(abs(exp["alpha"] - 3) < 0.0001)