ream: guideline

Raphael Hartmann and Matthew Murrow

2024-09-25

R Package to Calculate the PDF and CDF of Evidence Accumulation Models

Overview

In the package the following functions are implemented:

All these functions can be used for the following evidence accumulation models (EAM):

Model name Short model name
Simple DDM SDDM
Leaky integration model LM
Linear threshold model LTM
Exponential threshold model ETM
Weibull threshold model WTM
Urgency gating model UGM
Diffusion model of conflict DMC
Revised diffusion model of conflict RDMC
Shrinking spotlight model SSP
Dual process model DPM
Dual-stage two-process model DSTP

In this guideline the following use cases are discussed:

  1. Sampling from a specific EAM
  2. Calculating the PDF and CDF of a specific EAM (and plotting the PDF)
  3. Fitting a specific EAM to data
  4. Implement your custom EAM

0 Preparation

The package can be installed from CRAN directly or through GitHub using these function calls:

devtools::install_github("RaphaelHartmann/ream")  # installing from GitHub

Load the package as you typically would:

library(ream)

1 Random Sampling

For taking random samples (first-passage/response times and thresholds/responses) you can use the same logic as for R-native random sampling functions like rnorm(); Just prepend an r in front of the short model name (see table above; e.g., for the DMC it would be rDMC()). The random sampling function has three arguments: n the number of random samples, phi the vector of model parameters in a specific order which can be found on its help file, and dt the step size of the time in the trajectory.

For example, sampling from the diffusion model of conflict (DMC) one would call the function rDMC() like this:

?rDMC # check the help file
(samp <- rDMC(n = 10, phi = c(0.3, 0.5, -1, 0.2, 0.05, 2.5, 3, 1, 0.5, 0, 0, 1), dt = 1e-4))
## $rt
##  [1] 0.3444762 0.6515139 0.5390954 0.5391878 0.4224804 0.5194858 0.4760704
##  [8] 0.4459530 0.6134381 0.4717932
## 
## $resp
##  [1] "upper" "upper" "upper" "upper" "upper" "upper" "upper" "upper" "upper"
## [10] "upper"

2 Calculating the PDF and CDF

For calculating the PDF and CDF you can also use the same logic as for native PDFs and CDFs like dnorm() and pnorm(), respectively; Just prepend a d or p, respectively, in front of the short model name (see table above; e.g., for the DMC it would be dDMC() and pDMC(), respectively). These functions have five arguments: rt the response times (in seconds) as a vector, resp the responses ("upper" and "lower"), phi the vector of model parameters in a specific order which can be found on its help file, x_res the spatial/evidence resolution ("A", "B", "C", or "D"), and t_res the time resolution ("A", "B", "C", or "D").

For example, calculating the PDF and CDF of the above samples samp for the diffusion model of conflict (DMC) would look like this:

?dDMC # check the help file
(PDF <- dDMC(rt = samp$rt, resp = samp$resp, phi = c(0.3, 0.5, -1, 0.2, 0.05, 2.5, 3, 1, 0.5, 0, 0, 1), x_res = "default", t_res = "default"))
## $sum_log_pdf
## [1] 7.279195
## 
## $pdf
##  [1] 0.9479148 4.0056953 1.9278749 0.4998694 1.9300591 2.4637257 4.9919678
##  [8] 4.1669263 0.7742009 5.1740108
## 
## $log_pdf
##  [1] -0.0534907  1.3877172  0.6564183 -0.6934084  0.6575506  0.9016747
##  [7]  1.6078302  1.4271787 -0.2559239  1.6436482
(CDF <- dDMC(rt = samp$rt, resp = samp$resp, phi = c(0.3, 0.5, -1, 0.2, 0.05, 2.5, 3, 1, 0.5, 0, 0, 1), x_res = "default", t_res = "default"))
## $sum_log_pdf
## [1] 7.278755
## 
## $pdf
##  [1] 0.9478176 4.0048859 1.9278608 0.5003681 1.9301252 2.4628866 4.9899620
##  [8] 4.1672255 0.7739269 5.1733127
## 
## $log_pdf
##  [1] -0.05359319  1.38751509  0.65641099 -0.69241123  0.65758489  0.90133408
##  [7]  1.60742829  1.42725046 -0.25627789  1.64351323

The functions calculate the PDF (or CDF), their log values, and the sum of the log values all together and save them in a list.

GRID

3 Fitting a Model to Data

For fitting an EAM to data one can use the optim() function with the "L-BFGS-B" method (or nlminb() function). The "L-BFGS-B" method allows for defining lower and upper bounds for the parameters to be fitted.

Before we fit a model, we generate some data. For simplicity we simulate only one subject with 100 congruent and 100 incongruent trials.

N <- 100
con <- rDMC(n = N, phi = c(0.3, 0.5, 1, 0.2, 0.05, 2.5, 3, 1, 0.5, 0, 0, 1), dt = 1e-05)
incon <- rDMC(n = N, phi = c(0.3, 0.5, -1, 0.2, 0.05, 2.5, 3, 1, 0.5, 0, 0, 1), dt = 1e-05)
data <- data.frame(congruency = rep(1:2, each = N), rt = c(con$rt, incon$rt), resp = c(con$resp, incon$resp))

We now want to fit the DMC to the data. We start with defining the function to maximize.

deviation <- function(pars, data) {
  ind_con <- which(data$congruency==1)
  ind_incon <- which(data$congruency==2)
  ls_con <- dDMC(rt = data$rt[ind_con], resp = data$resp[ind_con], phi = c(pars[1:2], 1, pars[3:6], 1, 0.5, 0, 0, 1), x_res = "higher", t_res = "higher")$sum_log_pdf
  ls_incon <- dDMC(rt = data$rt[ind_incon], resp = data$resp[ind_incon], phi = c(pars[1:2], -1, pars[3:6], 1, 0.5, 0, 0, 1), x_res = "higher", t_res = "higher")$sum_log_pdf
  return(-2*(ls_con+ls_incon))
}

Now we only need to call the optim() function.

set.seed(3210)
(start_pars <- c(runif(1, .2, .6), runif(1, .3, .7), runif(1, .1, .6), runif(1, 0, .1), runif(1, 1.5, 4.5), runif(1, 0, 5)))
## [1] 0.3223039 0.4210046 0.1069716 0.0691608 4.0316741 3.2418832
optim(par = start_pars, fn = deviation, method = "L-BFGS-B", lower = c(.2, .1, .1, 0.001, 1, 0.001), upper = c(.6, .9, .6, .1, 5, 5), data = data)
## $par
## [1] 0.28784300 0.43358114 0.27300557 0.03695141 4.04280133 3.26391350
## 
## $value
## [1] -430.0266
## 
## $counts
## function gradient 
##       28       28 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Keep in mind, that more data with more experimental conditions might be needed for better estimations.

4. Implement your custom EAM

4.1 Get the source code

Download from GitHub or from CRAN the source code and extract it.

4.2 Get your system ready

Every operating has some special requirements for building R packages. You can find them, for example, in the book R Packages by Wickham and Bryan.

4.3 What source files to modify

Depending on your needs, you can choose between three cases:

  • Case 1: drift rate function only dependent on time t, but not on evidence state x, or dependent on neither of them
  • Case 2: drift rate function dependent on time t and evidence state x
  • Case 3: drift rate and diffusion rate function dependent on time t and some weighting w, but not on evidence state x

The last case is rather a rare case. If you do not want the drift rate to have any dependency on time or evidence state, you pick the first case. If your drift rate does depend on time but not on evidence, choose also the first case. And so on.

Depending on the case at hand, different source files should be modified in the source code under ream/src/. Modify the following class methods depending on the case at hand:

  • Case 1: In models_t.h -> class CSTM_T
  • Case 2: In models_tx.h -> class CSTM_TX
  • Case 3: In models_tw.h -> class CSTM_TW

These are already prepared classes for custom models and have corresponding R functions (e.g., dCSTM_T()). The only thing to do is to change the methods of the classes and to recompile the R package (probably with a different version number).

4.4 How to modify the necessary class methods

Do not change the arguments!

Only change the content of the function. For example, if you want to change the threshold to be dependent on time in an exponential way (like the ETM) change the following method (e.g., inside the CSTM_T class) from

double upper_threshold(const double phi[100], double t) const override {
  return phi[4];
}

to

double upper_threshold(const double phi[100], double t) const override {
  double b = phi[4];
  double tau = pow(10.0, phi[5]);
  double thres = b*exp(-t/tau);
  return thres;
}

In this example, one would also change the lower threshold to:

double lower_threshold(const double phi[100], double t) const override {
  return -upper_threshold(phi, t);
}

Make sure to get the indexing in phi[i] correct. In our example, phi[5] is already used by the next function contamination_strength(), and should therefore be changed. The easiest way is to go from the first method to the last method in the class with increasing index of phi[i] (see the other classes).

4.5 Recompile the R package

In the Terminal (of RStudio) navigate to the folder on top of ream and execute:

  1. R CMD build rtmpt
  2. R CMD check rtmpt_<version_number>.tar.gz --as-cran
  3. R CMD INSTALL rtmpt_<version_number>.tar.gz

where <version_number> is the version number from the downloaded package or the one you changed the package to in DESCRIPTION.

Or, if you have set up RStudio correctly and made an R-project of the package you can also click the following instead of using the terminal:

  1. Build>Load All (Ctrl+Shift+L)
  2. Build>Check Package (Ctrl+Shift+E)
  3. Build>Install Package (Ctrl+Shift+B)

You should now be ready to use the R functions CSTM_T(), CSTM_TX(), and/or CSTM_TW() with your custom model(s) once you have loaded the correct version of the package.

4.6 Some notes

If you have two or more custom models of the same case, we suggest copying the source package multiple times, make each copy a different version, and do the above steps for each copy. Note that the package/project folder should always be called ream. Therefore, copy it to different top-level folders (e.g., .../REAM1/ream/ and .../REAM2/ream/).

References