---
title: "How to Convert DMC's parameters"
output: rmarkdown::html_vignette
bibliography: REFERENCES.bib
vignette: >
  %\VignetteIndexEntry{How to Convert DMC's parameters}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, echo = FALSE, message=FALSE}
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]

```{r}
# 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: 

```{r}
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. 

```{r}
coef(dmc_s)
coef(dmc_ms)
```


## References


```{r, echo = F}
# "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)
```