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:
The package can be installed from CRAN directly or through GitHub using these function calls:
Load the package as you typically would:
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"
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
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.
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.
Depending on your needs, you can choose between three cases:
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:
models_t.h
-> class
CSTM_T
models_tx.h
-> class
CSTM_TX
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).
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
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).
In the Terminal (of RStudio) navigate to the folder on top of
ream
and execute:
R CMD build rtmpt
R CMD check rtmpt_<version_number>.tar.gz --as-cran
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:
Build>Load All
(Ctrl+Shift+L)Build>Check Package
(Ctrl+Shift+E)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.
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/
).