Function dfddm() evaluates the density function (or probability density function, PDF) for the Ratcliff diffusion decision model (DDM) using different methods for approximating the full PDF, which contains an infinite sum. An overview of the mathematical details of the different approximations is provided in the Math Vignette. An empirical validation of the implemented methods is provided in the Validity Vignette. Examples of using dfddm() for parameter estimation are provided in the Example Vignette.

Our implementation of the DDM has the following 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), \(sv \in (0, \infty)\) (inter-trial-variability of drift), and \(\sigma \in (0, \infty)\) (diffusion coefficient of the underlying Wiener Process).

Introduction

Since the DDM is widely used in parameter estimation usually involving numerical optimization, significant effort has been put into making the evaluation of its density as fast as possible. However, the density function for the DDM is notorious for containing an unavoidable infinite sum; hence, the literature has produced a few different methods of approximating the density. This vignette is designed to compare the speed of the currently available implementations that calculate the lower probability density of the DDM. To this aim, we present benchmark data of currently available approximation methods from the literature, our streamlined implementations of those approximation methods, and our own novel approximation method. The first section will record the benchmark data of our implementations of the density function approximations and show the results through a series of visualizations. The second section will record the benchmark data of parameter estimation that uses the density function approximations in the optimization process for fitting to real-world data; this section will also include visualizations to illustrate the differences among the density function approximations. Please note that we will not be testing the functions for accuracy or consistency in this vignette, as that is covered in the Validity Vignette.

Benchmarking the Density Function Approximations

We want to determine the performance of each implementation across a variety of parameters in order to identify any slow areas for individual implementations. To achieve this rigor, we define a parameter space (in a code chunk below) and loop through each combination of parameters. For the preliminary benchmark testing, we will input the response times as a vector to each function instead of inputting individual response times for two reasons: first, this is the most common way to input the response time data; and second, this allows for the R-based implementation to exploit vectorization to make that implementation as fast as possible. For each combination of parameters, we run the microbenchmark::microbenchmark() function from the package microbenchmark 1000 times for each implementation and only save the median benchmark time of these 1000. We will refer to these medians as simply the “benchmark times” for the remainder of this vignette. Running the benchmark tests this many times for each set of parameters generates a sufficient amount of benchmark data so that the median of these times is unlikely to be an outlier, either faster or slower than a typical run time.

Upon collecting the benchmark data, we perform a preliminary analysis of the implementations’ performance by showing the distributions of benchmark times as side-by-side violin plots for each implementation. Then we cull the implementations and only test further on the fastest and most widely used implementations. To elucidate the areas of the parameter space where the different implementations excel and struggle, we plot the median benchmark times as a function of the response time input. This plot will reveal the efficiencies of the various implementations that we are testing.

Generating Benchmark Data

Before analyzing how the different density function approximations perform, we must first generate the benchmark data. In addition to testing all of the implementations available in dfddm(), we also include three functions that are considered to be current and widely used. First, we include the function rtdists::ddiffusion() from the package rtdists as it is well known for being not only user friendly and feature-rich but also designed specifically for handling data with regard to distributions for response time models. Second, we also test the RWiener::dwiener() function from the RWiener package, which is mainly aimed at providing an R language interface to calculate various functions from the Wiener process. Third, we include some raw R code that is bundled with the paper written by Gondan, Blurton, and Kesselmeier (2014) about improving the approximation to the density function. As this code was not yet available in an R package, it is part of fddm and can be accessed after running source on the corresponding file as shown below.

To capture the behavior of the density function approximations, we test each one in a granular parameter space that includes the parameter values most likely used in practice. This parameter space can be found fully defined in the code chunk later in this section. Note that we are only testing the “lower” density function as the conversion from the “lower” density function to the “upper” density function involves only the mappings \(v \to -v\) and \(w \to 1-w\). Since our parameter space includes the negatives of all positive values of \(v\) and the complements of all values of \(w\), testing both the “lower” and “upper” density function approximations would be redundant.

In addition, we include two functions for benchmarking the performance of the implementations: one where the response times are input as a vector, and one where the response times are input one at a time to the density function. Inputting the response times as a vector is the most common way of inputting data to the density function and allows the R code from Gondan, Blurton, and Kesselmeier (2014) to exploit R’s powerful vectorization, operating at its highest efficiency. We also include a benchmark function that inputs the response times individually as this allows us to see the impact that the size of the input response time has on the evaluation time of the density function approximations.

At each point in this parameter space, we run the microbenchmark::microbenchmark() function 1000 times for each implementation and save the median of these 1000 benchmark times. Using the median of these data should reduce the number of outliers in the data by limiting the influence of computational noise on the results. Note that for expositional purposes we will only say “benchmark times” to mean the collection of medians of the 1000 microbenchmark::microbenchmark() measurements recorded at each point in the parameter space. The following code chunk defines the functions that we will use to benchmark the implementations; we will be using the microbenchmark package to time the evaluations.

library("fddm")
library("rtdists")
library("RWiener")
source(system.file("extdata", "Gondan_et_al_density.R",
                   package = "fddm", mustWork = TRUE))
library("microbenchmark")

rt_benchmark_vec <- function(RT, resp, V, A, t0 = 1e-4, W = 0.5,
                             err_tol = 1e-6, times = 100, unit = "ns") {

  fnames <- c("fs_SWSE_17", "fs_SWSE_14", "ft_SWSE_17", "ft_SWSE_14",
              "fb_SWSE_17", "fb_SWSE_14",
              "fs_Gon_17", "fs_Gon_14", "fb_Gon_17", "fb_Gon_14",
              "fs_Nav_17", "fs_Nav_14", "fb_Nav_17", "fb_Nav_14",
              "fl_Nav_09", "RWiener", "Gondan", "rtdists")
  nf <- length(fnames) # number of functions being benchmarked
  nV <- length(V)
  nA <- length(A)
  nW <- length(W)
  resp <- rep(resp, length(RT)) # for RWiener

  # Initialize the dataframe to contain the microbenchmark results
  mbm_res <- data.frame(matrix(ncol = 3+nf, nrow = nV*nA*nW))
  colnames(mbm_res) <- c('V', 'A', 'W', fnames)
  row_idx <- 1

  # Loop through each combination of parameters and record microbenchmark results
  for (v in 1:nV) {
    for (a in 1:nA) {
      for (w in 1:nW) {
        mbm <- microbenchmark(
        fs_SWSE_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                           t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                           switch_mech = "small", n_terms_small = "SWSE",
                           summation_small = "2017"),
        fs_SWSE_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                           t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                           switch_mech = "small", n_terms_small = "SWSE",
                           summation_small = "2014"),
        ft_SWSE_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                           t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                           switch_mech = "eff_rt", switch_thresh = 0.8,
                           n_terms_small = "SWSE", summation_small = "2017"),
        ft_SWSE_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                           t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                           switch_mech = "eff_rt", switch_thresh = 0.8,
                           n_terms_small = "SWSE", summation_small = "2014"),
        fb_SWSE_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                           t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                           switch_mech = "terms_large", switch_thresh = 0.8,
                           n_terms_small = "SWSE", summation_small = "2017"),
        fb_SWSE_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                           t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                           switch_mech = "terms_large", switch_thresh = 0.8,
                           n_terms_small = "SWSE", summation_small = "2014"),
        fs_Gon_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                          switch_mech = "small", n_terms_small = "Gondan",
                          summation_small = "2017"),
        fs_Gon_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                          switch_mech = "small", n_terms_small = "Gondan",
                          summation_small = "2014"),
        fb_Gon_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                          switch_mech = "terms", n_terms_small = "Gondan",
                          summation_small = "2017"),
        fb_Gon_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                          switch_mech = "terms", n_terms_small = "Gondan",
                          summation_small = "2014"),
        fs_Nav_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                          switch_mech = "small", n_terms_small = "Navarro",
                          summation_small = "2017"),
        fs_Nav_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                          switch_mech = "small", n_terms_small = "Navarro",
                          summation_small = "2014"),
        fb_Nav_17 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                          switch_mech = "terms", n_terms_small = "Navarro",
                          summation_small = "2017"),
        fb_Nav_14 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                          switch_mech = "terms", n_terms_small = "Navarro",
                          summation_small = "2014"),
        fl_Nav_09 = dfddm(rt = RT, response = resp, a = A[a], v = V[v],
                          t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                          switch_mech = "large"),
        RWiener = dwiener(RT, resp = resp, alpha = A[a],
                          delta = V[v], tau = t0, beta = W[w],
                          give_log = FALSE),
        Gondan = fs(t = RT-t0, a = A[a], v = V[v],
                    w = W[w], eps = err_tol), # only "lower" resp
        rtdists = ddiffusion(RT, resp, a = A[a], v = V[v],
                              t0 = t0, z = W[w]*A[a]),
        times = times, unit = unit)
        # add the v, a, and w values to the dataframe
        mbm_res[row_idx, 1] <- V[v]
        mbm_res[row_idx, 2] <- A[a]
        mbm_res[row_idx, 3] <- W[w]
        # add the median microbenchmark results to the dataframe
        for (i in 1:nf) {
          mbm_res[row_idx, 3+i] <- median(mbm[mbm[,1] == fnames[i],2])
        }
        # iterate start value
        row_idx <- row_idx + 1
      }
    }
  }
  return(mbm_res)
}

rt_benchmark_ind <- function(RT, resp, V, A, t0 = 1e-4, W = 0.5,
                             err_tol = 1e-6, times = 100, unit = "ns") {
  fnames <- c("fs_SWSE_17", "fs_SWSE_14", "ft_SWSE_17", "ft_SWSE_14",
              "fb_SWSE_17", "fb_SWSE_14",
              "fs_Gon_17", "fs_Gon_14", "fb_Gon_17", "fb_Gon_14",
              "fs_Nav_17", "fs_Nav_14", "fb_Nav_17", "fb_Nav_14",
              "fl_Nav_09", "RWiener", "Gondan", "rtdists")
  nf <- length(fnames) # number of functions being benchmarked
  nRT <- length(RT)
  nV <- length(V)
  nA <- length(A)
  nW <- length(W)

  # Initialize the dataframe to contain the microbenchmark results
  mbm_res <- data.frame(matrix(ncol = 4+nf, nrow = nRT*nV*nA*nW))
  colnames(mbm_res) <- c('RT', 'V', 'A', 'W', fnames)
  row_idx <- 1

  # Loop through each combination of parameters and record microbenchmark results
  for (rt in 1:nRT) {
    for (v in 1:nV) {
      for (a in 1:nA) {
        for (w in 1:nW) {
          mbm <- microbenchmark(
          fs_SWSE_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                             t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                             switch_mech = "small", n_terms_small = "SWSE",
                             summation_small = "2017"),
          fs_SWSE_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                             t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                             switch_mech = "small", n_terms_small = "SWSE",
                             summation_small = "2014"),
          ft_SWSE_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                             t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                             switch_mech = "eff_rt", switch_thresh = 0.8,
                             n_terms_small = "SWSE", summation_small = "2017"),
          ft_SWSE_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                             t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                             switch_mech = "eff_rt", switch_thresh = 0.8,
                             n_terms_small = "SWSE", summation_small = "2014"),
          fb_SWSE_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                             t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                             switch_mech = "terms_large", switch_thresh = 0.8,
                             n_terms_small = "SWSE", summation_small = "2017"),
          fb_SWSE_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                             t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                             switch_mech = "terms_large", switch_thresh = 0.8,
                             n_terms_small = "SWSE", summation_small = "2014"),
          fs_Gon_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                            t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                            switch_mech = "small", n_terms_small = "Gondan",
                            summation_small = "2017"),
          fs_Gon_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                            t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                            switch_mech = "small", n_terms_small = "Gondan",
                            summation_small = "2014"),
          fb_Gon_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                            t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                            switch_mech = "terms", n_terms_small = "Gondan",
                            summation_small = "2017"),
          fb_Gon_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                            t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                            switch_mech = "terms", n_terms_small = "Gondan",
                            summation_small = "2014"),
          fs_Nav_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                            t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                            switch_mech = "small", n_terms_small = "Navarro",
                            summation_small = "2017"),
          fs_Nav_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                            t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                            switch_mech = "small", n_terms_small = "Navarro",
                            summation_small = "2014"),
          fb_Nav_17 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                            t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                            switch_mech = "terms", n_terms_small = "Navarro",
                            summation_small = "2017"),
          fb_Nav_14 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                            t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                            switch_mech = "terms", n_terms_small = "Navarro",
                            summation_small = "2014"),
          fl_Nav_09 = dfddm(rt = RT[rt], response = resp, a = A[a], v = V[v],
                            t0 = t0, w = W[w], err_tol = err_tol, log = FALSE,
                            switch_mech = "large"),
          RWiener = dwiener(RT[rt], resp = resp, alpha = A[a],
                            delta = V[v], tau = t0, beta = W[w],
                            give_log = FALSE),
          Gondan = fs(t = RT[rt]-t0, a = A[a], v = V[v],
                      w = W[w], eps = err_tol), # only "lower" resp
          rtdists = ddiffusion(RT[rt], resp, a = A[a], v = V[v],
                                t0 = t0, z = W[w]*A[a]),
          times = times, unit = unit)
          # add the v, a, and w values to the dataframe
          mbm_res[row_idx, 1] <- RT[rt]
          mbm_res[row_idx, 2] <- V[v]
          mbm_res[row_idx, 3] <- A[a]
          mbm_res[row_idx, 4] <- W[w]
          # add the median microbenchmark results to the dataframe
          for (i in 1:nf) {
            mbm_res[row_idx, 4+i] <- median(mbm[mbm[,1] == fnames[i],2])
          }
          # iterate start value
          row_idx <- row_idx + 1
        }
      }
    }
  }
  return(mbm_res)
}

Control loops in R have a lot of overhead and are comparatively slow compared to loops in C++; thus, we will not run rt_benchmark_vec() in this vignette, and we will load pre-run benchmark data instead. The pre-run data was collected using the same functions and parameter space as defined in this vignette. Since this benchmark function contains one fewer loop because of the vectorized response times, we increase the number of samples that the microbenchmark::microbenchmark() function considers when timing the various implementations. To show comparisons across the implementations with respect to the input response time, we will again load pre-run benchmark data obtained through using rt_benchmark_ind() as this takes even longer to evaluate than the benchmark with vectorized response times. The following code chunk defines the parameter space to be explored throughout the benchmark testing and shows an example of how to run the benchmark functions defined in the previous code chunk.

# Define parameter space
RT <- seq(0.1, 2, by = 0.1)
A <- seq(0.5, 3.5, by = 0.5)
V <- c(-5, -2, 0, 2, 5)
t0 <- 1e-4 # must be nonzero for RWiener
W <- seq(0.3, 0.7, by = 0.1)
err_tol <- 1e-6 # this is the setting from rtdists

# Run benchmark tests
bm_vec <- rt_benchmark_vec(RT = RT, resp = "lower", V = V, A = A, t0 = t0,
                           W = W, err_tol = err_tol,
                           times = 1000, unit = "ns")
bm_ind <- rt_benchmark_ind(RT = RT, resp = "lower", V = V, A = A, t0 = t0,
                           W = W, err_tol = err_tol,
                           times = 100, unit = "ns")

Analysis of Benchmark Results

The first step in analyzing the benchmark data is a rough visualization of the results in side-by-side violin plots that show the distribution of microbenchmark::microbenchmark() timings. We use the data that we just generated in the previous code chunk to plot a density violin for each implementation. Each density violin is overlaid with a boxplot showing the median and the first and third quartiles, in addition to a dashed line representing the mean. We also zoom in to the area of the plot where most of the fddm implementations lie; depending on your system, this may obscure the results from rtdists as they are typically much higher than that of fddm.

For this visualization we will use the benchmark data where the response times were handed over as a vector to each fddm call. As discussed above, running the benchmark on individual response times is considerably slower as the loop over the response times resides in R and not in C++.

library("reshape2")
library("ggplot2")
library("ggforce")

# load data, will be in the variable 'bm_vec'
load(system.file("extdata", "dfddm_density", "bm_vec_0-2.Rds",
                 package = "fddm", mustWork = TRUE))

t_idx <- match("W", colnames(bm_vec))
bm_vec[, -seq_len(t_idx)] <- bm_vec[, -seq_len(t_idx)]/1000 # convert to microseconds
mbm_vec <- melt(bm_vec, measure.vars = -seq_len(t_idx),
                variable.name = "FuncName", value.name = "time")

Names_vec <- c("ft_SWSE_17", "ft_SWSE_14", "fb_SWSE_17", "fb_SWSE_14",
               "fb_Gon_17", "fb_Gon_14", "fb_Nav_17", "fb_Nav_14",
               "fs_SWSE_17", "fs_SWSE_14",
               "fs_Gon_17", "fs_Gon_14", "fs_Nav_17", "fs_Nav_14",
               "fl_Nav_09", "RWiener", "Gondan", "rtdists")
Color_vec <- c("#92c639", "#d3e8b0", "#729438", "#d2e3b5",
               "#b3724d", "#e0c7b8", "#4da7b3", "#b8dce0",
               "#5cc639", "#bee8b0",
               "#b34d4d", "#e0b8b8", "#4d80b3", "#b8cce0",
               "#dcdca3", "#deccba", "#c5a687", "#ac8053")
Outline_vec <- c("#92c639", "#92c639", "#729438", "#729438",
                 "#b3724d", "#b3724d", "#4da7b3", "#4da7b3",
                 "#5cc639", "#5cc639",
                 "#b34d4d", "#b34d4d", "#4d80b3", "#4d80b3",
                 "#dcdca3", "#deccba", "#c5a687", "#ac8053")

mi <- min(bm_vec[, -seq_len(t_idx)])
ma <- max(bm_vec[, (t_idx+1):(ncol(bm_vec)-4)])

ggplot(mbm_vec, aes(x = factor(FuncName, levels = Names_vec), y = time,
                    color = factor(FuncName, levels = Names_vec),
                    fill = factor(FuncName, levels = Names_vec))) +
  geom_violin(trim = TRUE, alpha = 0.5) +
  scale_color_manual(values = Outline_vec, guide = "none") +
  scale_fill_manual(values = Color_vec, guide = "none") +
  geom_boxplot(width = 0.15, fill = "white", alpha = 0.5) +
  stat_summary(fun = mean, geom = "errorbar",
               aes(ymax = after_stat(y), ymin = after_stat(y)),
               width = .35, linetype = "dashed") +
  scale_x_discrete(labels = c(
    bquote(f[t] ~ SWSE[17]), bquote(f[t] ~ SWSE[14]),
    bquote(f[c] ~ SWSE[17]), bquote(f[c] ~ SWSE[14]),
    bquote(f[c] ~ Gon[17]), bquote(f[c] ~ Gon[14]),
    bquote(f[c] ~ Nav[17]), bquote(f[c] ~ Nav[14]),
    bquote(f[s] ~ SWSE[17]), bquote(f[s] ~ SWSE[14]),
    bquote(f[s] ~ Gon[17]), bquote(f[s] ~ Gon[14]),
    bquote(f[s] ~ Nav[17]), bquote(f[s] ~ Nav[14]),
    bquote(f[l] ~ "Nav"), "RWiener", "Gondan", "rtdists")) +
  facet_zoom(ylim = c(mi, ma)) +
  labs(x = "Implementation", y = "Time (microseconds)") +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 16, angle = 90,
                                   vjust = 0.5, hjust = 1),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20,
                                    margin = margin(10, 0, 0, 0)),
        axis.title.y = element_text(size = 20,
                                    margin = margin(0, 10, 0, 0)),
        legend.position = "none")

Depending on your system, the results might differ from the ones shown here, but they should replicate the general pattern. The small-time implementations (\(f_s\)) have long tails due to very few outliers, but their central tendencies are rather fast. The implementations combining small-time and large-time (\(f_c\)) approximations show an almost normal distribution with only a short tail; they also have rather fast central tendencies. Of the fddm implementations, the large-time implementation (\(f_\ell\)) is the slowest with a long tail. The function from the RWiener package combines small-time and large-time approach and also shows an almost normal distribution, but its distribution of median benchmark times is shifted higher compared to those of fddm. The other two implementations are clearly slower; the function provided by Gondan, Blurton, and Kesselmeier (2014) is purely in R so it is constrained by the speed of R’s vectorization. Similarly constrained by R, the function from the rtdists package performs quite a few computations in R that make it very user-friendly but at the cost of having the slowest central tendencies.

To visualize how the response time input can affect the efficiency of each implementation, we plot the benchmark times as a function of the effective response time. We define the effective response time as \(\frac{t}{a^2}\) because the two model parameters \(t\) and \(a\) always appear inside the infinite sums in this relationship. In the following plot, we will vary the effective response time and show the mean of the benchmark times, the 10% quantiles and 90% quantiles of the benchmark times, as well as the minimum and maximum benchmark times. For each value of the effective response time, we aggregate the benchmark data for all values of the other parameter to generate the statistics plotted. Note that because we vary the response time in this plot, we use the benchmark data generated from rt_benchmark_ind). As this benchmark function records the evaluation time of approximating the density for only one response time in contrast to a vector of response times, we expect the typical benchmark times produced by rt_benchmark_ind to be faster than that of rt_benchmark_vec.

To avoid too many repetitious plots, we will only include six of the implementations from the previous plot that illustrate the different variations and implementations of the density function approximations. We include: the fastest small-time implementation in dfddm() (“\(f_s \text{SWSE}_{17}\)”), the one large-time implementation of the large-time in dfddm (“\(f_\ell\) Nav”), the fastest implementation in dfddm() that combines the small-time and large-time (“\(f_c \text{SWSE}_{17}\)”), and the three implementations from R packages and the literature as described above.

The following code chunk loads the pre-run benchmark data, prepares it, and plots the benchmark times as the effective response time varies.

# load data, will be in the variable 'bm_ind'
load(system.file("extdata", "dfddm_density", "bm_ind_0-2.Rds",
                 package = "fddm", mustWork = TRUE))
bm_ind[["RTAA"]] <- bm_ind[["RT"]] / bm_ind[["A"]] / bm_ind[["A"]]
bm_ind <- bm_ind[, c(1, 2, ncol(bm_ind), 3:(ncol(bm_ind)-1)) ]

t_idx <- match("W", colnames(bm_ind))
bm_ind[,-seq_len(t_idx)] <- bm_ind[, -seq_len(t_idx)]/1000 # convert to microseconds
mbm_ind <- melt(bm_ind, measure.vars = -seq_len(t_idx),
                variable.name = "FuncName", value.name = "time")

Names_meq <- c("ft_SWSE_17", "fs_SWSE_17", "fl_Nav_09",
               "RWiener", "Gondan", "rtdists")
Color_meq <- c("#92c639", "#5cc639", "#dcdca3",
               "#deccba", "#c5a687", "#ac8053")

mbm_meq <- subset(mbm_ind, FuncName %in% Names_meq)

my_labeller <- as_labeller(c(ft_SWSE_17 = "f[t] ~ SWSE[17]",
                             fs_SWSE_17 = "f[s] ~ SWSE[17]",
                             fl_Nav_09 = "f[l] ~ Nav",
                             RWiener = "RWiener",
                             Gondan = "Gondan",
                             rtdists = "rtdists"),
                           default = label_parsed)

ggplot(mbm_meq, aes(x = RTAA, y = time,
                    color = factor(FuncName, levels = Names_meq),
                    fill = factor(FuncName, levels = Names_meq))) +
  stat_summary(fun.min = min, fun.max = max,
               geom = "ribbon", color = NA, alpha = 0.1) +
  stat_summary(fun.min = function(z) { quantile(z, 0.1) },
               fun.max = function(z) { quantile(z, 0.9) },
               geom = "ribbon", color = NA, alpha = 0.2) +
  stat_summary(fun = mean, geom = "line") +
  scale_x_log10(breaks = c(0.001, 0.1, 10, 1000),
                labels = as.character(c(0.001, 0.1, 10, 1000))) +
  scale_color_manual(values = Color_meq) +
  scale_fill_manual(values = Color_meq) +
  labs(subtitle = paste(
         "The darker shaded regions represent the 10% and 90% quantiles",
         "The lighter shaded regions represent the min and max times",
         sep = ";\n"),
       x = bquote(frac(t, a^2) ~ ", effective response time, " ~ log[10]),
       y = "Time (microseconds)") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        plot.subtitle = element_text(size = 16,
                                     margin = margin(0, 0, 15, 0)),
        axis.text.x = element_text(size = 16, angle = 90,
                                   vjust = 0.5, hjust = 1),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20,
                                    margin = margin(10, 0, 0, 0)),
        axis.title.y = element_text(size = 20,
                                    margin = margin(0, 10, 0, 0)),
        strip.text = element_text(size = 16),
        strip.background = element_rect(fill = "white"),
        legend.position = "none") +
  facet_wrap(~ factor(FuncName, levels = Names_meq), scales = "free_y",
             labeller = my_labeller)

This plot confirms that the aptly named small-time and large-time density functions are generally more efficient at calculating the density for small response times and large response times, respectively. Of the three implementations used from the literature, the RWiener and rtdists functions use a switching mechanism between the small-time and large-time density functions, and the implementation from Gondan, Blurton, and Kesselmeier (2014) only uses the small-time density function. Importantly, the implementations that utilize a switching mechanism between the small-time and large-time density function approximations have much flatter plots and generally perform well across the \(\frac{t}{a^2}\) parameter space; this confirms that using a combination of both timescales results in a very stable implementation.

Benchmarking Model Fitting to Real-World Data

Having shown the performance of each implementation of the Ratcliff DDM density function approximations, we now examine how these implementations compare in practice. To determine how each implementation behaves in a practical setting, we will use each implementation as the basis for fitting the DDM to real-world data. The primary analysis of the implementations’ performance will come from benchmarking the time taken for the parameter estimation using the different implementations. In addition, a secondary analysis will compare the performance of the different implementations by using the number of function evaluations used by the optimizer during the optimization process.

The parameters that we will fit are: \(a\) (threshold separation), \(v\) (drift rate), \(t_0\) (non-decision time/response time constant), \(w\) (relative starting point), and \(sv\) (inter-trial-variability of drift). Since the example data we are using consists of two different item types that each require a different correct response (i.e., for one item type the correct response is mapped to the lower boundary and for the other item type the correct response is mapped to the upper boundary), the model includes two different versions of \(v\) (drift rate): \(v_\ell\) for fitting to the truthful lower boundary, and \(v_u\) for fitting to the truthful upper boundary. After using the various density function approximations in fitting the DDM to real-world data, we then validate that the produced parameter estimates are consistent across the various implementations within a given error tolerance.

We will not test for accuracy or consistency in this vignette because that has already been covered in the Validity Vignette and the optimization is deterministic so it will return the same result given the same inital setup.

Generating Benchmark Data for Parameter Estimation

The following subsections will define all of the functions used to generate the fittings, store the benchmark results, and provide the code to run the full fitting. Since running the full fitting for all of the individuals in the data takes a long time, we will forgo actually running the fitting in this vignette and instead load pre-fit parameter estimates that used the provided code.

While we can use every implementation in dfddm(), we cannot use the implementation from the RWiener package nor the R implementation from Gondan, Blurton, and Kesselmeier (2014) because neither of these two implementations contains inter-trial variability of the drift rate (i.e., the model parameter \(sv\)) in their estimation of the DDM density function. We can, however, use the implementation from rtdists (rtdists::ddiffusion()) as this implementation does include the variable drift rate versions of the density functions.While we can convert the constant drift rate density to the variable drift rate density using a multiplicative factor, the densities are still potentially calculated incorrectly. For more details on the differences between the constant drift rate density functions and their variable drift rate counterparts, see the Math Vignette.

The optimization routine that we use is the R function stats::nlminb() that uses the maximum likelihood method and a standard gradient-descent optimization algorithm. Since stats::nlminb() minimizes the likelihood function instead of maximizing it, we define a family of likelihood functions that return the negative sum of the log-likelihoods of the data. In particular, each likelihood function in this family uses a different implementation to evaluate the log-likelihoods of the data. This way we can use the same optimization process for each implementation, and the only difference in benchmark data arises from the underlying implementations that approximate the DDM density function.

First we load the necessary packages.

library("fddm")
library("rtdists")
library("microbenchmark")

Log-Likelihood Functions

This code chunk defines the log-likelihood functions used in the optimization algorithm; all of the log-likelihood functions are the same as in the Validity Vignette. The log-likelihood functions are fairly straightforward and split the responses and associated response times by the true item status (i.e., the correct response) to enable fitting distinct drift rates (\(v_\ell\) for the items in which the correct response is the lower boundary and \(v_u\) for the items for which the correct response is the upper boundary). In addition, the log-likelihood functions heavily penalize any combination of parameters that returns a log-density of \(-\infty\) (equivalent to a regular density of \(0\)).

ll_ft_SWSE_17 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "eff_rt", switch_thresh = 0.8,
                n_terms_small = "SWSE", summation_small = "2017")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_ft_SWSE_14 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "eff_rt", switch_thresh = 0.8,
                n_terms_small = "SWSE", summation_small = "2014")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fb_SWSE_17 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "terms_large", switch_thresh = 0.8,
                n_terms_small = "SWSE", summation_small = "2017")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fb_SWSE_14 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "terms_large", switch_thresh = 0.8,
                n_terms_small = "SWSE", summation_small = "2014")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fb_Gon_17 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "terms", n_terms_small = "Gondan",
                summation_small = "2017")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fb_Gon_14 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "terms", n_terms_small = "Gondan",
                summation_small = "2014")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fb_Nav_17 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "terms", n_terms_small = "Navarro",
                summation_small = "2017")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fb_Nav_14 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "terms", n_terms_small = "Navarro",
                summation_small = "2014")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fs_SWSE_17 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "small", n_terms_small = "SWSE",
                summation_small = "2017")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fs_SWSE_14 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "small", n_terms_small = "SWSE",
                summation_small = "2014")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fs_Gon_17 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "small", n_terms_small = "Gondan",
                summation_small = "2017")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fs_Gon_14 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "small", n_terms_small = "Gondan",
                summation_small = "2014")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fs_Nav_17 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "small", n_terms_small = "Navarro",
                summation_small = "2017")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fs_Nav_14 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "small", n_terms_small = "Navarro",
                summation_small = "2014")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_fl_Nav_09 <- function(pars, rt, resp, truth, err_tol) {
  v <- numeric(length(rt))
  v[truth == "upper"] <- pars[[1]]
  v[truth == "lower"] <- pars[[2]]
  dens <- dfddm(rt = rt, response = resp, a = pars[[3]], v = v,
                t0 = pars[[4]], w = pars[[5]], sv = pars[[6]], err_tol = 1e-6,
                log = TRUE, switch_mech = "large")
  return( ifelse(any(!is.finite(dens)), 1e6, -sum(dens)) )
}

ll_RTDists <- function(pars, rt, resp, truth) {
  rtu <- rt[truth == "upper"]
  rtl <- rt[truth == "lower"]
  respu <- resp[truth == "upper"]
  respl <- resp[truth == "lower"]

  densu <- ddiffusion(rtu, respu, a = pars[[3]], v = pars[[1]],
                      z = pars[[5]]*pars[[3]], t0 = pars[[4]], sv = pars[[6]])
  densl <- ddiffusion(rtl, respl, a = pars[[3]], v = pars[[2]],
                      z = pars[[5]]*pars[[3]], t0 = pars[[4]], sv = pars[[6]])

  densities <- c(densu, densl)
  if (any(densities <= 0)) return(1e6)
  return(-sum(log(densities)))
}

Fitting Function

As encountering local optima is a known problem for deterministic optimization algorithms, we run stats::nlminb() with eleven different sets of initial values. In a real-life situation we would probably use random initial parameter values to avoid local minima; however, we use a fixed set of initial parameter values to verify that the different underlying implementations produce the same results for different starting points.

There are two restrictions to these initial values that we must address (we would also have to do so if we picked the initial values randomly). First, the parameter \(t_0\) must be less than the response time, so we set the initial values for \(t_0\) to be strictly less than the minimum response time for to each individual in the input data frame. Second, we must place a lower bound on \(a\) that is necessarily greater than zero because the optimization algorithm occasionally evaluates the log-likelihood functions (and thus the underlying density function approximations) using values of \(a\) equal to its bounds. In the case where optimization algorithm evaluates using \(a = 0\), both the density function from the rtdists package and every implementation from fddm do not evaluate. In common use this is not an issue because very small values of \(a\) do not make any sense with regard to the psychological interpretation of the parameter, but this issue can arise in an exploratory optimization environment.

The following code chunk defines the function that will run the optimization and produce the fitted parameter estimates for \(v_u\), \(v_\ell\), \(a\), \(t_0\), \(w\), and \(sv\). As discussed above, we only use a selection of implementations from those available in dfddm(). This fitting function will run the optimization for each set of initial parameter values and store the following results from the optimization: 1) resulting convergence code (either \(0\) indicating successful convergence or \(1\) indicating unsuccessful convergence); 2) minimized value of the objective log-likelihood function; 3) the number of iterations used by the optimizer during the optimization process; 4) the number of calls to the log-likelihood function during the optimization process; and 5) the median benchmark time of the optimization using each implementation.

As in the previous section, we have the microbenchmark::microbenchmark() function run the data fitting process multiple times for each implementation at each set of initial parameter values, and we save the median of these benchmark times. However, we only use a smaller number of repetitions here because a full optimization run entails several calls to the underlying density function and takes much longer to run. Using the median of these data reduces the impact of computational noise on the benchmark results, and again we will refer to these medians as simply the “benchmark times”. We will be using the microbenchmark package to time the evaluations.

rt_fit <- function(data, id_idx = NULL, rt_idx = NULL, response_idx = NULL,
                   truth_idx = NULL, response_upper = NULL, err_tol = 1e-6,
                   ctrl_list = list(eval.max = 1000, iter.max = 1000),
                   times = 100, unit = "ns") {

  # Format data for fitting
  if (all(is.null(id_idx), is.null(rt_idx), is.null(response_idx),
      is.null(truth_idx), is.null(response_upper))) {
    df <- data # assume input data is already formatted
  } else {
    if(any(data[,rt_idx] < 0)) {
      stop("Input data contains negative response times; fit will not be run.")
    }
    if(any(is.na(data[,response_idx]))) {
      stop("Input data contains invalid responses (NA); fit will not be run.")
    }

    nr <- nrow(data)
    df <- data.frame(id = character(nr),
                     rt = double(nr),
                     response = character(nr),
                     truth = character(nr),
                     stringsAsFactors = FALSE)

    if (!is.null(id_idx)) { # relabel identification tags
      for (i in seq_along(id_idx)) {
        idi <- unique(data[,id_idx[i]])
        for (j in seq_along(idi)) {
          df[["id"]][data[,id_idx[i]] == idi[j]] <- paste(
            df[["id"]][data[,id_idx[i]] == idi[j]], idi[j], sep = " ")
        }
      }
      df[["id"]] <- trimws(df[["id"]], which = "left")
    }

    df[["rt"]] <- as.double(data[,rt_idx])

    df[["response"]] <- "lower"
    df[["response"]][data[,response_idx] == response_upper] <- "upper"

    df[["truth"]] <- "lower"
    df[["truth"]][data[,truth_idx] == response_upper] <- "upper"
  }

  # Preliminaries
  ids <- unique(df[["id"]])
  nids <- max(length(ids), 1) # if inds is null, there is only one individual

  init_vals <- data.frame(v1 = c( 0,  10, -.5,  0,  0,  0,  0,  0,  0,   0,  0),
                          v0 = c( 0, -10,  .5,  0,  0,  0,  0,  0,  0,   0,  0),
                          a  = c( 1,   1,   1, .5,  5,  1,  1,  1,  1,   1,  1),
                          t0 = c( 0,   0,   0,  0,  0,  0,  0,  0,  0,   0,  0),
                          w  = c(.5,  .5,  .5, .5, .5, .5, .5, .2, .8,  .5, .5),
                          sv = c( 1,   1,   1,  1,  1,  1,  1,  1,  1, .05,  5))
  ninit_vals <- nrow(init_vals)

  algo_names <- c("ft_SWSE_17", "ft_SWSE_14", "fb_SWSE_17", "fb_SWSE_14",
                  "fb_Gon_17", "fb_Gon_14", "fb_Nav_17", "fb_Nav_14",
                  "fs_SWSE_17", "fs_SWSE_14", "fs_Gon_17", "fs_Gon_14",
                  "fs_Nav_17", "fs_Nav_14", "fl_Nav_09", "rtdists")
  nalgos <- length(algo_names)
  ni <- nalgos*ninit_vals

  # Initilize the result dataframe
  cnames <- c("ID", "Algorithm", "Convergence", "Objective", "Iterations",
              "FuncEvals", "BmTime")
  res <- data.frame(matrix(ncol = length(cnames),
                           nrow = nids * ninit_vals * nalgos))
  colnames(res) <- cnames

  # label the result dataframe
  res[["ID"]] <- rep(ids, each = ni) # label individuals
  res[["Algorithm"]] <- rep(algo_names, each = ninit_vals) # label algorithms

  # Loop through each individual
  for (i in 1:nids) {
    # extract data for id i
    dfi <- df[df[["id"]] == ids[i], ]
    rti <- dfi[["rt"]]
    respi <- dfi[["response"]]
    truthi <- dfi[["truth"]]

    # starting value for t0 must be smaller than the smallest rt
    min_rti <- min(rti)
    t0_lo <- 0.01*min_rti
    t0_me <- 0.50*min_rti
    t0_hi <- 0.99*min_rti
    init_vals[["t0"]] <- c(rep(t0_me, 5), t0_lo, t0_hi, rep(t0_me, 4))

    # loop through all of the starting values
    for (j in 1:ninit_vals) {
      # get number of evaluations
      temp <- nlminb(init_vals[j, ], ll_ft_SWSE_17,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+0*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+0*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+0*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+0*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_ft_SWSE_14,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+1*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+1*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+1*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+1*ninit_vals+j] <- temp[["evaluations"]][[1]]
      temp <- nlminb(init_vals[j, ], ll_fb_SWSE_17,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+2*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+2*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+2*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+2*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fb_SWSE_14,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+3*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+3*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+3*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+3*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fb_Gon_17,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+4*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+4*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+4*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+4*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fb_Gon_14,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+5*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+5*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+5*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+5*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fb_Nav_17,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+6*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+6*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+6*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+6*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fb_Nav_14,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+7*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+7*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+7*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+7*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fs_SWSE_17,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+8*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+8*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+8*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+8*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fs_SWSE_14,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+9*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+9*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+9*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+9*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fs_Gon_17,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+10*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+10*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+10*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+10*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fs_Gon_14,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+11*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+11*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+11*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+11*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fs_Nav_17,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+12*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+12*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+12*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+12*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fs_Nav_14,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+13*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+13*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+13*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+13*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_fl_Nav_09,
                     rt = rti, resp = respi, truth = truthi, err_tol = err_tol,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+14*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+14*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+14*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+14*ninit_vals+j] <- temp[["evaluations"]][[1]]

      temp <- nlminb(init_vals[j, ], ll_RTDists,
                     rt = rti, resp = respi, truth = truthi,
                     # limits:   vu,   vl,   a,      t0, w,  sv
                     lower = c(-Inf, -Inf, .01,       0, 0,   0),
                     upper = c( Inf,  Inf, Inf, min_rti, 1, Inf),
                     control = ctrl_list)
      res[["Convergence"]][(i-1)*ni+15*ninit_vals+j] <- temp[["convergence"]]
      res[["Objective"]][(i-1)*ni+15*ninit_vals+j] <- temp[["objective"]]
      res[["Iterations"]][(i-1)*ni+15*ninit_vals+j] <- temp[["iterations"]]
      res[["FuncEvals"]][(i-1)*ni+15*ninit_vals+j] <- temp[["evaluations"]][[1]]

      # microbenchmark
      mbm <- microbenchmark(
        ft_SWSE_17 = nlminb(init_vals[j,], ll_ft_SWSE_17, err_tol = err_tol,
                            rt = rti, resp = respi, truth = truthi,
                            # limits:   vu,   vl,   a,      t0, w,  sv
                            lower = c(-Inf, -Inf, .01,       0, 0,   0),
                            upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        ft_SWSE_14 = nlminb(init_vals[j,], ll_ft_SWSE_14, err_tol = err_tol,
                            rt = rti, resp = respi, truth = truthi,
                            # limits:   vu,   vl,   a,      t0, w,  sv
                            lower = c(-Inf, -Inf, .01,       0, 0,   0),
                            upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fb_SWSE_17 = nlminb(init_vals[j,], ll_fb_SWSE_17, err_tol = err_tol,
                            rt = rti, resp = respi, truth = truthi,
                            # limits:   vu,   vl,   a,      t0, w,  sv
                            lower = c(-Inf, -Inf, .01,       0, 0,   0),
                            upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fb_SWSE_14 = nlminb(init_vals[j,], ll_fb_SWSE_14, err_tol = err_tol,
                            rt = rti, resp = respi, truth = truthi,
                            # limits:   vu,   vl,   a,      t0, w,  sv
                            lower = c(-Inf, -Inf, .01,       0, 0,   0),
                            upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fb_Gon_17 = nlminb(init_vals[j,], ll_fb_Gon_17, err_tol = err_tol,
                           rt = rti, resp = respi, truth = truthi,
                           # limits:   vu,   vl,   a,      t0, w,  sv
                           lower = c(-Inf, -Inf, .01,       0, 0,   0),
                           upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fb_Gon_14 = nlminb(init_vals[j,], ll_fb_Gon_14, err_tol = err_tol,
                           rt = rti, resp = respi, truth = truthi,
                           # limits:   vu,   vl,   a,      t0, w,  sv
                           lower = c(-Inf, -Inf, .01,       0, 0,   0),
                           upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fb_Nav_17 = nlminb(init_vals[j,], ll_fb_Nav_17, err_tol = err_tol,
                           rt = rti, resp = respi, truth = truthi,
                           # limits:   vu,   vl,   a,      t0, w,  sv
                           lower = c(-Inf, -Inf, .01,       0, 0,   0),
                           upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fb_Nav_14 = nlminb(init_vals[j,], ll_fb_Nav_14, err_tol = err_tol,
                           rt = rti, resp = respi, truth = truthi,
                           # limits:   vu,   vl,   a,      t0, w,  sv
                           lower = c(-Inf, -Inf, .01,       0, 0,   0),
                           upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fs_SWSE_17 = nlminb(init_vals[j,], ll_fs_SWSE_17, err_tol = err_tol,
                            rt = rti, resp = respi, truth = truthi,
                            # limits:   vu,   vl,   a,      t0, w,  sv
                            lower = c(-Inf, -Inf, .01,       0, 0,   0),
                            upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fs_SWSE_14 = nlminb(init_vals[j,], ll_fs_SWSE_14, err_tol = err_tol,
                            rt = rti, resp = respi, truth = truthi,
                            # limits:   vu,   vl,   a,      t0, w,  sv
                            lower = c(-Inf, -Inf, .01,       0, 0,   0),
                            upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fs_Gon_17 = nlminb(init_vals[j,], ll_fs_Gon_17, err_tol = err_tol,
                           rt = rti, resp = respi, truth = truthi,
                           # limits:   vu,   vl,   a,      t0, w,  sv
                           lower = c(-Inf, -Inf, .01,       0, 0,   0),
                           upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fs_Gon_14 = nlminb(init_vals[j,], ll_fs_Gon_14, err_tol = err_tol,
                           rt = rti, resp = respi, truth = truthi,
                           # limits:   vu,   vl,   a,      t0, w,  sv
                           lower = c(-Inf, -Inf, .01,       0, 0,   0),
                           upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fs_Nav_17 = nlminb(init_vals[j,], ll_fs_Nav_17, err_tol = err_tol,
                           rt = rti, resp = respi, truth = truthi,
                           # limits:   vu,   vl,   a,      t0, w,  sv
                           lower = c(-Inf, -Inf, .01,       0, 0,   0),
                           upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fs_Nav_14 = nlminb(init_vals[j,], ll_fs_Nav_14, err_tol = err_tol,
                           rt = rti, resp = respi, truth = truthi,
                           # limits:   vu,   vl,   a,      t0, w,  sv
                           lower = c(-Inf, -Inf, .01,       0, 0,   0),
                           upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        fl_Nav_09 = nlminb(init_vals[j,], ll_fl_Nav_09, err_tol = err_tol,
                           rt = rti, resp = respi, truth = truthi,
                           # limits:   vu,   vl,   a,      t0, w,  sv
                           lower = c(-Inf, -Inf, .01,       0, 0,   0),
                           upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        rtdists = nlminb(init_vals[j,], ll_RTDists,
                         rt = rti, resp = respi, truth = truthi,
                         # limits:   vu,   vl,   a,      t0, w,  sv
                         lower = c(-Inf, -Inf, .01,       0, 0,   0),
                         upper = c( Inf,  Inf, Inf, min_rti, 1, Inf)),
        times = times, unit = unit
      )
      for (k in 1:nalgos) {
        res[["BmTime"]][(i-1)*ni+(k-1)*ninit_vals+j] <- median(
          mbm[mbm[["expr"]] == algo_names[k], 2])
      }
    }
  }
  return(res)
}

Running the Fitting

As an example dataset, we use the fddm::med_dec data that comes with fddm. This dataset contains the accuracy condition reported in Trueblood et al. (2018), which investigates medical decision making among medical professionals (pathologists) and novices (i.e., undergraduate students). The task of participants was to judge whether pictures of blood cells show cancerous cells (i.e., blast cells) or non-cancerous cells (i.e., non-blast cells). The dataset contains 200 decision per participant, based on pictures of 100 true cancerous cells and pictures of 100 true non-cancerous cells. We remove the trials without response (indicated by rt \(< 0\) in the data) before fitting.

Having set up the fitting functions in the above chunks of code, we could pass the fddm::med_dec data to this function to get the parameter estimates. However, as this takes a relatively long time (around 30 minutes), we skip the fitting in this vignette and instead load the pre-fit parameter estimates in the next section.

data(med_dec, package = "fddm")
med_dec <- med_dec[which(med_dec[["rt"]] >= 0), ]
fit <- rt_fit(med_dec, id_idx = c(2,1), rt_idx = 8, response_idx = 7,
              truth_idx = 5, response_upper = "blast", err_tol = 1e-6,
              ctrl_list = list(eval.max = 10000, iter.max = 10000),
              times = 5, unit = "ns")

Analysis of Benchmark Results

We begin by loading the required packages for plotting the results.

library("reshape2")
library("ggplot2")

Before visualizing the results of the data fitting benchmark tests, we need to provide some context for each set of parameter estimates since each run of the optimization does not necessarily converge. To add this context, we first find the minimum value of the objective function (i.e., the log-likelihood function) attained for each set of initial parameter values for each individual. Then we calculate how the objective value of each optimization run compares to that minimum; we will refer to this quantity as the difference in minimized log-likelihood. The following code chunk defines a function to perform this calculation and some other minor operations that will be helpful when plotting the results. This code chunk also applies this prep function to the pre-run benchmark data and defines some constants that we will use throughout the visualizations.

fit_prep <- function(fit, eps = 1e-4) {
  nr <- nrow(fit)
  fit[["Obj_diff"]] <- rep(0, nr)

  ids <- unique(fit[["ID"]])
  nids <- length(ids)
  algos <- unique(fit[["Algorithm"]])
  nalgos <- length(algos)

  ninit <- nrow(fit[fit[["ID"]] == ids[1] & fit[["Algorithm"]] == algos[1], ])
  for (i in 1:nids) {
    for (j in 1:ninit) {
      idx <- which(fit[["ID"]] == ids[i])[ninit*(0:(nalgos-1)) + j]
      objs <- fit[idx, "Objective"]
      min_obj <- min(objs)
      abs_min_obj <- abs(min_obj)
      obj_diffs <- objs - min(objs)
      fit[idx, "Obj_diff"] <- ifelse(obj_diffs <= eps*abs_min_obj, 0,
        ifelse(obj_diffs > eps*abs_min_obj & obj_diffs <= 2*abs_min_obj, 1, 3))
    }
  }

  fit[["BmTime"]] <- fit[["BmTime"]]*1e-6 # convert to milliseconds
  fit[["Convergence"]] <- ifelse(fit[["Convergence"]] < 1, 0, 1)

  return(fit)
}

obj_diff_label <- function(y, df, col_name, mult = 1.15, upper_limit = NULL) {
  if (is.null(upper_limit)) {
    upper_limit <- max(df[[as.character(col_name)]])
  }
  return(
    data.frame(
      y = mult * upper_limit,
      label = paste(sum(y > 0, na.rm = TRUE))
    )
  )
}

Names <- c("ft_SWSE_17", "ft_SWSE_14", "fb_SWSE_17", "fb_SWSE_14",
           "fb_Gon_17", "fb_Gon_14", "fb_Nav_17", "fb_Nav_14",
           "fs_SWSE_17", "fs_SWSE_14", "fs_Gon_17", "fs_Gon_14",
           "fs_Nav_17", "fs_Nav_14", "fl_Nav_09", "rtdists")
Color <- c("#92c639", "#d3e8b0", "#92c639", "#d3e8b0",
           "#b3724d", "#e0c7b8", "#4da7b3", "#b8dce0",
           "#5cc639", "#bee8b0", "#b34d4d", "#e0b8b8",
           "#4d80b3", "#b8cce0", "#dcdca3", "#ac8053")
Outline <- c("#92c639", "#92c639", "#92c639", "#92c639",
             "#b3724d", "#b3724d", "#4da7b3", "#4da7b3",
             "#5cc639", "#5cc639", "#b34d4d", "#b34d4d",
             "#4d80b3", "#4d80b3", "#dcdca3", "#ac8053")
Shape <- c(21, 25)
Sizes <- c(0, 3, 3)
Stroke <- c(0, 1, 1)
Fills <- c("#ffffff00", "#ffffff00", "#80808099")

# load data, will be in the variable 'fit'
load(system.file("extdata", "dfddm_density", "bm_fit.Rds",
                 package = "fddm", mustWork = TRUE))
fit <- fit_prep(fit)

We will now present two plots to illustrate how the various implementations perform when used in fitting Ratcliff’s DDM to real-world data. The measures of interest for each implementation are the microbenchmark::microbenchmark() timings and the number of function evaluations in the optimization process. Each one of the two plots will focus on one of these measures by showing a distribution of values of the measure with violin plots and box plots. The distribution shown is the distribution across participants and starting values for each of the different implementations of approximating the diffusion model; that is, for each of the 37 participants, we present results based on the 11 different starting values for each implementation (407 data points in total).

To help explain problematic fits, the plots will show any data point that represents an optimization run for which the difference in minimized log-likelihood is greater than our predefined tolerance of \(\epsilon = 0.0001\). More specifically, the difference in minimized log-likelihood is comparatively small (i.e., smaller than 2 on the log-likelihood scale) for the points shown with a completely transparent fill color; contrastingly, points shown with a gray fill color have a difference from the optimum of larger than 2 points on the log-likelihood scale (2 points on the log-likelihood scale corresponds to 1 point on the deviance scale). Data points with difference smaller than \(\epsilon\) will not be shown on the plot. In addition, each of these visible data points will have a shape to indicate its convergence code: a circle for successful convergence, and a triangle for failed convergence.

The first plot will show the densities of the microbenchmark::microbenchmark() timings for each implementation.

fit_mbm <- melt(fit, id.vars = c("Algorithm", "Convergence", "Obj_diff"),
                measure.vars = "BmTime", value.name = "BmTime")[,-4]

mi <- min(fit[fit[["Algorithm"]] != "rtdists", "BmTime"])
ma <- max(fit[fit[["Algorithm"]] != "rtdists", "BmTime"])

ggplot(fit_mbm, aes(x = factor(Algorithm, levels = Names),
                    y = BmTime)) +
  geom_violin(trim = TRUE, alpha = 0.5,
              aes(color = factor(Algorithm, levels = Names),
                  fill = factor(Algorithm, levels = Names))) +
  geom_boxplot(width = 0.2, outlier.shape = NA,
               fill = "white", alpha = 0.4,
               aes(color = factor(Algorithm, levels = Names))) +
  stat_summary(fun = mean, geom = "errorbar",
               aes(ymax = after_stat(y), ymin = after_stat(y)),
               width = .5, linetype = "dashed",
               color = Color) +
  stat_summary(aes(y = Obj_diff, color = factor(Algorithm, levels = Names)),
               fun.data = obj_diff_label,
               fun.args = list(fit, "BmTime", 1.075, ma),
               geom = "label",
               hjust = 0.5,
               vjust = 0.9) +
  scale_x_discrete(labels = c(
    bquote(f[t] ~ SWSE[17]), bquote(f[t] ~ SWSE[14]),
    bquote(f[c] ~ SWSE[17]), bquote(f[c] ~ SWSE[14]),
    bquote(f[c] ~ Gon[17]), bquote(f[c] ~ Gon[14]),
    bquote(f[c] ~ Nav[17]), bquote(f[c] ~ Nav[14]),
    bquote(f[s] ~ SWSE[17]), bquote(f[s] ~ SWSE[14]),
    bquote(f[s] ~ Gon[17]), bquote(f[s] ~ Gon[14]),
    bquote(f[s] ~ Nav[17]), bquote(f[s] ~ Nav[14]),
    bquote(f[l] ~ "Nav"), "rtdists")) +
  coord_cartesian(ylim = c(mi, ma*1.05)) +
  scale_color_manual(values = Outline, guide = "none") +
  scale_fill_manual(values = Color, guide = "none") +
  scale_shape_manual(values = Shape,
                     name = "Convergence Code",
                     breaks = c(0, 1),
                     labels = c("Success", "Failure")) +
  scale_size_manual(values = Sizes, guide = "none") +
  scale_discrete_manual(aesthetics = "stroke", values = Stroke, guide = "none") +
  ggnewscale::new_scale_fill() +
  scale_fill_manual(values = Fills,
                    name = paste("Difference in", "Log-likelihood", "from MLE",
                                 sep = "\n"),
                    breaks = c(1, 2, 3),
                    labels = c("< 2", "NA", "> 2")) +
  geom_point(aes(color = factor(Algorithm, levels = Names),
                 shape = factor(Convergence, levels = c(0, 1)),
                 size = factor(Obj_diff, levels = c(0, 1, 3)),
                 stroke = factor(Obj_diff, levels = c(0, 1, 3)),
                 fill = factor(Obj_diff, levels = c(0, 1, 3)))) +
  labs(x = "Implementation", y = "Time (milliseconds)") +
  guides(shape = guide_legend(order = 1,
                              override.aes = list(size = Sizes[c(2, 3)])),
         fill = guide_legend(order = 2,
                             override.aes = list(size = Sizes[c(2, 3)],
                                                 shape = c(21, 21),
                                                 fill = Fills[c(2, 3)]))) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 16, angle = 90,
                                   vjust = 0.5, hjust = 1),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20,
                                    margin = margin(15, 0, 0, 0)),
        axis.title.y = element_text(size = 20,
                                    margin = margin(0, 10, 0, 0)),
        legend.position = "right",
        legend.box = "vertical",
        legend.direction = "vertical",
        legend.background = element_rect(fill = "transparent"),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 13))

This plot vaguely resembles the microbenchmark violin plot from the previous section. Immediately we observe that fitting to data using the rtdists function is noticeably slower than fitting to data using the functions provided in fddm, but the magnitude of the difference is much smaller. The large-time implementation has many problematic data points, further emphasizing its shortcomings when handling small effective response times.

Upon closer inspection, the means and medians of the benchmark times are slightly lower for the small-time implementations than for the implementations that combine both the small-time and large-time density function approximations. However, these two seemingly faster small-time implementations have more problematic data points in the optimization process, due to the poor support for large effective response times that we discussed in a previous section of this vignette.

We also see that the two combined-time implementations only show relatively very few cases where the optimization algorithm did not reach the optimum for one participant and starting value. These issues can be avoided by using multiple sets of initial values and choosing the best optimization result from that set. Despite these few problematic fits, the combined-time implementations are very stable and fast.

If the data is known to contain only small response times, then using a small-time implementation may yield fast and unproblematic fitting results. However, instability may still arise from the optimization routine testing small values of \(a\), as this would create large effective response times and cause the small-time aproximations to reveal their vulnerability. For this reason, we recommend using an implementation that combines the small-time and large-time approximations instead of an implementation that uses only one timescale.

The second plot will show the distribution of the number of function evaluations in the optimization process for each implementation.

fit_fev <- melt(fit, id.vars = c("Algorithm", "Convergence", "Obj_diff"),
                measure.vars = "FuncEvals", value.name = "FuncEvals")[,-4]

ggplot(fit_fev, aes(x = factor(Algorithm, levels = Names),
                              y = FuncEvals)) +
  geom_violin(trim = TRUE, alpha = 0.5,
              aes(color = factor(Algorithm, levels = Names),
                  fill = factor(Algorithm, levels = Names))) +
  geom_boxplot(width = 0.2, outlier.shape = NA,
               fill = "white", alpha = 0.4,
               aes(color = factor(Algorithm, levels = Names))) +
  stat_summary(fun = mean, geom = "errorbar",
               aes(ymax = after_stat(y), ymin = after_stat(y)),
               width = .5, linetype = "dashed",
               color = Color) +
  stat_summary(aes(y = Obj_diff, color = factor(Algorithm, levels = Names)),
               fun.data = obj_diff_label,
               fun.args = list(fit, "FuncEvals", 1.075),
               geom = "label",
               hjust = 0.5,
               vjust = 0.9) +
  scale_x_discrete(labels = c(
    bquote(f[t] ~ SWSE[17]), bquote(f[t] ~ SWSE[14]),
    bquote(f[c] ~ SWSE[17]), bquote(f[c] ~ SWSE[14]),
    bquote(f[c] ~ Gon[17]), bquote(f[c] ~ Gon[14]),
    bquote(f[c] ~ Nav[17]), bquote(f[c] ~ Nav[14]),
    bquote(f[s] ~ SWSE[17]), bquote(f[s] ~ SWSE[14]),
    bquote(f[s] ~ Gon[17]), bquote(f[s] ~ Gon[14]),
    bquote(f[s] ~ Nav[17]), bquote(f[s] ~ Nav[14]),
    bquote(f[l] ~ "Nav"), "rtdists")) +
  scale_color_manual(values = Outline, guide = "none") +
  scale_fill_manual(values = Color, guide = "none") +
  scale_shape_manual(values = Shape,
                     name = "Convergence Code",
                     breaks = c(0, 1),
                     labels = c("Success", "Failure")) +
  scale_size_manual(values = Sizes, guide = "none") +
  scale_discrete_manual(aesthetics = "stroke", values = Stroke, guide = "none") +
  ggnewscale::new_scale_fill() +
  scale_fill_manual(values = Fills,
                    name = paste("Difference in", "Log-likelihood", "from MLE",
                                 sep = "\n"),
                    breaks = c(1, 2, 3),
                    labels = c("< 2", "NA", "> 2")) +
  geom_point(aes(color = factor(Algorithm, levels = Names),
                 shape = factor(Convergence, levels = c(0, 1)),
                 size = factor(Obj_diff, levels = c(0, 1, 3)),
                 stroke = factor(Obj_diff, levels = c(0, 1, 3)),
                 fill = factor(Obj_diff, levels = c(0, 1, 3)))) +
  labs(x = "Implementation", y = "Number of function evaluations") +
  guides(shape = guide_legend(order = 1,
                              override.aes = list(size = Sizes[c(2, 3)])),
         fill = guide_legend(order = 2,
                             override.aes = list(size = Sizes[c(2, 3)],
                                                 shape = c(21, 21),
                                                 fill = Fills[c(2, 3)]))) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text.x = element_text(size = 16, angle = 90,
                                   vjust = 0.5, hjust = 1),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20,
                                    margin = margin(15, 0, 0, 0)),
        axis.title.y = element_text(size = 20,
                                    margin = margin(0, 10, 0, 0)),
        legend.position = "right",
        legend.box = "vertical",
        legend.direction = "vertical",
        legend.background = element_rect(fill = "transparent"),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 13))

Here we see that the implementations that use both the small-time and large-time density function approximations appear to be the most stable. They have relatively short tails, suggesting that the optimization routine has little trouble in finding the optimum parameter estimates. Moreover, these implementations have few data points where the objective difference is high or the convergence code indicates unsuccessful convergence. In contrast, the large-time implementation suffers greatly from its inability to handle small effective response times. The rtdists function performs well in this test relative to the previous tests, but it is still hindered by longer benchmark times.

In addition to being the fastest implementation available, the default implementation for fddm, “\(f_t \text{SWSE}_{17}\)”, appears to be the most stable implementation for fitting the DDM to real-world data. The optimization routine handles this implementation consistently well in that there are very few troublesome fits compared to that of the implementations which only implement the small-time density function approximation. The default implementation also outperforms its counterparts, “\(f_c \text{Gon}_{17}\)” and “\(f_c \text{Nav}_{17}\)”, because it uses a faster mechanism to choose between the two timescales and exploits the large-time implementation where it is most useful.

When comparing the implementations for benchmark time, the small-time implementations have a faster mean and median than the default implementation, but their lengthy tails indicate potential stability issues during data fitting. On balance, the default implementation, “\(f_t \text{SWSE}_{17}\)”, provides the best compromise between speed and stability of the currently available implementations for approximating the probability density functions of the diffusion decision model. For more information about the default implementation, see its section in the Math Vignette.

Benchmark Testing for the DDM Density Function

Kendal Foster and Henrik Singmann

July 01, 2024

Introduction

Benchmarking the Density Function Approximations

Generating Benchmark Data

Analysis of Benchmark Results

Benchmarking Model Fitting to Real-World Data

Generating Benchmark Data for Parameter Estimation

Log-Likelihood Functions

Fitting Function

Running the Fitting

Analysis of Benchmark Results

R Session Info

References