brms_example.Rmd

---
title: "5-Parameter DDM in `brms` (`Stan`)"
author: "Kendal Foster and Henrik Singmann"
date: '`r format(Sys.Date(), "%B %d, %Y")`'
bibliography: markdown_stuff/references.bib
output:
  rmarkdown::html_vignette:
    css: markdown_stuff/stile.css
    toc: false
    fig_width: 8
    fig_height: 6
vignette: >
  %\VignetteIndexEntry{5-Parameter DDM in `brms` (`Stan`)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
```{r echo=FALSE}
req_suggested_packages <- c("brms", "rtdists", "fddm", "RWiener", "bayestestR", "ggplot2")
pcheck <- lapply(req_suggested_packages, requireNamespace, 
                 quietly = TRUE)
if (any(!unlist(pcheck))) {
   message("Required package(s) for this vignette are not available/installed and code will not be executed.")
   knitr::opts_chunk$set(eval = FALSE)
}
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  error = TRUE,
  comment = "#>"
)
op <- options(width = 100, digits = 4)
```

<div id="TOC">
<ul>
  <li><a href="#intro">Introduction</a></li>
  <li><a href="#cusfam">Custom `brms` Family: `ddm`</a>
  <ul>
    <li><a href="#cusfam-5par">Custom 5-Parameter DDM Definition</a></li>
    <li><a href="#cusfam-4par">Custom 4-Parameter DDM Definition</a></li>
    <li><a href="#cusfam-priors">Priors for the Custom Families</a></li>
  </ul></li>
  <li><a href="#parrec">Parameter Recovery</a>
  <ul>
    <li><a href="#parrec-data">Generating Some Data for the Parameter Recovery</a></li>
    <li><a href="#parrec-samp">Estimation (i.e., Sampling)</a></li>
    <li><a href="#parrec-post">Posterior Predictive Analysis</a></li>
  </ul></li>
  <li><a href="#rewrld">Fitting Real-World Data</a>
  <ul>
    <li><a href="#rewrld-data">Load Real-World Data</a></li>
    <li><a href="#rewrld-samp">Estimation (i.e., Sampling)</a></li>
    <li><a href="#rewrld-post">Posterior Predictive Analysis</a></li>
  </ul></li>
  <li><a href="#references">References</a></li>
</ul>
</div>





This document illustrates how to use the 5-parameter Diffusion Decision Model (DDM) (@ratcliff1978theory) with `brms` (@brms1, @brms2) and `rstan` (@rstan).

This  implementation of the DDM has the following 5 parameters: $a \in (0, \infty)$ (threshold separation), $v \in (-\infty, \infty)$ (drift rate), $t_0 \in [0, \infty)$ (non-decision time/response time constant), $w \in (0, 1)$ (relative starting point), and $\eta \in (0, \infty)$ (inter-trial-variability of drift).
<br><br>





# Introduction {#intro}
<hr class="sec1">

The Ratcliff diffusion decision model (DDM) (@ratcliff1978theory; @ratcliff2008diffusion) is the most prominent evidence accumulation model for jointly modelling binary decision and associated response times. It assumes noisy information uptake; this noise is characterized by the Wiener process. The DDM is popular because (a) it utilizes a lot of information from the collected data, which in turn allows us to draw conclusions about the underlying cognitive processes that are assumed to underlie the empirical response time distribution (@voss2008fast) and (b) it provides a good fit to data across different empirical domains (@forstmann2016sequential).

However, this efficient use of information comes at a rather high computational cost. The probability density function (PDF) of the DDM is the density of first-passage times, that is, the time it takes for the Wiener process to first cross a fixed boundary. The computation of the PDF of the first-passage times is computationally expensive because it contains an infinite sum. Of course computers cannot actually evaluate an infinite sum, so several papers have developed approximation methods to this density function (@voss2008fast; @navarro2009fast; @blurton2012fast; @gondan2014even). The different approximation methods provide different truncation rules that guarantee the value of the approximated density function does not differ from its true value by more than a pre-specified error tolerance.

The analytic formulation of the diffusion model PDF originates from @feller1968introduction, who provides its seminal derivation. In this derivation from first principles (i.e., the Wiener process), there is a step that requires taking a limit. @feller1968introduction provides two different – but equivalent – limiting processes that yield two different – but equal – forms of the density function. Each of these forms contains an infinite sum, and they are known individually as the “large-time” density function and the “small-time” density function (@navarro2009fast). As their names suggest, the "large-time" density function is more efficient at handling large response times, and the "small-time" density function is more efficient at calculating small response times. Even though the two distinct forms of the DDM density function are mathematically equivalent, they can produce slightly different results when calculated numerically with a fixed and limited precision.

For the remainder of this section, we will refer to the 5-parameter DDM density function as $f(t | v, \eta, a, w)$. Note that we use $t := t - t_0$ to make the notation cleaner.



## Large-Time Density Function {#intro-large}

The 5-parameter "large-time" PDF is

\begin{equation} \label{eq:pdf-l}
f_\ell(t | v, \eta, a, w) = \frac{\pi}{a^2 \sqrt{1 + \eta^2 t}}
                            \exp{ \left( \frac{\eta^2 a^2 w^2 -2vaw -v^2 t}{2 (1 + \eta^2 t)} \right)}
														\sum_{j = 1}^{k_\ell} j \sin \left( j w \pi \right) \exp{ \left( -\frac{j^2 \pi^2 t}{2a^2} \right)}
\end{equation}

The first truncation rule for the infinite sum in this density function was provided by @navarro2009fast; it uses the first $k_\ell$ of the infinite sum. It's still the only one out there.



## Small-Time Density Function {#intro-small}

The 5-parameter "small-time" PDF is

\begin{equation} \label{eq:pdf-s}
f_s(t | v, \eta, a, w) = \frac{a}{\sqrt{2 \pi t^3 \left( 1 + \eta^2 t \right)}}
											 	 \exp{ \left( \frac{\eta^2 a^2 w^2 -2vaw -v^2 t}{2 (1 + \eta^2 t)} \right)}
											 	 \sum_{j = -k_s}^{k_s} (w + 2j) \exp{ \left( -\frac{a^2}{2t} \left( w + 2j \right)^2 \right)}
\end{equation}

@navarro2009fast also provided the first truncation rule for the infinite sum in this density function; the truncated sum runs from term $-k_s$ to term $k_s$, $2k + 1$ terms in total. @gondan2014even provided another truncation rule for this infinite sum that follows the same premise, using the first $2k + 1$ terms. In their paper, @gondan2014even mention that the mathematical properties of the sum allow it to potentially be truncated sooner, but then the number of terms required for the truncated sum cannot be pre-calculated.



## Combining the Large-time and Small-time Density Functions {#intro-combined}

Along with their "large-time" and "small-time" truncation rules, @navarro2009fast also propose a mechanism to switch between the "large-time" and "small-time" PDFs so that the more efficient one is used. Given a pre-defined error tolerance, the mechanism calculates the number of terms required for the "large-time" PDF, $k_\ell$, and the "small-time" PDF, $k_s$; then it uses whichever PDF requires fewer terms. Even though @gondan2014even mentioned exploiting the mathematical properties of the infinite sum in the "small-time" PDF to speed up the approximation, they favored pre-calculating $k_s$ so that the same mechanism could directly compare the efficiency of the "large-time" and "small-time" PDFs.

Instead of precalculating $k_s$, we proposed using a heuristic to represent $k_s$; this way we could save computation time from not calculating $k_s$ and we could still see how efficient the "large-time" PDF could be. Our mechanism for choosing between the "large-time" and "small-time" PDF treats $k_s = 2$, so that if $k_\ell$ is calculated to be $0$ or $1$, then the "large-time" PDF is used; otherwise if $k_\ell > 1$, then the "small-time" PDF is used. We determined the heuristic value $1$ from empirical testing; the results of this testing can be seen in our arXiv paper (@foster2021approximation). Formulaically, it looks like

\begin{equation} \label{eq:pdf-c}
  f(t | v, \eta, a, w) = \begin{cases}
                           f_\ell(t | v, \eta, a, w) & \text{ if } k_\ell \le 1\\
                           f_s(t | v, \eta, a, w) & \text{ if } k_\ell > 1
                         \end{cases}
\end{equation}

For a pre-defined error tolerance, this method of combining the "large-time" and "small-time" is not only the fastest in terms of computation time, but it is also the most stable algorithm to use in optimization routines (e.g., fitting the model to data). You can see these results in our arXiv paper (@foster2021approximation). This method is the default in our `fddm` package and is the basis of the `ddm_lpdf()` function that we define in the `Stan` "functions" block (loaded into `rstan` via `brms` through the object `stan_funs` defined [below](#stan-func-def)).





# Custom `brms` Family: `ddm` {#cusfam}
<hr class="sec1">

As `rstan` is a great tool for performing Bayesian analysis, we thought that it should have access to the 5-parameter DDM. Moreover, the popular `R` package `brms` is a widely used interface to `rstan` that natively handles many popular distributions, but it too is missing the 5-parameter DDM. This section will detail how to define the 5-parameter DDM for use in `brms` (and also `rstan`).

Let's begin by loading the packages that we'll use.

```{r load-pkgs}
library("brms")
library("rtdists")
library("fddm")
library("RWiener")
library("bayestestR")
library("ggplot2")
library("cowplot")
```

The package `brms` contains many distributions, including the Wiener distribution, which is equivalent to the DDM; however, `brms` only natively supports the 4-parameter Wiener distribution (it excludes the inter-trial variability in the drift rate, $sv$). To use the 5-parameter DDM with `brms`, we must define a "custom family" that contains all of the necessary information about the distribution to pass to `Stan`. We will compare the custom 5-parameter variant with the 4-parameter variant that is native in `brms`. In addition, we will also define a custom family for the 4-parameter variant using the same algorithm from the 5-parameter variant; we will compare this variant to the native 4-parameter variant as well as our other custom family for the 5-parameter variant. After defining the custom family for the 5-parameter variant, we will define the custom family for the 4-parameter variant (the process will be nearly identical).



## Custom 5-Parameter DDM Definition {#cusfam-5par}

First, we define the custom family "ddm" for the 5-parameter variant.

```{r brms-family}
ddm <- custom_family(
  "ddm", dpars = c("mu", "a", "ndt", "w", "sv"),
  links = c("identity", "log", "log", "logit", "log"), 
  lb = c(NA, 0, 0, 0, 0), ub = c(NA, NA, NA, 1, NA),
  type = "real", vars = "resp[n]"
)
```

Note that we have changed the parameter $v$ to $mu$ because `brms` demands that there be a parameter called "mu" in the family. We also changed the parameter `t0` to `ndt` because `brms` does not allow a parameter to end in a number. In addition, we changed the parameter $\eta$ to `sv` because that is how it is written in our software package. Each parameter in the family must have a mapping, defined in the `links` argument (in the order defined in `dpars` in the line above). `lb` and `ub` are vectors of the lower and upper bounds of the parameters (in the order defined in `dpars` two lines above), respectively. The argument `type` is the type of data that the distribution handles; in our case, the DDM handles `real` numbers (or `double`s in other programming languages) as opposed to integers (`int`s). The last argument, `vars`, allows us to include other data that are relative to the distribution; as the data are actually tuples of responses and associated response times, this arguments allows us to input a vector of responses as the `resp` variable.

We must also provide the probability density function for the underlying `Stan` sampler to use. This function is written in the `Stan` language, and we will input it to `brms` later. Note that this function accepts only single real values for each of the inputs (aside from `resp`, which is of the `int` type). For a version of this density function that provides support for vectors of response times (`rt`) and responses (`resp`), see the file `extra_stan_stuff/brms_ddm_stan_funcs.stan`.

<a id="stan-func-def"></a>

```{r brms-stanfuncs}
stan_funs <- "
real ddm_lpdf(real rt, real v, real a, real ndt, real w, real sv, int resp)
{
  // initialize output
  real logp = 0.0;
  
  // check parameter values
  if (a <= 0 || is_inf(a) || is_nan(a)) return negative_infinity();
  if (is_inf(v) || is_nan(v)) return negative_infinity();
  if (ndt < 0 || is_inf(ndt) || is_nan(ndt)) return negative_infinity();
  if (w <= 0 || w >= 1 || is_nan(w)) return negative_infinity();
  if (sv < 0 || is_inf(sv) || is_nan(sv)) return negative_infinity();
  
  // adjust parameters
  {
  real t = rt - ndt;
  real a_i = a;
  real sv_i = sv;
  real v_i = (resp == 1) ? v : -v;
  real w_i = (resp == 1) ? w : 1 - w;
  real taa = t / (a_i * a_i);
  
  // define variables for use later (have to do it early because of Stan)
  // constants
  real log_2_pi_2 = 0.5 * log(2 * pi());
  real sv_thresh = 0.0;
  // used in both large-time and small-time
  real mult = 0.0;
  real exp_err = 0.0;
  real gamma = -0.5 * pi() * pi() * taa;
  real summ = 0.0;
  int max_int = 1000;
  // used only in large-time
  int kl = 0;
  real k_dbl = 0.0;
  real bc = 1 / (pi() * sqrt(taa));
  // used only in small-time
  real minterms = 0.5 * sqrt(taa) - 0.5 * w_i; // min number of terms, truncates toward 0
  real term = 0.0;
  real rj = 0.0;
  int js = 0;
  
  // check large-time number of terms
  if (sv_i < sv_thresh) {
    mult = - v_i * a_i * w_i - 0.5 * v_i*v_i * t - 2 * log(a_i);
  } else {
    mult = (sv_i*sv_i * a_i*a_i * w_i*w_i - 2 * v_i * a_i * w_i - v_i*v_i * t)
    / (2 + 2 * sv_i*sv_i * t) - 0.5 * log(1 + sv_i*sv_i * t) - 2 * log(a_i);
  }
  exp_err = 0.000001 * exp(-mult);
  if (bc > max_int) { // boundary condition
    k_dbl = max_int;
  } else {
    if (exp_err * pi() * taa < 1) { // error threshold is low enough
      k_dbl = sqrt(-2 * log(pi() * taa * exp_err) / (pi() * pi() * taa));
      if (k_dbl > max_int) {
        k_dbl = max_int;
      } else {
        k_dbl = fmax(k_dbl, bc);
      }
    } else { // threshold not low enough, so set to boundary condition
      k_dbl = bc;
    }
  }
  while (kl < k_dbl) kl += 1;
  
  // compare large-time and small-time
  if (kl <= 1) { // do large-time
    for(jl in 1:kl) {
      summ += jl * sin(jl * w_i * pi()) * exp(gamma * jl * jl);
    }
    logp += log(pi());
  } else { // do small-time
    if (sv_i < sv_thresh) {
      mult = log(a_i) - log_2_pi_2 - 1.5 * log(t) - v_i * a_i * w_i
      - 0.5 * v_i*v_i * t;
    } else {
      mult = log(a_i) - 1.5 * log(t) - log_2_pi_2
      - 0.5 * log(1 + sv_i*sv_i * t)
      + (sv_i*sv_i * a_i*a_i * w_i*w_i - 2 * v_i * a_i * w_i
      - v_i*v_i * t) / (2 + 2 * sv_i*sv_i * t);
    }
    exp_err = 0.000001 * exp(-mult);
    gamma = -1 / (2 * taa);
    summ = w_i * exp(gamma * w_i*w_i); // initialize with j=0 term
    while (js < minterms) {
      js += 1;
      rj = 2 * js - w_i;
      summ -= rj * exp(gamma * rj*rj);
      rj = 2 * js + w_i;
      summ += rj * exp(gamma * rj*rj);
    }
    js += 1;
    rj = 2 * js - w_i;
    term = rj * exp(gamma * rj*rj);
    summ -= term;
    while (term > exp_err) {
      rj = 2 * js + w_i;
      term = rj * exp(gamma * rj*rj);
      summ += term;
      if (term <= exp_err) break;
      js += 1;
      rj = 2 * js - w_i;
      term = rj * exp(gamma * rj*rj);
      summ -= term;
    }
    
    // check summ and chuck everything in logp
    if (summ >= 0) {  // if result is negative, don't add to logp
      logp += mult + log(summ);
    }
  }
  }
  
  return logp;
}
"
```

We'll save this in the object `stanvars`, to which we'll add later.

```{r brms-stanvars}
stanvars <- stanvar(scode = stan_funs, block = "functions")
```

The last bit of prep work that we need to do in order to run `Stan` sampling is to define the priors on the model parameters. These priors were chosen to be mildly informative, as we have an idea of the reasonable range of the parameters. These priors may seem strange in that they support the negative real numbers when the parameters cannot be negative. In the definition of the custom family, we included "link" functions which will transform the parameters and their priors; more details about this can be found [below](#cusfam-priors). Note that the prior for the non-decision time (`ndt`) technically extends to positive infinity, but in practice the non-decision time is bounded above by the smallest response time in the dataset.

```{r brms-priors}
priors <- c(
  set_prior("normal(0, 2.5)", class = "Intercept"), # for mu (which is v)
  set_prior("logistic(0.75, 0.5)", class = "a"),
  set_prior("logistic(-1, 0.5)", class = "ndt"),
  set_prior("logistic(0, 0.67)", class = "w"),
  set_prior("logistic(-0.5, 0.5)", class = "sv")
)
```

<a id="brms-pp"></a>

Although we could begin with `Stan` sampling now, we'll instead define two functions that will later help us with the posterior analysis.

First, we define the log-likelihood function (to be used in `loo`). Note the particular naming scheme used for the function name and the specific arguments the function expects. This log-likelihood function requires the calculation of the 5-parameter DDM probability density function; for this calculation, we use the package `fddm` as it contains the most recent implementation of the 5-parameter DDM density function.

```{r brms-loglik}
log_lik_ddm <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  a <- brms::get_dpar(prep, "a", i = i)
  w <- brms::get_dpar(prep, "w", i = i)
  ndt <- brms::get_dpar(prep, "ndt", i = i)
  sv <- brms::get_dpar(prep, "sv", i = i)
  y <- prep$data$Y[i]
  resp <- prep$data$resp[i]
  fddm::dfddm(y, resp, a, mu, ndt, w, sv, log = TRUE)
}
```

Second, we define the posterior predictive function (to be used in `pp_check`). Again, note the particular naming scheme used for the function name and the specific arguments the function expects. We use the package `rtdists` to generate random responses and associated response times in this function.

```{r brms-postpred}
posterior_predict_ddm <- function(i, prep, negative_rt = FALSE, ...) {
  if (!require("rtdists")) {
    stop("Package 'rtdists' must be installed for the function 'posterior_predict_ddm'.")
  }
  
  v <- as.numeric(brms::get_dpar(prep, "mu", i = i))
  a <- as.numeric(brms::get_dpar(prep, "a", i = i))
  w <- as.numeric(brms::get_dpar(prep, "w", i = i))
  t0 <- as.numeric(brms::get_dpar(prep, "ndt", i = i))
  sv <- as.numeric(brms::get_dpar(prep, "sv", i = i))
  n <- max(lengths(list(v, a, w, t0, sv)), 1) # I'm pretty sure it's always 1
  
  out <- rtdists::rdiffusion(
    n = n,
    a = a,
    v = v,
    t0 = t0,
    z = w * a,
    sv = sv
  )
  colnames(out) <- c("q", "resp")

  if (negative_rt) { # code "lower" responses as negative RTs
    out[["q"]] <- out[["q"]] * ifelse(out[["resp"]] == "upper", 1, -1)
  }
  
  return(out[["q"]])
}
```



## Custom 4-Parameter DDM Definition {#cusfam-4par}

Next we will define the custom family "ddm4" for the 4-parameter variant that uses the PDF algorithm from the `fddm` package. This custom family will act as a direct check with the native Wiener family, as well as providing a reference for the time it takes for `Stan` to sample. Defining this 4-parameter custom family is nearly identical to what we just did for the "ddm" custom family above, but we add the suffix "4" to distinguish it.

```{r brms-ddm4}
ddm4 <- custom_family(
  "ddm4", dpars = c("mu", "a", "ndt", "w"),
  links = c("identity", "log", "log", "logit"), 
  lb = c(NA, 0, 0, 0), ub = c(NA, NA, NA, 1),
  type = "real", vars = "resp[n]"
)

stan_funs4 <- "
real ddm4_lpdf(real rt, real v, real a, real ndt, real w, int resp)
{
  // initialize output
  real logp = 0.0;
  
  // check parameter values
  if (a <= 0 || is_inf(a) || is_nan(a)) return negative_infinity();
  if (is_inf(v) || is_nan(v)) return negative_infinity();
  if (ndt < 0 || is_inf(ndt) || is_nan(ndt)) return negative_infinity();
  if (w <= 0 || w >= 1 || is_nan(w)) return negative_infinity();

  // adjust parameters
  {
  real t = rt - ndt;
  real a_i = a;
  real v_i = (resp == 1) ? v : -v;
  real w_i = (resp == 1) ? w : 1 - w;
  real taa = t / (a_i * a_i);
  
  // define variables for use later (have to do it early because of Stan)
  // constants
  real log_2_pi_2 = 0.5 * log(2 * pi());
  // used in both large-time and small-time
  real mult = 0.0;
  real exp_err = 0.0;
  real gamma = -0.5 * pi() * pi() * taa;
  real summ = 0.0;
  int max_int = 1000;
  // used only in large-time
  int kl = 0;
  real k_dbl = 0.0;
  real bc = 1 / (pi() * sqrt(taa));
  // used only in small-time
  real minterms = 0.5 * sqrt(taa) - 0.5 * w_i; // min number of terms, truncates toward 0
  real term = 0.0;
  real rj = 0.0;
  int js = 0;
  
  // check large-time number of terms
  mult = - v_i * a_i * w_i - 0.5 * v_i*v_i * t - 2 * log(a_i);
  exp_err = 0.000001 * exp(-mult);
  if (bc > max_int) { // boundary condition
    k_dbl = max_int;
  } else {
    if (exp_err * pi() * taa < 1) { // error threshold is low enough
      k_dbl = sqrt(-2 * log(pi() * taa * exp_err) / (pi() * pi() * taa));
      if (k_dbl > max_int) {
        k_dbl = max_int;
      } else {
        k_dbl = fmax(k_dbl, bc);
      }
    } else { // threshold not low enough, so set to boundary condition
      k_dbl = bc;
    }
  }
  while (kl < k_dbl) kl += 1;
  
  // compare large-time and small-time
  if (kl <= 1) { // do large-time
    for(jl in 1:kl) {
      summ += jl * sin(jl * w_i * pi()) * exp(gamma * jl * jl);
    }
    logp += log(pi());
  } else { // do small-time
    mult = log(a_i) - log_2_pi_2 - 1.5 * log(t) - v_i * a_i * w_i
           - 0.5 * v_i*v_i * t;
    exp_err = 0.000001 * exp(-mult);
    gamma = -1 / (2 * taa);
    summ = w_i * exp(gamma * w_i*w_i); // initialize with j=0 term
    while (js < minterms) {
      js += 1;
      rj = 2 * js - w_i;
      summ -= rj * exp(gamma * rj*rj);
      rj = 2 * js + w_i;
      summ += rj * exp(gamma * rj*rj);
    }
    js += 1;
    rj = 2 * js - w_i;
    term = rj * exp(gamma * rj*rj);
    summ -= term;
    while (term > exp_err) {
      rj = 2 * js + w_i;
      term = rj * exp(gamma * rj*rj);
      summ += term;
      if (term <= exp_err) break;
      js += 1;
      rj = 2 * js - w_i;
      term = rj * exp(gamma * rj*rj);
      summ -= term;
    }
    
    // check summ and chuck everything in logp
    if (summ >= 0) {  // if result is negative, don't add to logp
      logp += mult + log(summ);
    }
  }
  }
  
  return logp;
}
"

stanvars4 <- stanvar(scode = stan_funs4, block = "functions")

priors4 <- c(
  set_prior("normal(0, 2.5)", class = "Intercept"), # for mu (which is v)
  set_prior("gamma(3, 1.5)", class = "a"),
  set_prior("gamma(2, 4)", class = "ndt"),
  set_prior("beta(2, 2)", class = "w")
)

log_lik_ddm4 <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  a <- brms::get_dpar(prep, "a", i = i)
  w <- brms::get_dpar(prep, "w", i = i)
  ndt <- brms::get_dpar(prep, "ndt", i = i)
  y <- prep$data$Y[i]
  resp <- prep$data$resp[i]
  fddm::dfddm(y, resp, a, mu, ndt, w, 0, log = TRUE)
}

posterior_predict_ddm4 <- function(i, prep, negative_rt = FALSE, ...) {
  if (!require("rtdists")) {
    stop("Package 'rtdists' must be installed for the function 'posterior_predict_ddm'.")
  }
  
  v <- as.numeric(brms::get_dpar(prep, "mu", i = i))
  a <- as.numeric(brms::get_dpar(prep, "a", i = i))
  w <- as.numeric(brms::get_dpar(prep, "w", i = i))
  t0 <- as.numeric(brms::get_dpar(prep, "ndt", i = i))
  n <- max(lengths(list(v, a, w, t0)), 1) # I'm pretty sure it's always 1
  
  out <- rtdists::rdiffusion(
    n = n,
    a = a,
    v = v,
    t0 = t0,
    z = w * a
  )
  colnames(out) <- c("q", "resp")

  if (negative_rt) { # code "lower" responses as negative RTs
    out[["q"]] <- out[["q"]] * ifelse(out[["resp"]] == "upper", 1, -1)
  }
  
  return(out[["q"]])
}
```



## Priors for the Custom Families {#cusfam-priors}

We defined the priors for the 5-parameter and 4-parameter custom families in the previous sections, but they might appear strange at first glance since most of the parameters are transformed according to their "links" in the call to `brms::custom_family()`. In particular, the parameters `a`, `ndt`, and `sv` (where applicable) use the "log" link function because their valid range is the positive real numbers; the parameter `w` uses the "logit" link function because it must be in the range $(0, 1)$. The parameter `mu` (`v`) uses the identity link function, so it is not transformed.

Since the parameters `a`, `ndt`, `w`, and `sv` are transformed in the custom family definition, their priors are also transformed, and we must consider these transformations when choosing the priors. By inspecting the output of the `brms::make_stancode()` function, we see that the inverse transformation is applied to the priors. For example, consider the prior that we want to place on the parameter `a`; call it $p(a)$. Since `a` uses the "log" link function, $p(a)$ will get transformed, and the sampler will apply the inverse link function and use $exp(p(a))$ as the prior. Working with the exponentiated version of a distribution can be awkward and may require lots of fiddling to adjust near the prior that we want. However, we can make this guess-and-check process a bit easier if we cleverly choose a prior distribution that is both flexible and has a well-known log-version.

This clever choice of prior distribution makes our life easier in two ways. First, it helps with knowing the summary statistics of the distribution that will actually be used, such as the mean, median, mode, and variance of the log-distribution. Understanding how the distribution parameters affect these summary statistics remove a lot of "guessing" from trying to shape the prior. Second, it allows for easy conversion from the log-distribution parameters to the regular distribution parameters; this is necessary as we need to input the non-log version of the distribution into `brms`.

The priors that we chose for the parameters `a`, `ndt`, and `sv` employ this technique by exploiting the relationship between the logistic and log-logistic probability distributions. The logistic distribution is defined such that its CDF is the logistic function; consequently, its PDF is very similar to that of the normal distribution, but with heavier tails. Upon exponentiation, the logistic distribution becomes the log-logistic distribution (i.e., the log of a log-logistic random variable follows the logistic distribution), and we can use the summary statistics to guide our fiddling with the shape of the prior (e.g., the median is just the scale parameter). Once we are satisfied with the shape of the log-logistic prior, we must convert the hyperparameters (of the log-logistic distribution) to the logistic distribution; conveniently, there exists an alternate parameterization for the log-logistic distribution that relates the parameters back to the logistic distribution.

The prior for `w` gets transformed by the logistic function, so we just fiddle with plotting a prior under its inverse transform until we are satisfied with the shape.

Below are the (post-transformation) priors:

```{r brms-priors-plot}
n <- 1e5

p_a <- ggplot(data = data.frame(logis = exp(rlogis(n, location = 0.75, scale = 0.5)))) +
  geom_histogram(aes(x = logis), binwidth = 0.25, boundary = 0, alpha = 0.7) +
  coord_cartesian(xlim = c(0, 10)) +
  labs(x = "a, threshold separation", y = "Count") +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

p_ndt <- ggplot(data = data.frame(logis = exp(rlogis(n, location = -1, scale = 0.5)))) +
  geom_histogram(aes(x = logis), binwidth = 0.05, boundary = 0, alpha = 0.7) +
  coord_cartesian(xlim = c(0, 1)) +
  labs(x = "ndt, non-decision time", y = "Count") +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

p_w <- ggplot(data = data.frame(logis = plogis(rlogis(n, location = 0, scale = 2/3)))) +
  geom_histogram(aes(x = logis), binwidth = 0.05, boundary = 0, alpha = 0.7) +
  coord_cartesian(xlim = c(0, 1)) +
  labs(x = "w, initial bias", y = "Count") +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

p_sv <- ggplot(data = data.frame(logis = exp(rlogis(n, location = -0.5, scale = 0.5)))) +
  geom_histogram(aes(x = logis), binwidth = 0.1, boundary = 0, alpha = 0.7) +
  coord_cartesian(xlim = c(0, 2.5)) +
  labs(x = "sv, variability in the drift rate", y = "Count") +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

cowplot::plot_grid(p_a, p_ndt, p_sv, p_w, ncol = 2)
```





# Parameter Recovery {#parrec}
<hr class="sec1">

Now that we've got our custom families ready in `brms`, we'll do a basic parameter recovery for each of the custom families and compare it with the native 4-parameter family to show that everything works as expected.



## Generating Some Data for the Parameter Recovery {#parrec-data}

Using the package `rtdists` (@rtdists), we will generate some response times using the given set of parameters. For reproducibility, we will set the seed.

```{r pr-gendata, eval=FALSE, include=TRUE}
a <- 1
v <- -2
t0 <- 0.3
w <- 0.5 # rtdists uses z = w*a, defaults to a/2
sv <- 0.8
N <- 1000
set.seed(1234567)
par_rec_data <- rtdists::rdiffusion(n = N, a = a, v = v, t0 = t0, sv = sv)
```

```{r pr-savegendata, eval=FALSE, include=FALSE}
saveRDS(par_rec_data,
        file = file.path(getwd(), "fits", "parameter_recovery_data.RDS"))
```

```{r pr-readgendata, include=FALSE}
par_rec_data <- readRDS(file.path(getwd(), "fits",
                                  "parameter_recovery_data.RDS"))
```

Next, we need to update our `stanvars` object to incorporate the data we just generated

```{r pr-stanvars}
stanvars_pr <- stanvars +
  stanvar(x = as.integer(par_rec_data[["response"]]), 
          name = "resp", scode = "int resp[N];")
```

We must also update `stanvars` for our 4-parameter custom family.

```{r pr-stanvars4}
stanvars_pr4 <- stanvars4 +
  stanvar(x = as.integer(par_rec_data[["response"]]),
          name = "resp", scode = "int resp[N];")
```

If we desire, we may check the `Stan` code that `brms` will generate to do the sampling. This step is optional, but it is a good opportunity to catch errors that we made in the code.

```{r pr-check-stancode, eval=FALSE, include=TRUE}
make_stancode(rt ~ 1,
              data = par_rec_data,
              family = ddm,
              stanvars = stanvars_pr,
              prior = priors
              )
```



## Estimation (i.e., Sampling) {#parrec-samp}

Now we let `brms` call `Stan` to run its sampling algorithm and generate estimates for the model parameters. Estimation usually takes on the scale of 90 seconds per chain with 2000 iterations and the default control options (i.e., `max_treedepth = 10` and `adapt_delta = 0.8`). Because this process can take a little while, we can also load the fit from file.

```{r pr-sampling, eval=FALSE, include=TRUE}
fit_parrec_5par <- brm(rt ~ 1,
                       family = ddm,
                       prior = priors,
                       stanvars = stanvars_pr,
                       data = par_rec_data,
                       chains = 1,
                       iter = 2000,
                       warmup = 1000,
                       thin = 1,
                       cores = getOption("mc.cores", 1),
                       control = list(
                         max_treedepth = 10,
                         adapt_delta = 0.8
                         )
                       )
```

```{r pr-check-if-already-fit, include=FALSE}
if (file.exists(file.path(getwd(), "fits", "fit_parrec_5par.RDS"))) {
  fit_parrec_5par <- readRDS(file = file.path(getwd(), "fits",
                                              "fit_parrec_5par.RDS"))
} else { # run fit
  fit_parrec_5par <- brm(rt ~ 1,
                         family = ddm,
                         prior = priors,
                         stanvars = stanvars_pr,
                         data = par_rec_data,
                         chains = 1,
                         iter = 2000,
                         warmup = 1000,
                         thin = 1,
                         cores = getOption("mc.cores", 1),
                         control = list(
                           max_treedepth = 10,
                           adapt_delta = 0.8
                           )
                         )
  saveRDS(fit_parrec_5par,
        file = file.path(getwd(),"fits", "fit_parrec_5par.RDS"))
}
```

We will also save the fit in case we need it for future use.

```{r pr-save-fit, eval=FALSE, include=TRUE}
saveRDS(fit_parrec_5par,
        file = file.path(getwd(), "fits", "fit_parrec_5par.RDS"))
```

Similarly for the 4-parameter custom family:

```{r pr-sampling4, eval=FALSE, include=TRUE}
fit_parrec_4par <- brm(rt ~ 1,
                       family = ddm4,
                       prior = priors4,
                       stanvars = stanvars_pr4,
                       data = par_rec_data,
                       chains = 1,
                       iter = 2000,
                       warmup = 1000,
                       thin = 1,
                       cores = getOption("mc.cores", 1),
                       control = list(
                         max_treedepth = 10,
                         adapt_delta = 0.8
                         )
                       )
```

```{r pr-check-if-already-fit4, include=FALSE}
if (file.exists(file.path(getwd(), "fits", "fit_parrec_4par.RDS"))) {
  fit_parrec_4par <- readRDS(file = file.path(getwd(), "fits",
                                              "fit_parrec_4par.RDS"))
} else { # run fit
  fit_parrec_4par <- brm(rt ~ 1,
                         family = ddm4,
                         prior = priors4,
                         stanvars = stanvars_pr4,
                         data = par_rec_data,
                         chains = 1,
                         iter = 2000,
                         warmup = 1000,
                         thin = 1,
                         cores = getOption("mc.cores", 1),
                         control = list(
                           max_treedepth = 10,
                           adapt_delta = 0.8
                           )
                         )
  saveRDS(fit_parrec_4par,
        file = file.path(getwd(), "fits", "fit_parrec_4par.RDS"))
}
```

```{r pr-save-fit4, eval=FALSE, include=TRUE}
saveRDS(fit_parrec_4par,
        file = file.path(getwd(), "fits", "fit_parrec_4par.RDS"))
```

And again, we estimate using the native 4-parameter Wiener distribution. Note that sampling a native distribution is more straightforward than a custom family. We only need to define the priors on the model parameters and relabel the generated data.

```{r pr-prepw}
priorsw <- c(
  set_prior("normal(0, 2.5)", class = "Intercept"), # for mu (which is v)
  set_prior("logistic(0.75, 0.5)", class = "bs"),
  set_prior("logistic(-1, 0.5)", class = "ndt"),
  set_prior("logistic(0, 0.65)", class = "bias")
)

par_rec_dataw <- par_rec_data
par_rec_dataw[["response2"]] <- as.integer(par_rec_dataw[["response"]]) - 1
```

Then we can proceed with sampling.

```{r pr-samplingw, eval=FALSE, include=TRUE}
fit_parrec_wiener <- brm(rt | dec(response2) ~ 1,
                         family = wiener(link_bs = "log",
                                         link_ndt = "log",
                                         link_bias = "logit"),
                         prior = priorsw,
                         data = par_rec_dataw,
                         chains = 1,
                         iter = 2000,
                         warmup = 1000,
                         thin = 1,
                         cores = getOption("mc.cores", 1),
                         control = list(
                           max_treedepth = 10,
                           adapt_delta = 0.8
                           )
                         )
```

```{r pr-check-if-already-fitw, include=FALSE}
if (file.exists(file.path(getwd(), "fits", "fit_parrec_wiener.RDS"))) {
  fit_parrec_wiener <- readRDS(file = file.path(getwd(), "fits",
                                                "fit_parrec_wiener.RDS"))
} else { # run fit
  fit_parrec_wiener <- brm(rt | dec(response2) ~ 1,
                           family = wiener(link_bs = "log",
                                           link_ndt = "log",
                                           link_bias = "logit"),
                           prior = priorsw,
                           data = par_rec_dataw,
                           chains = 1,
                           iter = 2000,
                           warmup = 1000,
                           thin = 1,
                           cores = getOption("mc.cores", 1),
                           control = list(
                             max_treedepth = 10,
                             adapt_delta = 0.8
                             )
                           )
  saveRDS(fit_parrec_wiener,
        file = file.path(getwd(), "fits", "fit_parrec_wiener.RDS"))
}
```

```{r pr-save-fitw, eval=FALSE, include=TRUE}
saveRDS(fit_parrec_wiener,
        file = file.path(getwd(), "fits", "fit_parrec_wiener.RDS"))
```

From the `Stan` output, we see that fitting the 5-parameter custom family DDM took 100 seconds. Contrastingly, it took about 40 seconds per chain to fit either the 4-parameter custom family DDM or the 4-parameter Wiener family native to `brms`.

Now we'll take a quick look at the summaries of the sampling.

```{r pr-print4}
print(fit_parrec_4par)
print(fit_parrec_wiener)
```

The 4-parameter custom family and native 4-parameter Wiener family yield extremely similar estimates and successfully recover three of the four parameters; the threshold separation ($a$), non-decision time ($t0$), and initial bias ($w$) are all recovered successfully. The estimates for the drift rate ($v$) are somewhat close to the underlying value in both families, but remain biased ($-1.82$ vs $-2.00$).

```{r pr-print5}
print(fit_parrec_5par)
```

The 5-parameter custom family successfully recovers four of the five parameters successfully, with the only exception being the variability in the drift rate ($sv$). The estimates for $sv$ are somewhat close to the underlying value, but there is still moderate uncertainty in this estimate. Importantly, however, the 5-parameter custom family is able to successfully recover the drift rate ($v$) without bias.

We can see the differences in the estimates for the drift rate by plotting histograms of the posterior draws from each family.

```{r pr-plot-v}
n <- 1000 # number of posterior draws (excluding warmup)
models <- c("Wiener", "4-par DDM", "5-par DDM")
colors <- c("#b34d4d", "#4da7b3", "#5cc639")

df_v <- data.frame(
  model = rep(models, each = n),
  v = c(
    rstan::extract(fit_parrec_wiener[["fit"]],
                   pars = "b_Intercept")[["b_Intercept"]],
    rstan::extract(fit_parrec_4par[["fit"]],
                   pars = "b_Intercept")[["b_Intercept"]],
    rstan::extract(fit_parrec_5par[["fit"]],
                   pars = "b_Intercept")[["b_Intercept"]]
  )
)

ggplot(data = df_v,
       aes(x = v,
           color = factor(model, levels = models),
           fill = factor(model, levels = models))) +
  geom_density(alpha = 0.5) +
  scale_color_manual(values = colors, guide = "none") +
  scale_fill_manual(values = colors,
                    name = "Model",
                    labels = c("4-par Wiener", "4-par DDM", "5-par DDM")) +
  labs(x = "v, drift rate", y = "Posterior Density") +
  guides(fill = guide_legend(override.aes = list(color = colors))) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        legend.position = c(0.2, 0.82),
        legend.box = "vertical",
        legend.direction = "vertical",
        legend.background = element_rect(fill = "white"),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 13))
```

The main difference between the 4-parameter and 5-parameter DDM families (aside from the fitting times) is that the 5-parameter DDM captures variability in the drift rate, whereas the 4-parameter DDM attempts to estimate the variable drift rate by adding a bias. There is no difference in the estimates of the other parameters (threshold separation, non-decision time, initial bias) between the 4-parameter and 5-parameter DDM.



## Posterior Predictive Analysis {#parrec-post}

We will use the two functions that [we defined previously](#brms-pp) (the log-likelihood function and the posterior predictive function) to aid `brms` with further analysis of our custom family.

The `loo` function can compare two models; in our case, we will use it to perform three pairwise comparisons between the three models that we are testing. First, we'll compare the two custom `brms` families:

```{r pr-loo-custom}
loo(fit_parrec_5par, fit_parrec_4par)
```

Second, we'll compare the two 4-parameter variants:

```{r pr-loo-4par}
loo(fit_parrec_4par, fit_parrec_wiener)
```

Third, we'll compare the custom 5-parameter variant to the native 4-parameter Wiener distribution:

```{r pr-loo-5cuswiener}
loo(fit_parrec_5par, fit_parrec_wiener)
```

Based on the output of the `loo` function, we can see that the 5-parameter model has an ever so slightly higher estimate for ELPD than the 4-parameter models. However, the standard error on the difference between the ELPD estimates is much larger than the actual difference between the ELPD estimates. Therefore, `loo` suggests that these models perform equally well for the generated data. The main difference between the 5-parameter and 4-parameter models is potential biasing of the drift rate $v$ if its variability is not included in the model (parameter $sv$).

Last, we can use the posterior predictive function that [we defined earlier](#brms-pp) to run posterior-predictive checking for each of the three models

```{r pr-pp-5par}
pp_check(fit_parrec_5par)
```

```{r pr-pp-4par}
pp_check(fit_parrec_4par)
```

```{r pr-pp-4parw}
pp_check(fit_parrec_wiener)
```




# Fitting Real-World Data {#rewrld}
<hr class="sec1">

Now we'll try fitting some actual data, because this is how the software is intended to be used.



## Load Real-World Data {#rewrld-data}

Load the speed-accuracy data from the `rtdists` package and only consider the frequency "high" for words and non-words. Also update the stanvars to include this dataset.

```{r rw-data}
data(speed_acc, package = "rtdists")
speed_acc <- droplevels(speed_acc[!speed_acc$censor,]) # remove extreme RTs
speed_acc <- droplevels(speed_acc[ speed_acc$frequency %in% 
                                     c("high", "nw_high"),])

stanvars_rw <- stanvars +
  stanvar(x = as.integer(speed_acc[["response"]]), 
          name = "resp", scode = "int resp[N];")
```



## Estimation (i.e., Sampling) {#rewrld-samp}

We need to define the modelling formula and set the priors. We also set the contrast to orthonormal (which is why the `bayestestR` package is loaded) and define a function to randomly generate initial values for the sampling.

```{r rw-prep}
options(contrasts = c('contr.orthonorm', 'contr.poly'))

formula <- brmsformula(rt ~ condition*frequency,
                       a ~ 0 + Intercept + condition,
                       ndt ~ 0 + Intercept + condition,
                       w ~ 0 + Intercept + condition,
                       sv ~ 1,
                       center = FALSE
                       )

priors <- c(
  set_prior("student_t(3, 0, 2.5)", class = "b", coef = "Intercept"), # for mu (which is v)
  set_prior("logistic(0.75, 0.5)", coef = "Intercept", dpar = "a"),
  set_prior("logistic(-1, 0.5)", coef = "Intercept", dpar = "ndt"),
  set_prior("logistic(0, 0.67)", coef = "Intercept", dpar = "w"),
  set_prior("logistic(-0.5, 0.5)", class = "Intercept", dpar = "sv"),
  set_prior("student_t(4, 0, 0.5)", class = "b"),
  set_prior("student_t(4, 0, 0.5)", class = "b", dpar = "a"),
  set_prior("student_t(4, 0, 0.1)", class = "b", dpar = "ndt"),
  set_prior("student_t(4, 0, 0.5)", class = "b", dpar = "w")
)

tmp_dat <- make_standata(formula,
                         family = ddm,
                         data = speed_acc,
                         prior = priors)

init_fun <- function() {
  list(
    b = rnorm(tmp_dat$K, 0, 0.25),
    b_a = c(rnorm(1, log(1), 0.1), rnorm(tmp_dat$K_a-1, 0, 0.1)),
    b_ndt = c(runif(1, log(0.1), log(0.2)), rnorm(tmp_dat$K_ndt-1, 0, 0.01)),
    b_w = c(runif(1, logit_scaled(0.4), logit_scaled(0.6)), rnorm(tmp_dat$K_w-1, 0, 0.01)),
    Intercept_sv = rnorm(1, log(0.1), 0.01)
  )
}
```

Now we run the sampling.

```{r rw-sampling5, eval=FALSE, include=TRUE}
fit_rwdata_5par <- brm(formula,
                       family = ddm,
                       prior = priors,
		       inits = init_fun,
                       stanvars = stanvars_rw,
                       data = speed_acc,
                       chains = 4,
                       iter = 1000,
                       warmup = 500,
                       thin = 1,
                       cores = 4,
                       control = list(
                         max_treedepth = 10,
                         adapt_delta = 0.8
                       )
)
```

```{r rw-check-if-already-fit5, include=FALSE}
if (file.exists(file.path(getwd(), "fits", "fit-rwdata-5par_prior-on-intercepts-and-bs_2-chains.RDS"))) {
  fit_rwdata_5par <- readRDS(file = file.path(getwd(), "fits",
                                              "fit-rwdata-5par_prior-on-intercepts-and-bs_2-chains.RDS"))
} else { # run fit
  fit_rwdata_5par <- brm(formula,
                         family = ddm,
                         prior = priors,
                         stanvars = stanvars_rw,
                         data = speed_acc,
                         chains = 4,
                         iter = 1000,
                         warmup = 500,
                         thin = 1,
                         cores = 4,
                         control = list(
                           max_treedepth = 10,
                           adapt_delta = 0.8
                         )
)
  saveRDS(fit_rwdata_5par,
        file = file.path(getwd(),"fits", "fit_rwdata_5par.RDS"))
}
```

```{r rw-save-fit5, eval=FALSE, include=TRUE}
saveRDS(fit_rwdata_5par,
        file = file.path(getwd(),"fits", "fit_rwdata_5par.RDS"))
```



## Posterior Predictive Checks {#rewrld-post}

LOO

```{r rw-loo}
loo(fit_rwdata_5par)
```

Posterior predictive check

```{r rw-pp}
pp_check(fit_rwdata_5par)
```





</div>
# {.unlisted .unnumbered}
#### R Session Info {.unlisted .unnumbered}
```{r session-info, collapse=TRUE}
sessionInfo()
```

```{r reset-options, include=FALSE}
options(op)  # reset options
```


# References