dRiftDM: Estimating (Time-Dependent) Drift Diffusion Models

Description

Fit and explore Drift Diffusion Models (DDMs), a common tool in psychology for describing decision processes in simple tasks. It can handle both time-independent and time-dependent DDMs. You either choose prebuilt models or create your own, and the package takes care of model predictions and parameter estimation. Model predictions are derived via the numerical solutions provided by Richter, Ulrich, and Janczyk (2023, doi:10.1016/j.jmp.2023.102756).

Author(s)

Maintainer: Valentin Koob v.koob@web.de [copyright holder]

Authors:

• Thomas Richter [copyright holder]

• Markus Janczyk

See Also

Useful links:

• https://github.com/bucky2177/dRiftDM

• https://bucky2177.github.io/dRiftDM/

• Report bugs at https://github.com/bucky2177/dRiftDM/issues

Convolute the First Passage Times with the Non-Decision Time Distribution

Description

Calls stats::convolve for the first passage times and the non-decision time distribution to derive the full distribution of response times. Before convolution, I add the robustness parameter.

Usage

add_residual(pdf_nt, pdf_u, pdf_l, dt, nt)

Arguments

Value

a list of PDFs for one condition "pdf_u" and "pdf_l"

Aggregate Statistics Across ID

Description

aggregate_stats is a (not exported) function to aggregate stats_dm objects across IDs. Since the column names may vary by the statistic type, the behavior of aggregate depends on the subclass of stats_dm.

Usage

aggregate_stats(stat_df)

Arguments

Details

aggregate_stats calls the internal_aggregate() with the relevant arguments

Value

Returns the statistics aggregated across the relevant cols.

See Also

new_stats_dm, calc_stats, internal_aggregate()

The Coding of the Boundaries

Description

Functions to get or set the "boundary coding" of an object.

Usage

b_coding(object, ...) <- value

## S3 replacement method for class 'drift_dm'
b_coding(object, ...) <- value

b_coding(object, ...)

## S3 method for class 'drift_dm'
b_coding(object, ...)

## S3 method for class 'fits_ids_dm'
b_coding(object, ...)

## S3 method for class 'fits_agg_dm'
b_coding(object, ...)

Arguments

Details

b_coding() is a generic accessor function, and ⁠b_coding<-()⁠ a generic replacement function. The default methods get and set the "boundary coding", which is an attribute of drift_dm model.

The boundary coding summarizes which response time belongs to which boundary and how the boundaries shall be "labeled". The list specifies three entries:

• column, contains a single character string, indicating which column in an observed data set codes the boundaries.

• u_name_value, contains a numeric or character vector of length 1. The name of this vector gives a label for the upper boundary, and the entry gives the value stored in obs_data[[column]] coding the upper boundary.

• l_name_value, contains a numeric or character vector of length 1. The name of this vector gives a label for the lower boundary, and the entry gives the value stored in obs_data[[column]] coding the lower boundary.

The package dRiftDM has a default boundary coding:

• column = "Error"

• u_name_value = c("corr" = 0)

• l_name_value = c("err" = 1)

Thus, per default, dRiftDM assumes that any observed data set has a column "Error", providing the values 0 and 1 for the upper and lower boundary, respectively. The upper and lower boundaries are labeled "corr" and "err", respectively. These labels are used, for example, when calculating statistics (see calc_stats).

When calling ⁠b_coding<-()⁠ with value = NULL, the default "accuracy" coding is evoked

Value

For b_coding() a list containing the boundary coding For ⁠b_coding<-()⁠ the updated drift_dm or fits_ids_dm object

See Also

drift_dm()

Examples

# show the default accuracy coding of dRiftDM
my_model <- ratcliff_dm() # get a pre-built model
b_coding(my_model)

# can be modified/replaced
b_coding(my_model)[["column"]] <- "Response"

# accessor method also available for fits_ids_dm objects
# get an exemplary fits_ids_dm object (see estimate_model_ids)
fits <- get_example_fits("fits_ids_dm")
names(b_coding(fits))

Constant Boundary

Description

Constant Boundary

Usage

b_constant(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

a vector of the same length as t_vec with the b value for each element of the vector.

Collapsing Boundary - Hyperbolic Ratio Function

Description

Provides the boundary, collapsing as a hyperbolic ratio function.

Usage

b_hyperbol(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

dt_b_hyperbol(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Details

b_hyperbol and dt_b_hyperbol provide the plain boundary values and the respective derivative, respectively.

Value

a vector of the same length as t_vec with the boundary values (or the deriviative) for each element of the vector.

Collapsing Boundary - Weibull Function

Description

Provides the boundary, collapsing in accordance with a Weibull function

Usage

b_weibull(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

dt_b_weibull(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Details

b_weibull and dt_b_weibull provide the plain boundary values and the respective derivative, respectively.

Value

a vector of the same length as t_vec with the boundary values (or the deriviative) for each element of the vector.

Calculate Basic Statistics for Response Times or Probability Densities

Description

Function that calls the underlying functions calc_basic_stats_obs and calc_basic_stats_pred. Handles input checks and data wrangling.

Usage

calc_basic_stats(
  pdf_u = NULL,
  pdf_l = NULL,
  rts_u = NULL,
  rts_l = NULL,
  one_cond,
  b_coding,
  t_vec = NULL,
  dt = NULL,
  skip_if_contr_low = NULL
)

Arguments

Details

• If pdf_u and pdf_l are provided, the function computes statistics for the probability densities.

• If rts_u and rts_l are provided, the function computes statistics for the observed RTs.

• If both sets of inputs are provided, both types of statistics are computed and returned.

Value

A data.frame with columns:

• Source: Indicates whether the statistics refer to observed (obs) or predicted (pred) data.

• Cond: A condition label.

• ⁠Mean_<u_label>⁠: Mean response time for the upper boundary.

• ⁠Mean_<l_label>⁠: Mean response time for the lower boundary.

• ⁠P_<u_label>⁠: Proportion of upper-boundary responses.

See Also

calc_basic_stats_obs, calc_basic_stats_pred, new_stats_dm

Calculate Basic Mean, Standard Deviations, and Percentages for Observed and Predicted Data

Description

Backend functions to compute the mean response times, standard deviations of response times, and response proportions; for both observed RTs or the predicted probability density functions.

Usage

calc_basic_stats_obs(rts_u, rts_l, one_cond)

calc_basic_stats_pred(
  pdf_u,
  pdf_l,
  one_cond,
  t_vec,
  dt,
  skip_if_contr_low = NULL
)

Arguments

Details

• For observed data, calculates mean RTs, standard deviations of RTs, and the proportion of upper responses.

• The same statistics are calculated for the probability density values (via simple numerical integration)

Value

A data.frame with columns:

• Cond: Condition label.

• Mean_U: Mean response time for the upper boundary

• Mean_L: Mean response time for the lower boundary

• SD_U: Mean response time for the upper boundary

• SD_L: Mean response time for the lower boundary

• P_U: Proportion of upper-boundary responses.

Calculate CAFs

Description

Function that calls the underlying CAF calculation functions calc_cafs_obs and calc_cafs_pred. Does input checks and the data wrangling

Usage

calc_cafs(
  pdf_u = NULL,
  pdf_l = NULL,
  t_vec = NULL,
  rts_u = NULL,
  rts_l = NULL,
  one_cond,
  n_bins = NULL,
  b_coding
)

Arguments

Details

if pdf_u and pdf_l are not NULL, returns CAFs of the densities

if rts_u and rts_l are not NULL, returns CAFs of the response times

if all are not NULL, returns both.

Value

a data.frame with "Source", "Cond", "Bin"s, "P_<u_label>" for the CAFs of type c("cafs", "sum_dist", "stats_dm", "data.frame")

See Also

new_stats_dm()

Calculate CAFs

Description

Backend functions to calculate conditional accuracy functions for RT vectors or pdfs

Usage

calc_cafs_obs(rts_u, rts_l, one_cond, n_bins)

calc_cafs_pred(pdf_u, pdf_l, t_vec, one_cond, n_bins)

Arguments

Details

for RTs: first elements are attributed to a bin (with bins calculated across all RTs using equally spaced quantiles), then accuracy per bin is calculated.

for Densities: Add density values, calculate a CDF and force it between 0 and

1. Then determine the indices that cut the CDF into bins by considering equally spaced quantiles. Then calculate the ratio of probability mass per bin.

Value

a data.frame with the "Cond" label, the "Bin"s and "P_U" for the CAFs

Calculate delta function(s)

Description

Given a dataset providing the quantiles (calc_quantiles), calculates delta function(s) for the character vectors minuends and subtrahends

Usage

calc_delta_funs(
  quantiles_dat,
  minuends = NULL,
  subtrahends = NULL,
  dvs = NULL,
  b_coding
)

Arguments

Details

Takes the quantile data_frame, stats::reshape it to wide, and then access the relevant dv columns, together with minuends and subtrahends to calculate the delta functions.

Value

a data.frame with columns "Source", "Prob", the "Quant_<u_label>", "Quant_<l_label". May have the following additional columns:

• if only one dv: as many Delta_<minuend_subtrahend> and Avg_<minuends_subtrahends> as minuends and subtrahends.

• if more than one dv: as many Delta_<u/l-label><minuend_subtrahend> and Avg<u/l-label>_<minuends_subtrahends> as minuends and subtrahends.

The data.frame is of type c("delta_funs", "sum_dist", "stats_dm", "data.frame")

See Also

new_stats_dm()

Calculate and Combine Density Estimates for Observed and Predicted Data

Description

Internal helper functions that return a data.frame summarizing density values for observed and predicted response times.

Usage

calc_dens_obs(
  rts_u,
  rts_l,
  one_cond,
  t_max = NULL,
  discr = NULL,
  scaling_factor = 1
)

calc_dens(
  pdf_u = NULL,
  pdf_l = NULL,
  t_vec = NULL,
  t_max = NULL,
  discr = NULL,
  rts_u = NULL,
  rts_l = NULL,
  one_cond,
  b_coding,
  scaling_factor = 1
)

Arguments

Details

calc_dens_obs() computes empirical histograms and kernel density estimates for a single condition based on observed RTs.

calc_dens() serves as a general interface that combines observed and predicted data into a single data.frame. Observed data (rts_u and rts_l) is passed to calc_dens_obs(). Predicted data (pdf_u and pdf_l) is wrapped into a data.frame that matches the structure returned by calc_dens_obs(). If both are provided, observed and predicted data are row-bound into a single data.frame.

These functions are used internally to support type = "density" in calc_stats(), providing a full distributional overview.

Value

A stats_dm object (via new_stats_dm()) containing a data.frame with columns:

• Source: indicates whether the row is from observed ("obs") or predicted ("pred") data.

• Cond: the condition label.

• Time: the time point corresponding to the density value.

• Stat: type of density summary—"hist" or "kde" (for observed data), or "pdf" (for predicted data).

• ⁠Dens_<U>⁠: density value for the upper response.

• ⁠Dens_<L>⁠: density value for the lower response.

The ⁠<U>⁠ and ⁠<L>⁠ placeholders are determined by the b_coding argument.

Calculate Fit Statistics

Description

Computes/Summarizes multiple fit statistics, inclduing Log-Likelihood, the Negative Log-Likelihood, the Akaike Information Criterion (AIC), the Bayesian Information Criterion (BIC), and the Root-Mean Squared-Error (RMSE) statistic.

Usage

calc_fit_stats(drift_dm_obj, k = 2, ...)

Arguments

Value

A custom object of class stats_dm (c("fit_stats", "stats_dm", "data.frame")). The columns are:

• Log_Like: the log-likelihood value

• Neg_Log_Like: the negative log-likelihood value

• AIC: the calculated AIC value

• BIC: the calculated BIC value

• RMSE_s: the root-mean-squared error (for RTs in seconds)

• RMSE_ms: the root-mean-squared error (for RTs in milliseconds) If a respective statistic cannot be calculated, the respective column contains NA.

See Also

new_stats_dm(), logLik.drift_dm

Calculate the Log-Likelihood

Description

Wrapper function around log_like_heart

Usage

calc_log_like(pdfs, t_vec, obs_data)

log_like_heart(pdf_u, pdf_l, t_vec, rts_u, rts_l)

Arguments

Details

calc_log_like

Iterates over all conditions, and passes forward the (unpacked) arguments to log_like_heart, adding each log-likelihood of a condition.

pdfs must be a list with entries named as the conditions, and then each condition being a list of the two PDFs (named pdf_u and pdf_l)

obs_data must be a list with entries "rts_u" and "rts_l", and then each rts_* entry being a named list with the RT values for each condition

log_like_heart

Gets the density values for RTs in rts_u/rts_l via stats::approx(), takes the log of that, and then sums across both. Wraps up the calculation in a tryCatch statement, throwing warnings when log_like_values can not be calculated

Value

a single value of the log-likelihood. If no data is provided, NULL is returned. If the calculation fails, -Inf is returned.

Calculate the PDFs

Description

This method takes a model, the time and space vectors, and the unpacked parameters for solving the PDF to derive the first passage time across all conditions. It is a wrapper around the cpp implementations and add_residual. Important: This function is used in the depths of the package and the generic method is not exported.

Usage

calc_pdfs(drift_dm_obj, x_vec, t_vec, prms_solve)

## S3 method for class 'drift_dm'
calc_pdfs(drift_dm_obj, x_vec, t_vec, prms_solve)

Arguments

Details

calc_pdfs is a generic method which dispatches the function call (not exported). Currently, the method only considers objects of type drift_dm.

calc_pdfs.drift_dm is the function that will be called for all models. It evaluates the different components of a model, and subsequently calls the cpp implementations for the KFE or integral method. It also calls the add_residual function to convolute the non-decision time to the first passage time.

The numerical methods for deriving the PDFs are based on the code provided by (Richter et al. 2023).

Value

a list of PDFs, with named entries for each condition. Each of this entry contains a list of vectors, named "pdf_u" and "pdf_l"

References

Richter T, Ulrich R, Janczyk M (2023). “Diffusion models with time-dependent parameters: An analysis of computational effort and accuracy of different numerical methods.” Journal of Mathematical Psychology, 114, 102756. doi:10.1016/j.jmp.2023.102756.

See Also

add_residual

Calculate Quantiles

Description

Function that calls the underlying quantile calculation functions calc_quantiles_obs() and calc_quantiles_pred(). Does input checks and the data wrangling

Usage

calc_quantiles(
  pdf_u = NULL,
  pdf_l = NULL,
  t_vec = NULL,
  dt = NULL,
  rts_u = NULL,
  rts_l = NULL,
  one_cond,
  probs = NULL,
  b_coding,
  skip_if_contr_low = NULL
)

Arguments

Details

if pdf_u and pdf_l are not NULL, returns quantiles for the densities

if rts_u and rts_l are not NULL, returns quantiles for the response times

if all are not NULL, returns both.

Value

a data.frame with "Source", "Cond", "Prob"s, "Quant_<u_label>", "Quant_<l_label>" of type c("quantiles", "sum_dist", "stats_dm", "data.frame")

See Also

new_stats_dm()

Calculate Quantiles

Description

Backend functions to calculate quantiles for RT vectors or pdfs

Usage

calc_quantiles_obs(rts_u, rts_l, one_cond, probs)

calc_quantiles_pred(
  pdf_u,
  pdf_l,
  t_vec,
  one_cond,
  probs,
  dt,
  skip_if_contr_low = NULL
)

Arguments

Details

for RTs: straightforward via stats::quantile.

for Densities: Calculate CDF (for each pdf separately here), and then map the desired probability level via the CDF (y-axis) to the time space (x-axis)

Value

a data.frame with the "Cond" label, the "Prob"s and "Quant_U" and "Quant_L" for the quantiles

RMSE Calculation and Helpers

Description

Internal helpers for computing the root mean squared error (RMSE) between predicted and observed quantiles and conditional accuracy functions.

Usage

calc_rmse_eval(pdfs, t_vec, dt, stats_agg, stats_agg_info, weight_err = 1.5)

calc_rmse(quants_pred, cafs_pred, quants_obs, cafs_obs, weight_err = 1.5)

Arguments

Details

• calc_rmse_eval() prepares observed and predicted quantiles/CAFs from PDFs and aggregated info, then calls calc_rmse().

• calc_rmse() computes the weighted RMSE given predicted and observed quantiles/CAFs.

Value

A single numeric RMSE value, or NULL if no observed stats were provided, or Inf if predictions failed (contain NA).

See Also

stats_from_pdfs_agg_info()

Calculate Statistics

Description

calc_stats provides an interface for calculating statistics/metrics on model predictions and/or observed data. Supported statistics include basic statistics on mean and standard deviation, Conditional Accuracy Functions (CAFs), Quantiles, Delta Functions, and fit statistics. Results can be aggregated across individuals.

Usage

calc_stats(object, type, ...)

## S3 method for class 'data.frame'
calc_stats(
  object,
  type,
  ...,
  conds = NULL,
  resample = FALSE,
  progress = 1,
  level = "individual",
  b_coding = NULL
)

## S3 method for class 'drift_dm'
calc_stats(object, type, ..., conds = NULL, resample = FALSE)

## S3 method for class 'fits_ids_dm'
calc_stats(
  object,
  type,
  ...,
  conds = NULL,
  resample = FALSE,
  progress = 1,
  level = "individual"
)

## S3 method for class 'fits_agg_dm'
calc_stats(
  object,
  type,
  ...,
  conds = NULL,
  resample = FALSE,
  progress = 1,
  level = "group",
  messaging = TRUE
)

## S3 method for class 'stats_dm'
print(
  x,
  ...,
  round_digits = NULL,
  print_rows = NULL,
  some = NULL,
  show_header = NULL,
  show_note = NULL
)

## S3 method for class 'stats_dm_list'
print(x, ...)

Arguments

Details

calc_stats is a generic function to handle the calculation of different statistics/metrics for the supported object types. Per default, it returns the requested statistics/metrics.

List of Supported Statistics

Basic Statistics

With "basic statistics", we refer to a summary of the mean and standard deviation of response times, including a proportion of response choices.

Conditional Accuracy Function (CAFs)

CAFs are a way to quantify response accuracy against speed. To calculate CAFs, RTs (whether correct or incorrect) are first binned and then the percent correct responses per bin is calculated.

When calculating model-based CAFs, a joint CDF combining both the pdf of correct and incorrect responses is calculated. Afterwards, this CDF is separated into even-spaced segments and the contribution of the pdf associated with a correct response relative to the joint CDF is calculated.

The number of bins can be controlled by passing the argument n_bins. The default is 5.

Quantiles

For observed response times, the function stats::quantile is used with default settings.

Which quantiles are calcuated can be controlled by providing the probabilites, probs, with values in [0, 1]. Default is seq(0.1, 0.9, 0.1).

Delta Functions

Delta functions calculate the difference between quantiles of two conditions against their mean:

• Delta_i = Q_{i,j} - Q_{i,k}

• Avg_i = 0.5 \cdot Q_{i,j} + 0.5 \cdot Q_{i,k}

With i indicating a quantile, and j and k two conditions.

To calculate delta functions, users have to specify:

• minuends: character vector, specifying condition(s) j. Must be in conds(drift_dm_obj).

• subtrahends: character vector, specifying condition(s) k. Must be in conds(drift_dm_obj)

• dvs: character, indicating which quantile columns to use. Default is "Quant_<u_label>". If multiple dvs are provided, then minuends and subtrahends must have the same length, and matching occurs pairwise. In this case, if only one minuend/subtrahend is specified, minuend and subtrahend are recycled to the necessary length.

• specifying probs is possible (see Quantiles)

Densities

With "densities", we refer to a summary of the distribution of observed or predicted data. For observed data, histogram values and kernel density estimates are provided. For predicted data, the model's predicted PDFs are provided.

Optional arguments are:

• discr: numeric, the band-width when calculating the histogram or the kernel density estimates. Defaults to 0.015 seconds

• t_max: numeric, the maximum time window when calculating the distribution summaries of observe data. Defaults to the longest RT (for observed data) or the maximum of the time domain of a model (which is the preferred choice, if possible). If necessary, t_max is slightly adjusted to match with discr.

• scale_mass: logical, only relevant if observed data is available. If TRUE, density masses are scaled proportional to the number of trials per condition.

Fit Statistics

Calculates the Log-Likelihood, Akaike and Bayesian Information Criteria, and root-mean squared-error statistic.

Optional arguments are:

• k: numeric, for penalizing the AIC statistic (see also stats::AIC and AIC.fits_ids_dm).

• n_bins, probs: numeric vectors, see the section on CAFs and Quantiles above

• weight_err: numeric scalar, determines how CAFs and quantiles are weighted. Default is 1.5.

Resampling

When resampling = TRUE, an uncertainty interval is provided via simulation. The default number of iterations is R = 100, which can be changed by passing the optional argument R.

If resampling is requested, the returned stats_dm object contains the column "Estimate", coding the interval. The interval width is controlled via the optional argument interval_level, a single numeric value between 0 and 1 (default: 0.95). The interpretation of this interval depends on the specific situation (see below).

Resampling at the Individual Level

If object is a drift_dm object (i.e., a single model created by drift_dm()), synthetic data are simulated under the model, and for each synthetic data set the requested statistic is calculated. The interval then reflects the range of these simulated statistics. To determine the number of trials for each synthetic data set, dRiftDM either uses the observed data attached to the model (if available) or the optional argument n_sim (passed to simulate_data()). Note that n_sim must be provided if no observed data are available, and that n_sim always has priority.

If object is a drift_dm object with attached observed data, resampling is also performed for the observed data. In this case, trials are bootstrapped, and for each bootstrap sample the requested statistic is calculated.

If object is a data.frame, fits_agg_dm, or fits_ids_dm object, resampling is performed for each individual if level = "individual". For both models and observed data, synthetic or bootstrapped data sets are generated as described above.

Resampling at the Group Level

Group-level resampling is possible only if object is a data.frame (with an "ID" column), fits_agg_dm, or fits_ids_dm object. To request this, set level = "group". Participants are then bootstrapped, and for each bootstrapped sample the aggregated statistic is calculated.

Interpretation of Intervals

For level = "group", intervals represent bootstrapped confidence intervals For level = "individual", intervals represent the variability in the statistic when data for a single participant are resampled or simulated under the model.

Note

For objects of type fits_agg_dm, which contain a mixture of group- and individual-level information, the level argument only affects resampling for the observed data. For the model itself, resampling is always performed under the fitted model, in the same way as for a drift_dm object.

Value

If type is a single character string, then a subclass of data.frame is returned, containing the respective statistic. Objects of type sum_dist will have an additional attribute storing the boundary encoding (see also b_coding). The reason for returning subclasses of data.frame is to provide custom plot() methods (e.g., plot.cafs). To get rid of the subclass label and additional attributes (i.e., to get just the plain underlying data.frame, users can use unpack_obj()).

If type contains multiple character strings (i.e., is a character vector) a subclass of list with the calculated statistics is returned. The list will be of type stats_dm_list (to easily create multiple panels using the respective plot.stats_dm_list() method).

The print methods print.stats_dm() and print.stats_dm_list() each invisibly return the supplied object x.

Note

When a model's predicted density function integrates to a value of less than drift_dm_skip_if_contr_low(), means and quantiles return the values NA. Users can alter this by explicitly passing the argument skip_if_contr_low when calling calc_stats() (e.g., calc_stats(..., skip_if_contr_low = -Inf))

Examples

# Example 1: Calculate CAFs and Quantiles from a model ---------------------
# get a model for demonstration purpose
a_model <- ssp_dm()
# and then calculate cafs and quantiles
some_stats <- calc_stats(a_model, type = c("cafs", "quantiles"))
print(some_stats)

# Example 2: Calculate a Delta Function from a data.frame ------------------
# get a data set for demonstration purpose
some_data <- ulrich_simon_data
conds(some_data) # relevant for minuends and subtrahends
some_stats <- calc_stats(
  a_model,
  type = "delta_funs",
  minuends = "incomp",
  subtrahends = "comp"
)
print(some_stats, print_rows = 5)


# Example 3: Calculate Quantiles from a fits_ids_dm object -----------------
# get an auxiliary fits_ids_dm object
all_fits <- get_example_fits("fits_ids_dm")
some_stats <- calc_stats(all_fits, type = "quantiles")
print(some_stats, print_rows = 5) # note the ID column

# one can also request that the statistics are averaged across individuals
print(
  calc_stats(all_fits, type = "quantiles", average = TRUE)
)

Calculate Statistics for Model Prediction and/or Observed Data

Description

This function derives statistics that can be calculated for model predictions and/or observed data. However, it does not calculate it, but rather calls the respective backend functions. Supported statistics currently include:

• Basic Summary Statistics (i.e., means and response percentages calc_basic_stats())

• Conditional Accuracy Functions (CAFs; calc_cafs())

• Quantiles (calc_quantiles())

• Delta Functions (calc_delta_funs()).

• Density Estimates (calc_dens()).

Usage

calc_stats_pred_obs(type, b_coding, conds, ..., scale_mass = FALSE)

Arguments

Details

When calling this function the arguments all_rts_u/all_rts_l and/or all_pdfs must always be specified (see re_evaluate_model, obs_data). Otherwise, the backend functions won't work properly. Further arguments are:

• for CAFS: n_bins controls the number of bins, with a default of 5.

• for Quantiles and Delta Functions: probs controls the quantiles to calculate. Default is seq(0.1, 0.9, 0.1) (see drift_dm_default_probs()).

• for basic summary satistics, Quantiles, and Delta Function: skip_if_contr_low controls if quantiles and means are calculated for PDFs with very small contribution (see also drift_dm_skip_if_contr_low()).

• for densities: discr controls the bin width for the observed data. Default is 0.015 seconds

This function gets called by calc_stats()

Value

A data frame with the calculated statistic across conds (ordered according to Source).

Metropolis Acceptance Step for Vectorized MCMC Sampling

Description

This internal function computes the Metropolis acceptance decision for a set of MCMC proposals, using either hierarchical or non-hierarchical posterior evaluation. It returns updated log-posterior and log-likelihood values, as well as a logical vector indicating accepted proposals.

Usage

call_log_posterior_m(
  proposal_mat,
  prev_prms_mat,
  prev_pis,
  prev_lls,
  level,
  re_eval,
  ...
)

Arguments

Details

This function implements a vectorized Metropolis acceptance step for multiple MCMC chains simultaneously. The posterior is calculated using either log_posterior_lower() for subject-level parameters or log_posterior_hyper() for group-level parameters.

Log-posterior and log-likelihood values are only updated where proposals were accepted. In cases where proposals yield invalid posteriors (i.e., NA), they are automatically rejected.

Value

A list with three elements:

pis

A numeric vector of updated log-posterior values.

lls

A numeric vector of updated log-likelihood values.

accept

A logical vector of length equal to the number of chains, indicating which proposals were accepted.

Check a B_Coding

Description

Checks if a list satisfies the requirements to be considered a valid b_coding

Usage

check_b_coding(b_coding)

Arguments

Details

Checks for...

• input being a list

• list is of length three and provides the names "column", "u_name_value", "l_name_value",

• if b_coding$column provides a single string

• if b_coding$u_name_value and b_coding$l_name_value provides a single named value of type character or numeric, and that both are of the same type

Value

the unmodified list for convenience

Checks if a variable/vector of length 1 is a number > 0 or 0/expression

Description

assumed that x is either a digit or an expression

Usage

check_digit_larger_0(x)

Arguments

Value

a single logical

Check time/space discretization via reference comparison

Description

check_discretization() helps you choose or check time (dt) and space (dx) discretization settings. It computes a high-precision reference solution of the model's PDFs with dt_ref/dx_ref, and then compares the reference PDFs to the discretization settings of the supplied object, using the Hellinger distance per condition. Smaller distances indicate closer agreement with the reference — i.e., a sufficiently fine grid.

There are not yet overall and officially published recommendations on how large the Hellinger distance can be without affecting model precision, and this very likely will depend on the model itself. Based on some preliminary simulations using dmc_dm(), we would recommend trying to keep the Hellinger Distance at best below 10 percent. However, we also observed for extreme parameter values that the Hellinger distance can be even larger without sacrificing the qualitative model behavior, and vice versa! It is thus best to iterate between plotting model predictions and calculating the Hellinger Distance, to ensure that you can best interpret this quantity for your model at hand. Furthermore, we recommend to run parameter recoveries using simulate_data() and estimate_dm(), to check if you can recover data generated under your model with fine discretization using that same model with coarse discretization.

Usage

check_discretization(object, ...)

## S3 method for class 'drift_dm'
check_discretization(
  object,
  ...,
  dt_ref = 0.001,
  dx_ref = 0.001,
  round_digits = 5
)

## S3 method for class 'fits_ids_dm'
check_discretization(object, ...)

## S3 method for class 'fits_agg_dm'
check_discretization(object, ...)

Arguments

Details

Under the hood, for each condition, we concatenate the lower- and upper- boundary PDFs (pdf_l, pdf_u), interpolate the model PDFs to a time space matching with the reference PDFs, and then compute the Hellinger distance: H(p,q) = \sqrt{1 - \int \sqrt{p(t)\,q(t)}\,dt}

There are not yet overall, officially published recommendations on how large the Hellinger distance can be without affecting model precision, and this may even depend on the specific model. Based on preliminary simulations, we recommend trying to keep the average Hellinger distance below 5\

The reference discretizations (dt_ref/dx_ref) must be at least as fine as the object's current discretization settings (dt_model/dx_model). If dt_model < dt_ref or dx_model < dx_ref, an error is raised because the “reference” would not be the finest solution.

Value

a named numeric vector of Hellinger distances (one per condition) if object is of type drift_dm or fits_agg_dm. A data.frame of Hellinger distances across IDs and conditions if object is of type fits_ids_dm. Hellinger distances are in ⁠[0, 1]⁠, where 0 means identical to the reference.

See Also

estimate_dm(), trapz()

Examples

# Example:
my_model <- ratcliff_dm()

# Assess current (dt=0.0075, dx=0.02) against a fine reference:
check_discretization(my_model)

# If distances are near zero across conditions, the current grid is adequate.

Check if Object is a Named Numeric Vector

Description

Validates that an object is a named numeric vector with specified attributes. Optionally checks specific names, length, and restrictions on label characters.

Usage

check_if_named_numeric_vector(
  x,
  var_name,
  labels = NULL,
  length = NULL,
  allow_non_word_chars = FALSE
)

Arguments

Details

Checks for:

• Numeric type of x with non-zero length

• Required length, if specified

• Unique, non-empty names for each entry in x

• Match of all names in x to labels, if labels is specified

• Absence of NA of Inf values in x

• Optional absence of non-word names if allow_non_word_chars is FALSE

Value

Throws an error if the conditions are not met. If all checks pass, no output is returned.

checks if all entries of internal_list are an expression or integer. Throws an error if not (internal docu)

Description

Checks also if there are entries for each parameter and condition

Usage

check_internal_list(internal_list, prm_names, cond_names)

Arguments

Value

the internal list for convenience

Check and Reduce the Observed Data

Description

Checks a data set that is considered an "observed data set". Used in the internals of dRiftDM. When calling this function, unncessary column names are stripped away.

Usage

check_reduce_raw_data(obs_data, b_coding_column, u_value, l_value)

Arguments

Details

Checks:

• if obs_data is a data.frame

• if "RT", b_coding_column, and "Cond" column are present

• when IDs are present, if each ID has values on each condition. At the same time unused factor levels are dropped drop_levels_ID_column

• If all columns are there, the data set is reduced to the relevant ones

• for missing Values, and drops rows with missing values

• if "Cond" is of type character, and if not casts it to character

• if RT is of type numeric, and of not casts it to numeric

• RTs are >= 0

• that the values in b_coding_column match with u_value and l_value (casts the column if necessary)

• if b_coding_column has only 1 or 2 unique values

Value

the obs_data for convenience (with edits as listed under Details).

Access Coefficients of a Model

Description

Extract or set the coefficients/parameters objects supported by dRiftDM.

Usage

coef(object, ...) <- value

## S3 replacement method for class 'drift_dm'
coef(object, ..., eval_model = FALSE) <- value

## S3 method for class 'drift_dm'
coef(object, ..., select_unique = TRUE, select_custom_prms = TRUE)

## S3 method for class 'fits_agg_dm'
coef(object, ...)

## S3 method for class 'fits_ids_dm'
coef(object, ...)

## S3 method for class 'mcmc_dm'
coef(object, ..., .f = mean, id = NULL)

## S3 method for class 'coefs_dm'
print(
  x,
  ...,
  round_digits = drift_dm_default_rounding(),
  print_rows = 10,
  some = FALSE,
  show_header = TRUE,
  show_note = TRUE
)

Arguments

Details

⁠coef.*()⁠ are methods for the generic stats::coef() function; ⁠coefs<-()⁠ is a generic replacement function, currently supporting objects of type drift_dm.

The argument value supplied to the ⁠coefs<-()⁠ function must match with the vector returned from ⁠coef(<object>)⁠. It is possible to update just part of the (unique) parameters.

Whenever the argument select_unique is TRUE, dRiftDM tries to provide unique parameter labels.

Value

For objects of type drift_dm, coefs() returns either a named numeric vector if select_unique = TRUE, or a matrix if select_unique = FALSE. If select_custom_prms = TRUE, custom parameters are added to the matrix.

For objects of type fits_ids_dm, coefs() returns a data.frame. If select_unique = TRUE, the columns will be the (unique, free) parameters, together with a column coding IDs. If select_unique = FALSE, the columns will be the parameters as listed in the columns of prms_matrix (see drift_dm), together with columns coding the conditions and IDs. If select_custom_prms = TRUE, the data.frame will also contain columns for the custom parameters. The returned data.frame has the class label coefs_dm to easily plot histograms for each parameter (see hist.coefs_dm).

For objects of type fits_agg_dm, returns the same as coef.drift_dm() (i.e., as if calling coef() with an object of type drift_dm)

For objects of type mcmc_dm, the return type depends on the model structure and the .f output:

• If the model is non-hierarchical or id is a single value (not NA), the function returns either a vector or a matrix, depending on whether .f returns a single value or a vector.

• In the hierarchical case, when id is a vector or NA, the function returns a data.frame. If .f returns a single value, the data.frame will contain one row per participant (with an ID column and one column per parameter). If .f returns a vector, the data.frame will include an additional column .f_out, coding the output of .f in long format.

See Also

drift_dm()

Examples

# get a pre-built model and a data set for demonstration purpose
# (when creating the model, set the discretization to reasonable values)
a_model <- dmc_dm()
coef(a_model) # gives the free and unique parameters
coef(a_model, select_unique = FALSE) # gives the entire parameter matrix

The Component Functions of A Model

Description

Functions to get or set the "component functions" of an object. The component functions are a list of functions providing the drift rate, boundary, starting point distribution, and non-decision time distribution They are at the heart of the package and shape the model's behavior.

Usage

comp_funs(object, ...) <- value

## S3 replacement method for class 'drift_dm'
comp_funs(object, ..., eval_model = FALSE) <- value

comp_funs(object, ...)

## S3 method for class 'drift_dm'
comp_funs(object, ...)

## S3 method for class 'fits_ids_dm'
comp_funs(object, ...)

## S3 method for class 'fits_agg_dm'
comp_funs(object, ...)

Arguments

Details

comp_funs() is a generic accessor function, and ⁠comp_funs<-()⁠ is a generic replacement function. The default methods get and set the "component functions". The component functions are a list of functions, with the following names (see also vignette("customize_ddms", "dRiftDM") for examples):

• mu_fun and mu_int_fun, provide the drift rate and its integral, respectively, across the time space.

• x_fun provides a distribution of the starting point across the evidence space.

• b_fun and dt_b_fun provide the values of the upper decision boundary and its derivative, respectively, across the time space. It is assumed that boundaries are symmetric.

• nt_fun provides a distribution of the non-decision component across the time space.

All of the listed functions are stored in the list comp_funs of the respective model (see also drift_dm()).

Each component function must take the model's parameters (i.e., one row of prms_matrix), the parameters for deriving the PDFs, the time or evidence space, a condition, and a list of optional values as arguments. These arguments are provided with values when dRiftDM internally calls them.

In order to work with dRiftDM, mu_fun, mu_int_fun, b_fun, dt_b_fun, and nt_fun must have the following declaration: ⁠my_fun = function(prms_model, prms_solve, t_vec, one_cond, ddm_opts⁠). Here, prms_model is one row of prms_matrix, prms_solve the parameters relevant for deriving the PDFs, t_vec the time space, going from 0 to t_max with length nt + 1 (see drift_dm), and one_cond a single character string, indicating the current condition. Finally dmm_opts may contain additional values. Each function must return a numeric vector of the same length as t_vec. For mu_fun, mu_int_fun, b_fun, dt_b_fun the returned values provide the respective boundary/drift rate (and their derivative/integral) at every time step t. For nt_fun the returned values provide the density of the non-decision time across the time space (which get convoluted with the pdfs when solving the model)

In order to work with dRiftDM, x_fun must have the following declaration: ⁠my_fun = function(prms_model, prms_solve, x_vec, one_cond, ddm_opts⁠). Here, x_vec is the evidence space, going from -1 to 1 with length nx + 1 (see drift_dm). Each function must return a numeric vector of the same length as x_vec, providing the density values of the starting points across the evidence space.

Drift rate and its integral:

The drift rate is the first derivative of the expected time-course of the diffusion process. For instance, if we assume that the diffusion process X is linear with a slope of v...

E(X) = v \cdot t

...then the drift rate at every time step t is the constant v, obtained by taking the derivative of the expected time-course with respect to t:

\mu(t) = v

Conversely, the integral of the drift rate is identical to the expected time-course:

\mu_{int}(t) = v \cdot t

For the drift rate mu_fun, the default function when calling drift_dm() is a numeric vector containing the number 3. Its integral counterpart mu_int_fun will return a numeric vector containing the values t_vec*3.

Starting Point Distribution:

The starting point of a diffusion model refers to the initial value taken by the evidence accumulation process at time t=0. This is a PDF over the evidence space.

The default function when calling drift_dm() will be a function returning a dirac delta on zero, meaning that every potential diffusion process starts at 0.

Boundary:

The Boundary refers to the values of the absorbing boundaries at every time step t in a diffusion model. In most cases, this will be a constant. For instance:

b(t) = b

In this case, its derivative with respect to t is 0.

The default function when calling drift_dm() will be function for b_fun returning a numeric vector of length length(t_vec) containing the number 0.5. Its counterpart dt_b will return a numeric vector of the same length containing its derivative, namely, 0.

Non-Decision Time:

The non-decision time refers to an additional time-requirement. Its distribution across the time space will be convoluted with the PDFs derived from the diffusion process.

In psychology, the non-decision time captures time-requirements outside the central decision process, such as stimulus perception and motor execution.

The default function when calling drift_dm() returns a dirac delta on t = 0.3.

Value

For comp_funs() the list of component functions.

For ⁠comp_funs<-()⁠ the updated drift_dm object.

Note

There is only a replacement function for drift_dm objects. This is because replacing the component functions after the model has been fitted (i.e., for a fits_ids_dm object) doesn't make sense.

See Also

drift_dm()

Examples

# get a pre-built model for demonstration
my_model <- ratcliff_dm()
names(comp_funs(my_model))

# direct replacement (see customize_ddms for a more information on
# how to write custom component functions)
# 1. Choose a uniform non-decision time from the pre-built component_shelf()
nt_uniform <- component_shelf()$nt_uniform
# swap it in
comp_funs(my_model)[["nt_fun"]] <- nt_uniform

# now update the flex_prms object to ensure that this model has the required
# parameters
prms <- c(muc = 3, b = 0.6, non_dec = 0.3, range_non_dec = 0.05)
conds <- "null"
new_flex_prms <- flex_prms(prms, conds = conds)
flex_prms(my_model) <- new_flex_prms

# accessor method also available for fits_ids_dm objects
# (see estimate_model_ids)
# get an exemplary fits_ids_dm object
fits <- get_example_fits("fits_ids_dm")
names(comp_funs(fits))

Evaluate all Component Functions

Description

Gets/calculates all values provided by the component functions of a drift_dm object

Usage

comp_vals(
  drift_dm_obj,
  x_vec = NULL,
  t_vec = NULL,
  nt = NULL,
  dt = NULL,
  nx = NULL,
  dx = NULL,
  prms_solve = NULL,
  solver = NULL,
  prms_matrix = NULL
)

Arguments

Details

arguments are optional, because they can be extracted from the model. However, supplying these are faster than creating them.

Value

If solver "kfe", a named list with entries "mu_vals", "x_vals", "b_vals", "dt_b_vals", "nt_vals".

If solver "im_zero", the returned list will also contain "mu_int_vals".

Diffusion Model Components

Description

This function is meant as a convenient way to access pre-built model component functions.

Usage

component_shelf()

Details

The function provides the following functions:

• mu_constant, provides the component function for a constant drift rate with parameter muc.

• mu_dmc, provides the drift rate of the superimposed diffusion process of DMC. Necessary parameters are muc (drift rate of the controlled process), a (shape..), A (amplitude...), tau (scale of the automatic process).

• mu_ssp, provides the drift rate for SSP. Necessary parameters are p (perceptual input of flankers and target), sd_0 (initial spotlight width), r (shrinking rate of the spotlight) and 'sign' (an auxiliary parameter for controlling the contribution of the flanker stimuli). Note that no mu_int_ssp exists.

• mu_int_constant, provides the complementary integral to mu_constant.

• mu_int_dmc, provides the complementary integral to mu_dmc.

• x_dirac_0, provides a dirac delta for a starting point centered between the boundaries (no parameter required).

• x_uniform, provides a uniform distribution for a start point centered between the boundaries. Requires a parameter range_start (between 0 and 2).

• x_beta, provides the function component for a symmetric beta-shaped starting point distribution with parameter alpha.

• b_constant, provides a constant boundary with parameter b.

• b_hyperbol, provides a collapsing boundary in terms of a hyperbolic ratio function with parameters b0 as the initial value of the (upper) boundary, kappa the size of the collapse, and t05 the point in time where the boundary has collapsed by half.

• b_weibull, provides a collapsing boundary in terms of a Weibull distribution with parameters b0 as the initial value of the (upper) boundary, lambda controlling the time of the collapse, k the shape of the collapse, and kappa the size of the collapse.

• dt_b_constant, the first derivative of b_constant.

• dt_b_hyperbol, the first derivative of b_hyperbol.

• nt_constant, provides a constant non-decision time with parameter non_dec.

• nt_uniform, provides a uniform distribution for the non-decision time. Requires the parameters non_dec and range_non_dec.

• nt_truncated_normal, provides the component function for a normally distributed non-decision time with parameters non_dec, sd_non_dec. The Distribution is truncated to [0, t_{max}].

• dummy_t a function that accepts all required arguments for mu_fun or mu_int_fun but which throws an error. Might come in handy when a user doesn't require the integral of the drift rate.

See vignette("customize_ddms", "dRiftDM") for more information on how to set/modify/customize the components of a diffusion model.

Value

A list of the respective functions; each entry/function can be accessed by "name" (see the Example and Details).

Examples

pre_built_functions <- component_shelf()
names(pre_built_functions)

The Conditions of an Object

Description

Extract the conditions from a (supported) object.

Usage

conds(object, ...) <- value

## S3 replacement method for class 'drift_dm'
conds(object, ..., eval_model = FALSE, messaging = TRUE) <- value

conds(object, ...)

## S3 method for class 'drift_dm'
conds(object, ...)

## S3 method for class 'fits_ids_dm'
conds(object, ...)

## S3 method for class 'fits_agg_dm'
conds(object, ...)

## S3 method for class 'data.frame'
conds(object, ...)

## S3 method for class 'traces_dm_list'
conds(object, ...)

Arguments

Details

conds() is a generic accessor function and ⁠conds<-()⁠ is a generic replacement function. The replacement method currently only supports drift_dm objects. The default methods get and set the conditions of an object.

When replacing the conditions of a drift_dm object, a new flex_prms object is created and then set to the model, resetting all parameter specifications and setting all parameter values to those of the previously first condition. In addition, if data was attached to the model, the data is removed. This is because there is no meaningful way for dRiftDM to know how the model should behave for the newly introduced condition(s), and how these new conditions relate to the old ones. Messages reminding the user of this behavior are displayed per default.

Value

For conds() NULL or a character vector with the conditions. NULL is given if the object has no conditions (e.g., when a data.frame has no Cond column).

For ⁠conds<-()⁠ the updated drift_dm object.

See Also

drift_dm()

Examples

# get a pre-built model to demonstrate the conds() function
my_model <- dmc_dm()
conds(my_model)

# accessor functions also work with other object types provided by dRiftDM
# (simulated traces; see the documentation of the respective function)
some_traces <- simulate_traces(my_model, k = 1)
conds(some_traces)

# get an exemplary fits_ids_dm object (see estimate_model_ids)
fits <- get_example_fits("fits_ids_dm")
conds(fits)

# also works with data.frames that have a "Cond" column
conds(dmc_synth_data)

Copy Class Attributes from One Object to Another

Description

This function transfers class attributes from an old object to a new object, ensuring that new inherits the class structure and missing attributes of old. The primary purpose is to enforce class consistency and restore any lost attributes when modifying or combining objects. It is used in the internals of the package and it is not exported.

Usage

copy_class_attributes(old, new)

## S3 method for class 'stats_dm'
copy_class_attributes(old, new)

Arguments

Details

The function assumes that all class attributes of new can be found in old. Note also, that the order of attributes is not ensured.

Value

The modified new object with attributes and class from old.

Access/Replace the Cost Function Label and Access the Cost Function Value

Description

Functions to access/replace the cost function label of a ⁠dRiftDM object⁠ and to access the current cost function value. The cost function label codes which cost function is used during estimation (e.g., the negative log-likelihood). The cost function value indicates the current value of the cost function given the current set of parameters and the data.

Usage

cost_function(object, ...) <- value

## S3 replacement method for class 'drift_dm'
cost_function(object, ..., eval_model = FALSE) <- value

cost_function(object, ...)

## S3 method for class 'drift_dm'
cost_function(object, ...)

## S3 method for class 'fits_ids_dm'
cost_function(object, ...)

## S3 method for class 'fits_agg_dm'
cost_function(object, ...)

cost_value(object, ...)

## S3 method for class 'drift_dm'
cost_value(object, ...)

## S3 method for class 'fits_ids_dm'
cost_value(object, ...)

## S3 method for class 'fits_agg_dm'
cost_value(object, ...)

Arguments

Value

• cost_function() returns a single character string, specifying the used cost function

• ⁠cost_function<-()⁠ returns the model object with the updated cost function.

• cost_value() returns a single numeric if object is of type drift_dm or fits_agg_dm. If there is no data attached to an object of type drift_dm, the function returns NULL. If object is of type fits_ids_dm, the function returns a data.frame with all cost values across participants.

See Also

drift_dm(), re_evaluate_model()

Examples

# get a pre-built model for demonstration purpose
a_model <- ratcliff_dm(obs_data = ratcliff_synth_data)
cost_function(a_model)
cost_value(a_model)

# switch the default cost function to rmse
cost_function(a_model) <- "rmse"
out <- estimate_dm(a_model, verbose = 0, messaging = FALSE)
# -> the model was estimated using the RMSE statistic

Count the number of digits > 0

Description

This function takes one entry of the internal_list (i.e., all conditions for one parameter) and counts how often there are unique digits > 0

Usage

count_unique_prms_one_internal_entry(one_internal_entry)

Arguments

Value

an integer number

Create a matrix for lower and upper

Description

Outsourced, deep inside the package function to avoid large nesting

Usage

create_matrix_smart(input, conds, prm_labels = NULL)

Arguments

Details

The goal of this function is to build up a matrix, serving as the upper or lower end of a parameter space (relevant when simulating data). The function gets called by get_parameters_smart().

It assumes the following: input is either a list or a numeric vector.

• The easiest case is when it is a numeric vector. In this case, the function builds a matrix with as many rows as entries in conds. The rows will also be labeled according to conds. The column names are either the names specified with the numeric vector, or the labels specified in prm_labels

• The less intuitive case is when input is a list. In this case, the list requires an entry called "default_values" which specifies the named or plain numeric vector as above. If the list only contains this entry, then the behavior is as if input was already a numeric vector. However, the input list can also have entries labeled as specific conditions, which contain named (!) numeric vectors with parameter labels. This will modify the value for the upper/lower parameter space with respect to the specified parameters in the respective condition.#'

Value

a matrix indicating either the upper or lower end of a parameter space. There will be as many rows as conds implies. The number of columns depend on input (matching its length if it is a vector, or matching the length of the entry "default_values" if it is a list). If input is NULL, then NULL is returned.

See Also

simulate_data(), simulate_values()

Create "Temperatures" for TIDE

Description

Create "Temperatures" for TIDE

Usage

create_temperatures(n_chains, sampler)

Arguments

Value

a numeric vector of length equal to n_chains. The returned values correspond to quantiles of a Beta(0.3, 1) distribution for sampler == "TIDE". Otherwise a numeric vector of 1s is returned.

Perform Crossover Between Chains

Description

This function dispatches to either full_crossover() or migration_crossover() depending on the which argument.

Usage

crossover(which, ...)

Arguments

Value

A list with the following components:

• new_prms_across_chains: The updated parameter matrix of shape ⁠p × n⁠.

• new_pis_across_chains: The updated vector of log-posterior values.

• new_log_likes_across_chains: The updated vector of log-likelihood values.

Default Prior for Group-Level (Hyper) Parameters

Description

These functions define and evaluate a default prior distribution for hyperparameters at the group level.

Usage

d_default_prior_hyper(x, mean, sd, lower, upper, shape, rate, log)

r_default_prior_hyper(n, mean, sd, lower, upper, shape, rate)

Arguments

Details

d_default_prior_hyper computes the (log) density of a prior for a two-element vector or a 2xN matrix, containing the mean and standard deviation (i.e., phi_j). The mean is modeled with a truncated normal distribution, and the standard deviation with a gamma distribution.

r_default_prior_hyper samples hyperparameter values from this prior.

the arguments mean, sd, lower, upper, shape, and rate are recycled if necessary with respect to the columns of x. For example, if x has two columns, then mean might provide two values.

Value

For d_default_prior_hyper, a numeric vector representing the (log) prior density value(s), with the simplifying assumption of independence of the mean and standard deviation.

For r_default_prior_hyper, a 2-row matrix with n columns. The first row contains sampled group means; the second row contains sampled standard deviations. Samples are drawn independently. If n is 1, then a named numeric vector is returned.

Optional Arguments for the Component Functions

Description

Functions to get or set the optional, user-defined R objects attached to a model object.

Usage

ddm_opts(object, ...) <- value

## S3 replacement method for class 'drift_dm'
ddm_opts(object, ..., eval_model = FALSE) <- value

ddm_opts(object, ...)

## S3 method for class 'drift_dm'
ddm_opts(object, ...)

## S3 method for class 'fits_agg_dm'
ddm_opts(object, ...)

Arguments

Details

When deriving model predictions, the model's component functions (see comp_funs()) are evaluated and the returned values are passed forward to dedicated numerical methods implemented in dRiftDM. To allow users to access arbitrary R objects within their custom component functions, models may contain a ddm_opts entry (see also drift_dm() and the end of vignette("customize_ddms", "dRiftDM") for an example).

ddm_opts() is a generic accessor function, and ⁠ddm_opts<-()⁠ is a generic replacement function. The default methods get and set the optional R object.

Value

For ddm_opts() the optional R object that was once supplied by the user, or NULL.

For ⁠ddm_opts<-()⁠ the updated drift_dm object.

See Also

drift_dm(), comp_funs()

Examples

# get a pre-built model for demonstration
a_model <- ratcliff_dm()
ddm_opts(a_model) <- "Hello World"
ddm_opts(a_model)

Default Values for the dRiftDM Package

Description

These functions provide default values for various settings in the dRiftDM package.

Usage

drift_dm_approx_error()

drift_dm_medium_approx_error()

drift_dm_small_approx_error()

drift_dm_rough_approx_error()

drift_dm_robust_prm()

drift_dm_default_rounding()

drift_dm_default_probs()

drift_dm_default_n_bins()

drift_dm_default_b_coding()

drift_dm_skip_if_contr_low()

drift_dm_n_id_trunc_warn()

Details

• drift_dm_approx_error(): Returns the default approximation error for precise calculations (1e-20).

• drift_dm_medium_approx_error(): Returns a 'medium' level of approximation error (1e-06).

• drift_dm_small_approx_error(): Returns a 'small' level of approximation error (.01).

• drift_dm_rough_approx_error(): Returns a rough level of approximation error (.1).

• drift_dm_robust_prm(): Returns a value that is added to the PDFs after convolution with the non-decision time to make parameter estimation and the evaluation of the log-likelihood more robust (1e-8).

• drift_dm_default_rounding(): Returns the default rounding precision for numerical outputs (3).

• drift_dm_default_probs(): Returns the default sequence of probabilities for quantiles (0.1, 0.2, ..., 0.9)

• drift_dm_default_n_bins(): Returns the default number of bins for a CAF (5)

• drift_dm_default_b_coding(): Returns the default boundary coding (list(column = "Error", u_name_value = c("corr" = 0), l_name_value = c("err" = 1))

• drift_dm_skip_if_contr_low(): returns the value 0.0001. If a PDF integrates to a value lower than that (i.e., if there is almost no contribution of a PDF; most likely this will be pdf_l), then summary functions returned by calc_stats() might contain the value NA for the respective PDF.

• drift_dm_n_id_trunc_warn(): returns 15. If there are warnings relevant to multiple participants, the printed IDs are truncated at 15.

Value

the respective values/lists as described in the Details section

Create the Diffusion Model for Conflict Tasks

Description

This function creates a drift_dm object that corresponds to the Diffusion Model for Conflict Tasks by Ulrich et al. (2015).

Usage

dmc_dm(
  var_non_dec = TRUE,
  var_start = TRUE,
  instr = NULL,
  obs_data = NULL,
  sigma = 1,
  t_max = 3,
  dt = 0.0075,
  dx = 0.02,
  b_coding = NULL
)

Arguments

Details

The Diffusion Model for Conflict Tasks is a model for describing conflict tasks like the Stroop, Simon, or flanker task.

It has the following properties (see component_shelf):

• a constant boundary (parameter b)

• an evidence accumulation process that results from the sum of two subprocesses:

  • a controlled process with drift rate muc

  • a gamma-shaped process with a scale parameter tau, a shape parameter a, and an amplitude A.

If var_non_dec = TRUE, a (truncated) normally distributed non-decision with mean non_dec and standard deviation sd_non_dec is assumed. If var_start = TRUE, a beta-shaped starting point distribution is assumed with shape and scale parameter alpha.

If var_non_dec = TRUE, a constant non-decision time at non_dec is set. If var_start = FALSE, a starting point centered between the boundaries is assumed (i.e., a dirac delta over 0).

Per default the shape parameter a is set to 2 and not allowed to vary. This is because the derivative of the scaled gamma-distribution function does not exist at t = 0 for a < 2. Currently, we recommend keeping a fixed to 2. If users decide to set a != 2, then a small value of tol = 0.001 (default) is added to the time vector t_vec before calculating the derivative of the scaled gamma-distribution as originally introduced by Ulrich et al. (2015). Users can control this value by passing a value via ddm_opts() (see the example below). Note, however, that varying a can lead to large numerical inaccuracies if a gets smaller.

The model assumes the amplitude A to be negative for incompatible trials. Also, the model contains the custom parameter peak_l, containing the peak latency ((a-2)*tau).

Value

An object of type drift_dm (parent class) and dmc_dm (child class), created by the function drift_dm().

Note

The scaling of the parameters in dRiftDM is different to Ulrich et al. (2015). This is because dRiftDM works in seconds and with a diffusion constant of 1, while the original DMC parameterization is in milliseconds and with a diffusion constant of 4. We describe how to convert the parameters on our website.

References

Ulrich R, Schröter H, Leuthold H, Birngruber T (2015). “Automatic and controlled stimulus processing in conflict tasks: Superimposed diffusion processes and delta functions.” Cognitive Psychology, 78, 148–174. doi:10.1016/j.cogpsych.2015.02.005.

Examples

# the model with default settings
my_model <- dmc_dm()

# the model with no variability in the starting point and a finer
# discretization
my_model <- dmc_dm(var_start = FALSE, dt = .005, dx = .01)

# we don't recommend this, but if you really want a != 2, just do...
# (see the Details for more warnings/information about this)
my_model <- dmc_dm(instr = "a ~!")
coef(my_model)["a"] <- 1.9
# -> if you want to control the small value that is added to t_vec when
# calculating the drift rate for a != 2, just use ...
ddm_opts(my_model) <- 0.0001 # ==> t_vec + 0.0001
ddm_opts(my_model) <- NULL # default ==> t_vec + 0.001

A synthetic data set with two conditions

Description

This dataset was simulated by using the Diffusion Model for Conflict tasks (see dmc_dm()) with parameter settings that are typical for a Simon task.

Usage

dmc_synth_data

Format

A data frame with 600 rows and 3 columns:

RT

Response Times

Error

Error Coding (Error Response = 1; Correct Response = 0)

Cond

Condition ('comp' and 'incomp')

Temporarily suppress new stats generation during resampling

Description

Internal helper that sets the skip_new_stats_dm and skip_validate_stats_dm options to TRUE before evaluating an expression, and resets them to NULL afterward (see also stats.options()). Intended to prevent the (unncessary) creation/checking of stats_dm objects during resampling.

Usage

do_resampling(x)

Arguments

Value

The result of evaluated x.

Draw Samples Using Inverse Transform Sampling

Description

draw_from_pdf generates samples from a given probability density function (PDF) using inverse transform sampling. This function takes in a numeric PDF vector and a corresponding domain vector, then returns a specified number of samples.

Usage

draw_from_pdf(a_pdf, x_def, k, seed = NULL, round_to = NULL, method = "discr")

Arguments

Details

This function implements inverse transform sampling by first constructing a cumulative distribution function (CDF) from the given PDF. Then k values between zero and one are sampled from a uniform distribution, and the corresponding values are mapped to x_def using linear interpolation.

Value

A numeric vector of length k containing the sampled values from the specified PDF. If k is 0, an empty numeric vector is returned.

Create a drift_dm object

Description

This function creates an object of type drift_dm, which serves as the parent class for all further created drift diffusion models (all of which have a child class label, e.g., dmc_dm). The objects created by drift_dm() are the backbone of the dRiftDM package. For a list of all pre-built models, see vignette("dRiftDM", "dRiftDM").

Usage

drift_dm(
  prms_model,
  conds,
  subclass,
  instr = NULL,
  obs_data = NULL,
  sigma = 1,
  t_max = 3,
  dt = 0.001,
  dx = 0.001,
  solver = "kfe",
  cost_function = "neg_log_like",
  mu_fun = NULL,
  mu_int_fun = NULL,
  x_fun = NULL,
  b_fun = NULL,
  dt_b_fun = NULL,
  nt_fun = NULL,
  b_coding = NULL
)

## S3 method for class 'drift_dm'
print(x, ..., round_digits = drift_dm_default_rounding())

Arguments

Details

To modify the entries of a model users can use the replacement methods and the modify_flex_prms() method (see also vignette("dRiftDM", "dRiftDM") and vignette("customize_ddms", "dRiftDM")).

Value

For drift_dm(), a list with the parent class label "drift_dm" and the child class label ⁠<subclass>⁠. The list contains the following entries:

• An instance of the class flex_prms for controlling the model parameters. Provides information about the number of parameters, conditions etc.

• Parameters used for deriving the model predictions, prms_solve, containing the diffusion constant (sigma), the maximum of the time space (t_max), the evidence and space discretization (dt and dx, respectively), and the resulting number of steps for the time and evidence space discretization (nt and nx, respectively).

• A character string solver, indicating the method for deriving the model predictions.

• A character string cost_function, indicating the cost function used for model estimation.

• A list of functions called comp_funs, providing the components of the diffusion model (i.e., mu_fun, mu_int_fun, x_fun, b_fun, dt_b_fun, nt_fun). These functions are called in the depths of the package and will determine the behavior of the model

If (optional) observed data were passed via obs_data(), the list will contain an entry obs_data. This is a (nested) list with stored response times for the upper and lower boundary and with respect to each condition. If the cost function is a summary statistic requiring quantiles, CAFs, etc., the model also contains the entries stats_agg and stats_agg_info. The former is a (nested) list with descriptive statistics. The latter contains information about the descriptive statistics (e.g., the quantile levels).

If the model has been evaluated (see re_evaluate_model()), the list will contain...

• ... the cost value; can be addressed via cost_value().

• ... the PDFs of the first passage time; can be addressed via pdfs().

If the model was estimated (which includes its evaluation), the list will contain estimate_info. This entry contains a convergence flag (conv_flag, logical) and the optimizer (a string).

Finally, if arbitrary R objects were passed via ddm_opts() (to access these objects when evaluating the component functions) the list will contain an entry ddm_opts.

Every model also has the attribute b_coding, which summarizes how the boundaries are labeled.

For print.drift_dm(), the supplied drift_dm object x (invisible return).

See Also

conds(), flex_prms(), prms_solve(), solver(), obs_data(), comp_funs(), b_coding(), coef(), pdfs()

Examples

# Plain call, with default component functions -----------------------------
# create parameter and condition vectors
prms <- c(muc = 4, b = 0.5)
conds <- c("one", "two")

# then call the backbone function (note that we don't provide any component
# functions, so dRiftDM uses the default functions as documented in
# comp_funs())
my_model <- drift_dm(prms_model = prms, conds = conds, subclass = "example")
print(my_model)

Available Cost Functions for Model Estimation

Description

Returns the names of implemented cost functions.

Usage

drift_dm_cost_functions()

Value

A character vector of cost function names.

Pre-built Drift Diffusion Models

Description

Returns the names of available pre-built DDMs in dRiftDM.

Usage

drift_dm_pre_built_models()

Value

A character vector of model names.

Available types of statistics

Description

Internal helper to return supported statistic types depending on the context (e.g., for observed data.frames or fitted model objects).

Usage

drift_dm_stats_types(context = NULL)

Arguments

Value

a character vector of valid statistic types for the given context.

Maybe droplevels of ID column

Description

This function takes a data frame with an ID colmumn, and drops the unused levels from the ID column if it is factor; in this case a warning is thrown

Usage

drop_levels_ID_column(some_data)

Arguments

Value

if the ID column is not of type factor, then the unmodified object is returned.

if the ID column is of type factor, droplevels is applied, and if levels were dropped, a warning is thrown

Derivative of a Constant Boundary

Description

Derivative of a Constant Boundary

Usage

dt_b_constant(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

a vector of the same length as t_vec only zeros.

Truncated Normal Density Function

Description

Computes the probability density function for the truncated normal distribution. This version supports both vector and matrix input for x.

Usage

dtnorm(x, mean = 0, sd = 1, lower = -Inf, upper = Inf, log = FALSE)

rtnorm(n, mean = 0, sd = 1, lower = -Inf, upper = Inf)

Arguments

Details

The function evaluates the normal density at x and scales it to reflect truncation to the interval (lower, upper). Values outside the truncation bounds are assigned a density of 0 (or -Inf on the log scale). Internally, stats::dnorm and stats::pnorm are used.

If x is a matrix, the result retains the same dimensions. All other arguments are recycled as needed. For example, if x has two rows and 5 columns, then mean might provide 2 values, so that the first/second row is evaluated under the first/second mean value.

Value

A numeric vector or matrix of the same shape as x, containing the (possibly log) densities.

(Hierarchical) Bayesian Estimation with Differential Evolution

Description

Estimate group-level and individual-level parameters with a hierarchical Bayesian approach using Differential Evolution MCMC (DE-MCMC) Turner et al. (2013). An approximation of the marginal likelihood to calculate Bayes Factors can be obtained with the Thermodynamic Integration via Differential Evolution (TIDE) sampler Evans and Annis (2019).

Usage

estimate_bayes_h(
  drift_dm_obj,
  obs_data_ids,
  sampler,
  n_chains,
  burn_in,
  samples,
  n_cores,
  prob_migration,
  prob_re_eval,
  verbose,
  seed = NULL,
  ...
)

estimate_bayes_one_subj(
  drift_dm_obj,
  sampler,
  n_chains,
  burn_in,
  samples,
  prob_migration,
  prob_re_eval,
  verbose,
  ...
)

Arguments

Details

The function estimate_bayes_h() handles the hierarchical case. The function estimate_bayes_one_subj() handles the case for estimating a single individual. The reason for writing two functions is that the hierarchical case has some unique tweaks to it that need to be considered ... and writing one function would be quite the mess.

Prior Settings: See the wrapper estimate_bayesian() and also get_default_prior_settings()

Value

A named ist containing posterior samples for group-level and individual-level parameters, log-posterior values, and log-likelihoods. Labels: phi, pis_phi, lls_phi, theta, pis_theta, lls_theta. The first three entries are only present in the hierarchical case.

The list also has an additional attribute named data_model. In the hierarchical case, the attribute contains a named list of model copies with all the individual data sets attached. The list is named according to the individual IDs in the argument obs_data_ids. In the non-hierarchical case, the attribute contains the model and its attached data.

References

Turner BM, Sederberg PB, Brown SD, Steyvers M (2013). “A method for efficiently sampling from distributions with correlated dimensions.” Psychological Methods, 18(3), 368–384. doi:10.1037/a0032222. Evans NJ, Annis J (2019). “Thermodynamic integration via differential evolution: A method for estimating marginal likelihoods.” Behavior Research Methods, 51, 930–947. doi:10.3758/s13428-018-1172-y.

(Hierarchical) Bayesian Estimation

Description

This function provides a wrapper around the implemented samplers for Bayesian inference in dRiftDM. For parameter estimation, Differential Evolution Markov-Chain Monte-Carlo (DE-MCMC) Turner et al. (2013) is used. An approximation of the marginal likelihood to calculate Bayes Factors can be obtained with the Thermodynamic Integration via Differential Evolution (TIDE) sampler Evans and Annis (2019). However, TIDE is not yet supported fully, and is at an experimental stage.

Usage

estimate_bayesian(
  drift_dm_obj,
  obs_data_ids = NULL,
  sampler,
  n_chains,
  burn_in,
  samples,
  prob_migration,
  prob_re_eval,
  verbose = NULL,
  ...
)

Arguments

Details

When a data.frame is supplied, a hierarchical approach to parameter estimation is done. In this case, the supplied data set must provide data for multiple individuals. To estimate the parameters for a single individual (i.e., pursue the non-hierarchical approach), then the supplied model drift_dm_obj must have data attached to it (see obs_data()).

Value

an object of type mcmc_dm containing posterior samples for parameters, log-posterior values, and log-likelihoods. In the hierarchical case, the respective values are available at both the group-level and the individual-level. The object contains two attributes: sampler and data_model. The former simply stores the type of sampler that was used and codes whether estimation was done in a hierarchical fashion or not. The latter either contains the model and the attached data (in the non-hierarchical case) or a named list of model copies with each individual's data attached.

References

Turner BM, Sederberg PB, Brown SD, Steyvers M (2013). “A method for efficiently sampling from distributions with correlated dimensions.” Psychological Methods, 18(3), 368–384. doi:10.1037/a0032222. Evans NJ, Annis J (2019). “Thermodynamic integration via differential evolution: A method for estimating marginal likelihoods.” Behavior Research Methods, 51, 930–947. doi:10.3758/s13428-018-1172-y.

See Also

summary.mcmc_dm(), estimate_bayes_h(), estimate_bayes_one_subj()

Estimate Parameters of a drift_dm Model via Classical Optimization

Description

estimate_classical() estimates the parameters of a drift_dm model by minimizing the model's cost function (e.g., RMSE or negative log-likelihood) using classical (non-Bayesian) optimization routines.

Available optimizers include:

• Nelder-Mead (bounded or unbounded): "Nelder-Mead", "nmkb" (via stats::optim() and dfoptim::nmkb(), respectively)

• BFGS and L-BFGS-B (via stats::optim())

• Differential Evolution (via DEoptim::DEoptim())

Usage

estimate_classical(
  drift_dm_obj,
  optimizer,
  start_vals = NULL,
  return_runs = NULL,
  lower = NULL,
  upper = NULL,
  verbose = NULL,
  de_n_cores = 1,
  control = list(),
  round_digits = NULL,
  seed = NULL,
  use_ez = NULL,
  n_lhs = NULL
)

Arguments

Details

Search space specification

lower and upper can be specified flexibly:

• As unnamed numeric vectors (not recommended unless you're sure of the parameter order)

• As named numeric vectors matching the parameter names of the model

• As lists with a default_values entry (plus optional condition-specific entries)

This design mirrors the structure used in simulate_data.drift_dm().

Optimization details

Some optimizers (i.e., "nmkb", "L-BFGS-B", "DEoptim") require both lower and upper bounds.

Differential Evolution (DEoptim) supports parallelization across cores via de_n_cores. If de_n_cores > 1, a parallel cluster is created and automatically closed after optimization.

The cost function being minimized depends on the cost_function() of the model.

During optimization, failed model evaluations yield a very high cost value (i.e., .Machine$double.xmax). In some cases, this ensures that the optimization doesn't crash, though, this is not guaranteed.

Value

The updated drift_dm_obj, with optimized parameters.

Internal wrapper for classical estimation of individuals

Description

This function wraps estimate_classical() to handle fitting multiple individuals in a consistent way. It prepares the data, distributes the estimation across individuals, and manages parallelization strategies and progress reporting. Unlike the deprecated estimate_model_ids(), this function no longer saves results to disk — instead, it directly returns an object of class fits_ids_dm.

Usage

estimate_classical_wrapper(
  drift_dm_obj,
  obs_data_ids,
  parallelization_strategy = NULL,
  progress = NULL,
  start_vals = NULL,
  optimizer,
  n_cores = NULL,
  seed = NULL,
  ...
)

Arguments

Details

Convergence issues are checked automatically. If one or more individuals fail to converge, a warning is issued with the corresponding IDs and messages returned by the optimizer.

Value

an object of class fits_ids_dm, which is a list with two components:

• drift_dm_fit_info — a list containing the model object, observed data, optimizer information, and convergence messages

• all_fits — a list of individual estimation results

See Also

estimate_classical(), estimate_dm()

Fit a DDM to Observed Data

Description

estimate_dm() is the main function to fit a drift diffusion model (DDM) in dRiftDM. Several ways of fitting a model are supported: fitting a single participant, fitting multiple participants separately or aggregated, and fitting a (hierarchical) Bayesian model. The particular way is controlled via the approach argument.

Usage

estimate_dm(
  drift_dm_obj,
  obs_data = NULL,
  approach = NULL,
  optimizer = NULL,
  control = list(),
  n_cores = 1,
  parallelization_strategy = NULL,
  lower = NULL,
  upper = NULL,
  start_vals = NULL,
  means = NULL,
  sds = NULL,
  shapes = NULL,
  rates = NULL,
  n_chains = 40,
  burn_in = 500,
  samples = 1000,
  prob_migration = 0.1,
  prob_re_eval = 1,
  messaging = TRUE,
  seed = NULL,
  ...
)

## S3 method for class 'fits_agg_dm'
print(x, ...)

## S3 method for class 'fits_ids_dm'
print(x, ...)

## S3 method for class 'mcmc_dm'
print(x, ..., round_digits = drift_dm_default_rounding())

Arguments

Details

Fitting Approaches

The function supports different "approaches" to fitting data.

• "sep_c": This means that data is always considered separately for each participant (if there are multiple participants) and that a classical approach to parameter optimization is used. This means that a standard cost_function is minimized (e.g., the negative log-likelihood). If users provide only a single participant or a data set without an ID column, then the model is fitted just once to that data set.

• "agg_c": This fits the model to aggregated data. For each individual in a data set, summary statistics (e.g., quantiles, accuracies) are calculated, and the model is fitted once to the average of these summary statistics.

• "sep_b": Similar to ⁠sep_b"⁠, although a Bayesian approach is used to sample from the posterior distribution.

• "hier_b": A hierarchical approach to parameter estimation. In this case all participants are considered simultaneously and samples are drawn both at the individual-level and group-level.

The optimizers "nmkb", "L-BFGS-B", and "DEoptim" (for classical parameter optimization) require the specification of the lower/upper arguments.

Fitting to Aggregated Data

For aggregated fits, aggregated statistics are set to the model and the cost function is switched to "rmse". If incompatible settings are requested, the function switches to a compatible configuration and informs the user with messages (these messages can be suppressed via the messaging argument).

Specifying lower/upper for Classical optimization

the function estimate_model_dm() provides a flexible way of specifying the optimization space; this is identical to specifying the parameter simulation space in simulate_data.drift_dm().

Users have three options to specify the search space (see also the examples below):

• Plain numeric vectors (not very much recommended). In this case, lower/upper must be sorted in accordance with the parameters in the underlying flex_prms object of drift_dm_obj that vary for at least one condition (call print(drift_dm_obj) and have a look at the columns of the ⁠Parameter Settings⁠ output; for each column that has a number > 0, specify an entry in lower/upper).

• Named numeric vectors. In this case lower/upper have to provide labels in accordance with the parameters that are considered "free" at least once across conditions (call coef(drift_dm_obj) and provide one named entry for each parameter; dRiftDM will try to recycle parameter values across conditions).

• The most precise way is when lower/upper are lists. In this case, the list requires an entry called "default_values" which specifies the named or plain numeric vectors as above. If the list only contains this entry, then the behavior is as if lower/upper were already numeric vectors. However, the lower/upper lists can also provide entries labeled as specific conditions, which contain named (!) numeric vectors with parameter labels. This will modify the value for the upper/lower parameter space with respect to the specified parameters in the respective condition.

Specifying Priors for Bayesian Estimation

(Default) Prior settings in the non-hierarchical case:

Let \theta^{(j)} indicate parameter j of a model (e.g., the drift rate). The prior on \theta^{(j)} is a truncated normal distribution:

\theta^{(j)} \sim NT(\mu^{(j)}, \sigma^{(j)}, l^{(j)}, u^{(j)})

With \mu^{(j)} and \sigma^{(j)} representing the mean and standard deviation of parameter j. l^{(j)} and u^{(j)} represent the lower and upper boundary. \mu^{(j)} is taken from the mean argument or the currently set model parameters (i.e., from coef(drift_dm_obj)) when calling the function. \sigma^{(j)} is, per default, equal to \mu^{(j)}. This can be changed by passing the sd argument. The lower and upper boundaries of the truncated normal are -Inf and Inf per default. This can be altered by passing the arguments lower and upper (see the examples below).

(Default) Prior settings in the hierarchical case:

Let \theta_i^{(j)} indicate parameter j for participant i (e.g., the drift rate estimated for individual i). The prior on \theta_i^{(j)} is a truncated normal distribution:

\theta_i^{(j)} \sim NT(\mu^{(j)}, \sigma^{(j)}, l^{(j)}, u^{(j)})

With \mu^{(j)} and \sigma^{(j)} representing the mean and standard deviation of parameter j at the group level. l^{(j)} and u^{(j)} represent the lower and upper boundary. The lower and upper boundaries of the truncated normal are -Inf and Inf per default. This can be altered by passing the arguments lower and upper.

For a group-level mean parameter, \mu^{(j)}, the prior is also a truncated normal distributions:

\mu^{(j)} \sim NT(M^{(j)}, SD^{(j)}, l^{(j)}, u^{(j)})

With M^{(j)} specified by the mean argument or the currently set model parameters. SD^{(j)} is, per default, equal to M^{(j)}. This can be changed by passing the sd argument.

For a group-level standard deviation parameter, \sigma^{(j)}, the prior is a gamma distribution:

\sigma^{(j)} \sim \Gamma(shape^{(j)},rate^{(j)})

With shape^{(j)} and rate^{(j)} being 1 by default. This can be changed by passing the arguments shape and rate.

Specifying Prior Settings/Arguments

Argument specification for mean, sd, lower, upper, shape and rate is conceptually identical to specifying lower/upper for the classical optimization approach (see the subsection above and the examples below).

Value

• If fitting a single individual: either a drift_dm object with fitted parameters and additional fit information (for the classical optimization framework) or an object of type mcmc_dm (for the Bayesian framework)

• If fitting multiple individuals separately: a fits_ids_dm object or a list of mcmc_dm objects, containing all the individual model fits.

• If fitting aggregated data: a fits_agg_dm object containing the model itself and the raw data.

• If fitting multiple individuals hierarchically: an object of type mcmc_dm.

Note

estimate_dm dispatches to underlying estimation routines that are not exported:

• Classical optimization of one individual via estimate_classical()

• Classical optimization of multiple individuals via estimate_classical_wrapper()

• Bayesian estimation via estimate_bayesian().

• Aggregated fitting is handled within estimate_dm() in combination with estimate_classical()

When fitting a model with optimizer = "DEoptim", the corresponding minimization routine always runs for 200 iterations by default, irrespective of whether a minimum has already been reached (see DEoptim::DEoptim.control). Therefore, with default optimization settings, estimate_dm() returns the convergence flag NA for optimizer = "DEoptim", because the termination of the routine does not necessarily indicate convergence. However, this is typically not an issue, as 200 iterations are generally sufficient for the algorithm to find the global minimum. If users explicitly define convergence criteria via the control argument of estimate_dm() (which is passed on to DEoptim::DEoptim.control), valid convergence messages and flags are returned.

See Also

estimate_classical(), estimate_bayesian(), estimate_classical_wrapper(), get_parameters_smart()

Examples

##########
# Note: The following examples were trimmed for speed to ensure they run
# within seconds. They do not always provide realistic settings.
##########

####
# Setup

# get a model for the examples (DMC with just two free parameters)
model <- dmc_dm(
  instr = "
   b <!>
   non_dec <!>
   sd_non_dec <!>
   tau <!>
   alpha <!>
   "
)

# get some data (the first two participants in the data set of Ulrich et al.)
data <- ulrich_flanker_data[ulrich_flanker_data$ID %in% 1:2, ]


####
# Fit a single individual (using unbounded Nelder-Mead)
fit <- estimate_dm(
  drift_dm_obj = model,
  obs_data = data[data$ID == 1, ],
  optimizer = "Nelder-Mead"
)
print(fit)




####
# Fit a single individual (using DEoptim)
l_u <- get_lower_upper(model)
set.seed(2)
fit <- estimate_dm(
  drift_dm_obj = model,
  obs_data = data[data$ID == 1, ],
  optimizer = "DEoptim",
  lower = l_u$lower, upper = l_u$upper,
  control = list(itermax = 5) # use default itermax in practice!
)
print(fit)


####
# Fit multiple individuals (separately; using bounded Nelder-Mead)
l_u <- get_lower_upper(model)
fit <- estimate_dm(
  drift_dm_obj = model,
  obs_data = data, # contains the data for two individuals
  optimizer = "nmkb",
  lower = l_u$lower, upper = l_u$upper,
)
print(fit)
coef(fit)


###
# Fit to aggregated data (using unbounded Nelder-Mead)
fit <- estimate_dm(
  drift_dm_obj = model,
  obs_data = data, # contains data for two individuals
  optimizer = "Nelder-Mead",
  approach = "agg_c"
)
print(fit)
coef(fit)


###
# EXPERIMENTAL
# Fit a single individual (using DE-MCMC; Bayesian; custom priors)
fit <- estimate_dm(
  drift_dm_obj = model,
  obs_data = data[data$ID == 1, ],
  approach = "sep_b",
  burn_in = 1, # higher in practice (e.g., 500)
  samples = 1, # higher in practice (e.g., 1000)
  n_chains = 5, # higher in practice (e.g., 40)
  mean = c(muc = 3, A = 0.9),
  sd = c(muc = 2, A = 0.8),
)
print(fit)
coef(fit)


###
# EXPERIMENTAL
# Fit multiple individuals (using DE-MCMC; hierarchical Bayesian)
fit <- estimate_dm(
  drift_dm_obj = model,
  approach = "hier_b",
  obs_data = data, # contains data for two individuals
  burn_in = 1, # higher in practice (e.g., 500)
  samples = 1, # higher in practice (e.g., 1000)
  n_chains = 5, # higher in practice (e.g., 40)
  n_cores = 1, # higher in practice (depending on your machine and data set)
)
print(fit)
coef(fit)

Estimate the Parameters of a drift_dm Model

Description

This function was deprecated in dRiftDM version v.0.3.0, please use the more general estimate_dm() function.

Old documentation: Find the 'best' parameter settings by fitting a drift_dm models' predicted probability density functions (PDFs) to the observed data stored within the respective object. The fitting procedure is done by minimizing the negative log-likelihood of the model.

Users have three options:

• Estimate the parameters via Differential Evolution (Default)

• Estimate the parameters via (bounded) Nelder-Mead

• Use Differential Evolution followed by Nelder-Mead.

See also vignette("dRiftDM", "dRiftDM")

Usage

estimate_model(
  drift_dm_obj,
  lower,
  upper,
  verbose = 0,
  use_de_optim = TRUE,
  use_nmkb = FALSE,
  seed = NULL,
  de_n_cores = 1,
  de_control = list(reltol = 1e-08, steptol = 50, itermax = 200, trace = FALSE),
  nmkb_control = list(tol = 1e-06)
)

Arguments

Details

Specifying lower/upper

the function estimate_model provides a flexible way of specifying the search space; identical to specifying the parameter simulation space in simulate_data.drift_dm.

Users have three options to specify the simulation space:

• Plain numeric vectors (not very much recommended). In this case, lower/upper must be sorted in accordance with the parameters in the flex_prms_obj object that vary for at least one condition (call print(drift_dm_obj) and have a look at the ⁠Parameter Settings⁠ output)

• Named numeric vectors. In this case lower/upper have to provide labels in accordance with the parameters that are considered "free" at least once across conditions.

• The most flexible way is when lower/upper are lists. In this case, the list requires an entry called "default_values" which specifies the named or plain numeric vectors as above. If the list only contains this entry, then the behavior is as if lower/upper were already numeric vectors. However, the lower/upper lists can also provide entries labeled as specific conditions, which contain named (!) numeric vectors with parameter labels. This will modify the value for the upper/lower parameter space with respect to the specified parameters in the respective condition.

Details on Nelder-Mead and Differential Evolution

If both use_de_optim and use_nmkb are TRUE, then Nelder-Mead follows Differential Evolution. Note that Nelder-Mead requires a set of starting parameters for which either the parameter values of drift_dm_obj or the estimated parameter values by Differential Evolution are used.

Default settings will lead DEoptim::DEoptim to stop if the algorithm is unable to reduce the negative log-likelihood by a factor of reltol * (abs(val) + reltol)after steptol = 50 steps, with reltol = 1e-8 (or if the default itermax of 200 steps is reached). Similarly, dfoptim::nmkb will stop if the absolute difference of the log-likelihood between successive iterations is below tol = 1e-6.See DEoptim::DEoptim.control and the details of dfoptim::nmkb for further information.

Value

the updated drift_dm_obj (with the estimated parameter values, log-likelihood, and probability density functions of the first passage time)

See Also

estimate_model_ids

Fit Multiple Individuals and Save Results

Description

This function was deprecated in dRiftDM version 0.3.0. Please use the more general estimate_dm() instead. NOTE: dRiftDM now supports multiple ways of estimating a model. To ensure a more consistent function interface, individual fits are no longer saved to disk when fitting multiple participants. Instead, estimate_dm() directly returns an object of type fits_ids_dm, which users can save manually if desired.

Old documentation: Provides a wrapper around estimate_model to fit multiple individuals. Each individual will be stored in a folder. This folder will also contain a file drift_dm_fit_info.rds, containing the main arguments of the function call. One call to this function is considered a "fit procedure". Fit procedures can be loaded via load_fits_ids.

Usage

estimate_model_ids(
  drift_dm_obj,
  obs_data_ids,
  lower,
  upper,
  fit_procedure_name,
  fit_path,
  fit_dir = "drift_dm_fits",
  folder_name = fit_procedure_name,
  seed = NULL,
  force_refit = FALSE,
  progress = 2,
  start_vals = NULL,
  ...
)

Arguments

Details

Examples and more information can also be found in vignette("dRiftDM", "dRiftDM").

When developing the fitting routine we had three levels of files/folders in mind:

• In a directory/folder named fit_dir multiple fitting routines can be stored (default is "drift_dm_fits")

• Each fitting routine has its own folder with a name as given by folder_name (e.g., "ulrich_flanker", "ulrich_simon", ...)

• Within each folder, a file called drift_dm_fit_info.rds contains the main information about the function call. That is, the time when last modifying/calling a fitting routine, the lower and upper parameter boundaries, the drift_dm_object that was fitted to each individual, the original data set obs_data_ids, and the identifier fit_procedure_name. In the same folder each individual has its own ⁠<individual>.rds⁠ file containing the modified drift_dm_object.

Value

nothing (NULL; invisibly)

See Also

load_fits_ids

Specify custom parameters

Description

This function takes a flex_prms_obj and adds or builds the entry cust_prms to allow for custom parameters. An examplary instruction is "peak_l = (a-1)*tau"

Usage

flex_cust_prm(flex_prms_obj, formula_instr)

Arguments

Details

The entry cust_prms is a list with entries expressions and values. Each of these is again a named list, that either contains the expression with instructions on how to calculate the custom parameter (e.g., "peak_l") or the respective values. Values are getting updated/calculated in update_special_values()

The cust_prms exists of two entries "expressions" and values". "expressions" contains a named list, with expressions referring to prms_matrix (see flex_prms) on how to calculate the custom parameter (across all conditions). The "values" contain a named list, with named numeric vectors (names are conditions, values the calculated custom parameter values)

Value

a modified flex_prms object with respect to the cust_prms entry

Exclude parameters from being modified (i.e., fix it; internal docu)

Description

This function modifies the (linear) internal list and sets the desired parameters (based on the instruction string) to 0. This indicates that this parameter is not altered within the function x2prms_vals (i.e., ' prm <!> conda')

Usage

flex_fix_prms(flex_prms_obj, formula_instr)

Arguments

Value

a modified flex_prms_obj with respect to the (linear) internal list

Flex_Prms

Description

Functions for creating, accessing replacing, or printing a flex_prms object. Any object of type flex_prms provides a user-friendly way to specify dependencies, parameter values etc. for a model.

Usage

flex_prms(object, ...) <- value

## S3 replacement method for class 'drift_dm'
flex_prms(object, ..., eval_model = FALSE) <- value

flex_prms(object, ...)

## S3 method for class 'numeric'
flex_prms(object, ..., conds, instr = NULL, messaging = NULL)

## S3 method for class 'flex_prms'
flex_prms(object, ...)

## S3 method for class 'drift_dm'
flex_prms(object, ...)

## S3 method for class 'flex_prms'
print(
  x,
  ...,
  round_digits = drift_dm_default_rounding(),
  dependencies = TRUE,
  cust_parameters = TRUE
)

Arguments

Details

Objects of type flex_prms can be modified using the generic modify_flex_prms() function and a corresponding set of "instructions" (see the respective function for more details).

flex_prms() is a generic function. If called with a named numeric vector, then this will create an object of type flex_prms (requires conds to be specified). If called with other data types, gives the respective flex_prms object

⁠flex_prms<-()⁠ is a generic replacement function. Currently this only supports objects of type drift_dm. It will replace/update the model with a new instance of type flex_prms.

Value

The specific value returned depends on which method is called

Creating an object of type flex_prms

Can be achieved by calling flex_prms() with a named numeric vector, thus when calling the underlying method flex_prms.numeric (see the example below). In this case a list with the class label "flex_prms" is returned. It containts three entries:

• A nested list internal_list. This list specifies the dependencies and restrains enforced upon the parameters across conditions. Integers >= 1 indicate that this parameter will be estimated for a specific condition, and conditions with the same number refer to a single parameter. Integers == 0 indicate thtat this parameter will not be esitmated for a specific condition (i.e., it is considered "fixed"). Expressions will be evaluated at run time and specify special dependencies among parameters.

• A nested list linear_internal_list. This list essentially contains the same information as internal_list, but the parameters are sorted so that they can be mapped to an integer vector (relevant only in the depths of the package for the minimization routines).

• A numeric matrix prms_matrix which contains the currently set values for each parameter across all conditions. Per default, the values of each parameter are set equal across all conditions. Additionally, each parameter is assumed to be restrained as equal across all conditions. The values for all parameters given a condition will be passed to the component functions (see comp_funs).

• (optional) A list of additional parameters cust_prms that are derived from the parameters in prms_matrix.

Accessing an object of type flex_prms

Users can access/get the flex_prms object when calling flex_prms() with an object of type drift_dm, fits_ids_dm (see estimate_model_ids()), or flex_prms. In this case, the stored flex_prms object is returned.

Replacing an object of type flex_prms

The flex_prms object stored within an object of type drift_dm can be replaced by calling the generic ⁠flex_prms<-⁠ replacement function. In this case, the modified drift_dm object is returned.

Printing an object of type flex_prms

The print.flex_prms() method invisibly returns the supplied flex_prms object.

Note

There is only a replacement function for drift_dm objects. This is because replacing the solver settings after the model has been fitted (i.e., for a fits_ids_dm object) doesn't make sense.

See Also

estimate_model_ids(), drift_dm(), summary.flex_prms(), modify_flex_prms()

Examples

# Create a flex_prms object -----------------------------------------------
conds <- c("one", "two")
prms <- c(muc = 3, b = 0.5)
one_instr <- "muc ~ one + two"
flex_prms_obj <- flex_prms(
  prms,
  conds = conds,
  instr = one_instr
)
print(flex_prms_obj)


# Access a flex_prms object of a model ------------------------------------
my_model <- ratcliff_dm() # the Ratcliff DDM comes with dRiftDM
print(flex_prms(my_model))


# Replace the flex_prms object of a model ---------------------------------
# create a new flex_prms object
conds <- c("one", "two")
prms <- c(muc = 3, b = 0.6, non_dec = 0.3)
new_flex_prms_obj <- flex_prms(
  prms,
  conds = conds
)

flex_prms(my_model) <- new_flex_prms_obj

# acess the new flex_prms object
print(flex_prms(my_model))


# Control the print method -------------------------------------------------
dmc_model <- dmc_dm() # another, more complex, model; comes with dRiftDM
print(flex_prms(dmc_model), round_digits = 1, cust_parameters = FALSE)

Set parameters as equal across conditions

Description

This function takes a flex_prms object and modifies the (linear) internal list so that a parameter is set as equal across multiple conditions, according to the instruction formula (i.e., ' prm ~! conda + condb)

Usage

flex_restrain_prms(flex_prms_obj, formula_instr)

Arguments

Value

a modified flex_prms object with an updated (linear) internal list

Set special dependencies (internal docu)

Description

Sets special dependencies so that parameters depend on other parameters. (i.e., 'prmX ~ conda == -(prmY ~ condb)')

Usage

flex_special_dependency(flex_prms_obj, formula_instr)

Arguments

Value

a modified flex_prms_object with a modified (linear) internal list and modified parameter and custom parameter matrices

Set a specific value to the parameter matrix (internal docu)

Description

This function takes a flex_prms_obj and sets certain values to the parameter matrix, based on the given instruction string (i.e., ' prm ~ conda => 0.3)

Usage

flex_specific_value(flex_prms_obj, formula_instr)

Arguments

Value

an updated flex_prms object with a modified prms_matrix object, and (if applicable) a modified cust_prms matrix

Allow parameters to vary

Description

This function takes an object of type flex_prms and a instruction string that refers to a "vary" instruction (i.e., ab ~ cd + ds). This string is broken down and unique parameters are introduced for the condition x parameter combinations

Usage

flex_vary_prms(flex_prms_obj, formula_instr)

Arguments

Value

an updated flex_prms_obj with an updated (linear) internal list

Perform a Full Crossover Step Using Differential Evolution

Description

This function updates each chain's parameters by proposing new values using a differential evolution strategy. For each chain k, two other chains m and n are randomly selected, and a proposal is generated via: prms_k + gamma * (prms_m - prms_n) + noise, where gamma = 2.38 / sqrt(2 * n_prms) and noise is uniform perturbation controlled by b. The proposal is accepted with Metropolis probability via call_log_posterior_m(), and accepted proposals replace the current values.

Usage

full_crossover(
  prms_across_chains,
  pis_across_chains,
  log_likes_across_chains,
  gamma = NULL,
  b = 0.001,
  ...
)

Arguments

Value

A list with the following components:

• new_prms_across_chains: The updated parameter matrix of shape ⁠p × n⁠.

• new_pis_across_chains: The updated vector of log-posterior values.

• new_log_likes_across_chains: The updated vector of log-likelihood values.

Compute average trials per condition across individuals

Description

internal helper — assumes each ID appears with the same set of conditions in column Cond. counts trials per ID×Cond and averages across ID.

Usage

get_avg_trials(obs_data_ids)

Arguments

Value

list with:

• N: number of unique individuals

• avg_trials: named numeric vector of average trials per condition

Get default/fall back component functions

Description

If arguments are provided that are not NULL, the respective argument is simply returned. If it is NULL, then a default/fall back component function is returned for the respective component. This function is called to fill up non-specified component functions when calling drift_dm().

Usage

get_default_functions(
  mu_fun = NULL,
  mu_int_fun = NULL,
  x_fun = NULL,
  b_fun = NULL,
  dt_b_fun = NULL,
  nt_fun = NULL
)

Arguments

Details

defaults...

• mu_fun -> constant drift rate of 3 (i.e., vector of 0s)

• mu_int_fun -> constant drift rate of 3 (i.e., vector of 3 times t_vec)

• x_fun -> dirac delta on zero x_dirac_0()

• b_fun -> constant boundary of 0.5 (i.e., vector of 0.5s)

• dt_b_fun -> derivate of constant boundary (i.e., vector of 0s). dt_b_constant()

• nt_fun -> constant non-decision time of 0.3 (i.e., vector for dirac delta on 0.5).

Value

a list of mu_fun, mu_int_fun, x_fun, b_fun, dt_b_fun, and nt_fun, with either the supplied component functions or the added/filled in default component functions (if an argument is NULL).

Generate Prior Functions for Model Parameters

Description

This function creates prior distribution functions for each model parameter in a drift diffusion model (DDM), depending on the specified hierarchical level. It returns both log-density functions and, where applicable, random-sample generators based on the user-defined prior settings.

Usage

get_default_prior_settings(
  drift_dm_obj,
  level,
  means = NULL,
  sds = NULL,
  lower = NULL,
  upper = NULL,
  shapes = NULL,
  rates = NULL
)

Arguments

Details

Each prior is parameter-specific and wrapped using purrr::partial() so that downstream sampling or density evaluation can be performed easily. At the hyper-level, the functions d_default_prior_hyper() and r_default_prior_hyper() are used. At the lower-level, the functions dtnorm() and rtnorm() are used.

The input arguments means, sds, lowers, uppers, shapes, and rates are handled by the function get_parameters_smart().

Value

A named list with two elements:

• log_dens_priors: A named list of functions. Each function returns the log-density for a parameter value, based on the chosen prior settings.

• r_priors: A named list of functions for sampling from the specified prior distributions.

See Also

get_parameters_smart(), dtnorm(), rtnorm(), d_default_prior_hyper(), r_default_prior_hyper()

Auxiliary Function to load a fits_ids_dm, fits_agg_dm, or mcmc_dm object

Description

The function is merely helper functions to create an object of type fits_ids_dm, fits_agg_dm, or mcmc_dm. It is used for example code.

Usage

get_example_fits(class, hierarchical = FALSE)

Arguments

Details

For "fits_ids_dm", the returned object comprises DMC (see dmc_dm()) fitted to three participants of the ulrich_flanker_data.

For "fits_agg_dm", the returned object comprises the Ratcliff model (see ratcliff_dm()) fitted to synthetic data of three participants.

For "mcmc_dm" and hierarchical = FALSE, the returned object comprises the Ratcliff model (see ratcliff_dm()) fitted to synthetic data of one participant.

For "mcmc_dm" and hierarchical = TRUE, the returned object comprises the Ratcliff model (see ratcliff_dm()) fitted to synthetic data of ten participants.

Value

An object of type fits_ids_dm, fits_agg_dm, or mcmc_dm, mimicking a result from calling estimate_dm().

Examples

get_example_fits(class = "fits_agg")

Compute EZ Diffusion parameters

Description

Internal helper that computes EZ diffusion model parameters for each condition in a drift_dm_obj. The computation is based on the equations from Wagenmakers et al. (2007) and estimates drift rate (muc), boundary separation (b), and non-decision time (non_dec).

Usage

get_ez_diffusion(drift_dm_obj)

Arguments

Details

If Pc equals 0, 0.5, or 1, small adjustments are applied to prevent numerical issues in the logit transformation.

Value

a matrix with rows muc, b, and non_dec

Get Default Parameter Ranges for a Model

Description

get_lower_upper() returns suggested default values for parameter bounds of a drift_dm model. The function inspects the model's component functions (e.g., drift, boundary, non-decision time, start) and provides heuristic defaults for some of the pre-built components. Only parameters that are currently considered free in the model are returned.

Usage

get_lower_upper(object, ...)

## S3 method for class 'drift_dm'
get_lower_upper(object, ..., warn = TRUE)

Arguments

Details

Supported components include: mu_constant, mu_dmc, mu_ssp, b_constant, x_uniform, x_beta, nt_constant, nt_uniform, nt_truncated_normal. For some defaults we use the model's discretization (dt, dx) to ensure sensible minima.

If a component is not recognized (or refers to currently unsupported components), no defaults are provided for that component. When warn = TRUE, a single warning lists components without defaults and any free parameters that remain unmatched. In this case, the user has to add the missing parameter ranges before attempting to fit the model.

The default ranges are heuristics intended to provide a reasonable starting point for new users. They are not guaranteed to be appropriate for every model or data set. Always review and, if needed, adjust the returned values as needed.

Value

a list with two named numeric vectors:

• lower — suggested lower bounds for free parameters

• upper — suggested upper bounds for free parameters

Examples

# get a model for the example
model <- dmc_dm(obs_data = dmc_synth_data)

# get the parameter ranges
lu <- get_lower_upper(model)
lu$lower
lu$upper

# then continue to estimate
# estimate_dm(model, lower = lu$lower, upper = lu$upper, optimizer = "nmkb")

Turn default/special parameter specifications to vectors

Description

The function is used in the depths to map parameter inputs to the parameters of a model. One application is to get the search space as a vector, matching with the free parameters of a model. Other applications map, for example, mean values to the free parameters of a model. Relevant when users use the "default parameters" approach where they only specify the parameter labels and assume the package figures out how each parameter relates across conditions (see simulate_data). This comes in handy, when freeing a parameter across conditions, while the search space remains the same (otherwise, a user would always have to adapt the vectors for lower/upper to match with x2prms_vals)

Usage

get_parameters_smart(
  drift_dm_obj,
  input_a,
  input_b = NULL,
  labels = TRUE,
  is_l_u = TRUE,
  fill_up_with = NULL
)

Arguments

Details

The function first gets all unique parameters across conditions using prms_cond_combo. The unique parameter labels are then forwarded to create_matrix_smart, together with all (!) the conditions in the model and the input_a/input_b arguments. Subsequently, the created matrices are wrangled into vectors in accordance with prms_cond_combo. The vectors are then passed back.

Value

a list with two entries named vec_a/vec_b. The length and names (if requested) matches with coef(model, select_unique = TRUE). When input_a and/or input_b is NULL, the respective entry for vec_a/vec_b will be NULL as well.

TITLE

Description

TITLE

Usage

get_starting_values(object, ...)

## S3 method for class 'drift_dm'
get_starting_values(
  object,
  ...,
  lower = NULL,
  upper = NULL,
  verbose = 0,
  use_ez = NULL,
  n_lhs = NULL
)

Arguments

Extract a Subset of MCMC Chains

Description

When calling estimate_bayesian(), the MCMC results are packed up as an mcmc_dm object. This function is used in the depths of dRiftDM to extract the relevant array of MCMC samples, depending on whether the model is hierarchical and whether a participant ID is provided.

Usage

get_subset_chains(chains_obj, id = NULL)

Arguments

Value

A 3D array of MCMC samples. The first dimension indicates parameters, the second dimension chains, and the third dimension iterations

Plot Parameter Distribution(s)

Description

This function creates a histogram for each parameter in a coefs_dm object, resulting from a call to coef.fits_ids_dm.

Usage

## S3 method for class 'coefs_dm'
hist(
  x,
  ...,
  conds = NULL,
  col = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = "value",
  ylab = NULL,
  bundle_plots = TRUE
)

Arguments

Details

The hist.coefs_dm function is designed for visualizing parameter distributions.

If multiple conditions are present, it overlays histograms for each condition with adjustable transparency.

When bundle_plots is set to TRUE, histograms for each parameter are displayed in a grid layout within a single graphics device.

This function has some customization options, but they are limited. If you want to have a highly customized histogram, it is best to create it on your own using R's graphics::hist() function (see the examples below).

Value

Nothing (NULL; invisibly)

Examples

# get an auxiliary fit procedure result (see the function load_fits_ids)
all_fits <- get_example_fits("fits_ids")
coefs <- coef(all_fits)
print(coefs)
hist(coefs, bundle_plots = FALSE) # calls hist.coefs_dm method of dRiftDM

# how to fall back to R's hist() function for heavy customization
coefs <- unpack_obj(coefs) # provides the plain data.frame
hist(coefs$muc, main = expression(mu[c])) # calls graphics::hist()

Aggregate Data Frame Columns by Group

Description

internal function to aggregate columns of a data frame across "ID" while considering a set of grouping columns. It retains the class and attriubtes of the input data.

Usage

internal_aggregate(data, group_cols)

Arguments

Details

internal_aggregate identifies DV columns as those not in group_cols or "ID". It then calculates the mean of these DV columns, grouped by the specified columns. Columns specified in group_cols that are not part of data are ignored silently.

Value

A data.frame containing the aggregated data.

See Also

aggregate_stats(), calc_stats(), new_stats_dm()

Turn an internal list to a matrix

Description

this function wrangles the internal list to character matrix to show how each parameter relates across conditions

Usage

internal_list_to_matrix(internal_list)

Arguments

Value

a character matrix with all parameter labels (columns) and conditions (rows)

Numerical integration using the trapezoidal rule

Description

These are internal helper functions to perform numerical integration via the trapezoidal rule. The workhorse is internal_trapz(), which computes the full integral or returns the cumulative integral.

Usage

internal_trapz(x, y, return_cumsum = FALSE)

cumtrapz(x, y)

trapz(x, y)

Arguments

Details

• internal_trapz(x, y, return_cumsum = FALSE): core implementation

• trapz(x, y) wrapper for internal_trapz(x, y, FALSE), returns the total integral.

• cumtrapz(x, y) wrapper for internal_trapz(x, y, TRUE), returns the cumulative integral.

Value

• trapz(): a single numeric value

• cumtrapz(): numeric vector of cumulative integrals (starting with 0)

• internal_trapz(): either of the above, depending on return_cumsum

Check if an object of length <= 1 is either empty or contains entry with ""

Description

Check if an object of length <= 1 is either empty or contains entry with ""

Usage

is_empty(x)

Arguments

Value

TRUE or FALSE otherwise

Check if an object is a valid numeric vector

Description

This function verifies whether the input is a numeric vector with no missing (NA, NaN) or infinite (Inf or -Inf) values.

Usage

is_numeric(x)

Arguments

Value

A logical value: TRUE if the input is a numeric vector without any missing or infinite values, otherwise FALSE.

Relabel the internal list

Description

The entries of the internal list are either digits (0-x) or expressions. To ensure a valid mapping of these values to an input vector (as done when an optimizer provides input parameters), we have to linearize the list (done whenever modifying the list, see the different flex_* functions )

Usage

linearize_internal_list(internal_list)

Arguments

Value

another list, but with remapped digits in increasing order, while leaving expressions or digits of 0 untouched.

See Also

flex_vary_prms(), flex_restrain_prms(), flex_fix_prms(), flex_special_dependency()

Load Estimates of a Fit Procedure

Description

This function was deprecated in dRiftDM version 0.3.0, because dRiftDM no longer saves model fits to disk when fitting multiple participants. When estimating multiple individuals with the new function estimate_dm(), an object of type fits_ids_dm is returned directly.

Usage

load_fits_ids(
  path = "drift_dm_fits",
  fit_procedure_name = "",
  detailed_info = FALSE,
  check_data = TRUE,
  progress = 2
)

Arguments

Details

Old documentation: This function loads the results of a fit procedure where a model was fitted to multiple individuals (see estimate_model_ids). It is also the function that creates an object of type fits_ids_dm.

with respect to the logic outlined in the details of estimate_model_ids on the organization of fit procedures, path could either point to a directory with (potentially) multiple fit routines or to a specific folder with the individual fits. In either case the intended location is recursively searched for files named drift_dm_fit_info.rds.

If the fit procedure was uniquely located, either because only one fit routine was found in the intended location or because only one drift_dm_fit_info.rds contains the optional identifier specified in fit_procedure_name, then all individual model fits including the information fit_procedure_name are loaded and returned.

In case multiple fit procedures are identified, the user is prompted with a utils::menu, listing information about the possible candidates. The intended fit procedure can then interactively be chosen by the user. The amount of displayed information is controlled via detailed_info.

The print() method for objects of type fits_ids_dm prints out basic information about the fit procedure name, the fitted model, time of (last) call, and the number of individual data sets.

Value

For load_fits_ids(), an object of type fits_ids_dm, which essentially is a list with two entries:

• drift_dm_fit_info, containing a list of the main arguments when estimate_model_ids was originally called, including a time-stamp.

• all_fits, containing a list of all the modified/fitted drift_dm objects. The list's entry are named according to the individuals' identifier (i.e., ID).

For print.fits_ids_dm(), the supplied fit_ids_dm object x (invisible return).

See Also

estimate_model_ids()

Extract Log-Likelihood for a drift_dm Object

Description

This method extracts the log-likelihood for a drift_dm object if possible.

Usage

## S3 method for class 'drift_dm'
logLik(object, ...)

Arguments

Value

A logLik object containing the log-likelihood value for the drift_dm object. This value has attributes for the number of observations (nobs) and the number of model parameters (df).

Returns NULL if the log-likelihood is not available (e.g., when the model has no observed data attached).

Examples

# get a pre-built model and a data set for demonstration purpose
# (when creating the model, set the discretization to reasonable values)
a_model <- dmc_dm()
obs_data(a_model) <- dmc_synth_data
logLik(a_model)

Extract Model Statistics for fits_ids_dm Object

Description

These methods are wrappers to extract specific model fit statistics (log-likelihood, AIC, BIC) for each model in a fits_ids_dm object.

Usage

## S3 method for class 'fits_ids_dm'
logLik(object, ...)

## S3 method for class 'fits_ids_dm'
AIC(object, ..., k = 2)

## S3 method for class 'fits_ids_dm'
BIC(object, ...)

Arguments

Details

Each function retrieves the relevant statistics by calling calc_stats with type = "fit_stats" and selects the columns for ID and the required statistic.

Value

An object of type fit_stats containing the respective statistic in one column (named Log_Like, AIC, or BIC) and a corresponding ID column. If any of the statistics can't be calculated, the function returns NULL.

See Also

stats::AIC(), stats::BIC(), logLik.drift_dm

Examples

# get an auxiliary fits_ids object for demonstration purpose;
# such an object results from calling load_fits_ids
all_fits <- get_example_fits("fits_ids_dm")

# AICs
AIC(all_fits)

# BICs
BIC(all_fits)

# Log-Likelihoods
logLik(all_fits)

# All unique and free parameters
coef(all_fits)

# Or all parameters across all conditions
coef(all_fits, select_unique = FALSE)

Conditional Log-Posterior Distributions for MCMC Sampling

Description

These functions compute conditional log-posterior distributions used in a (hierarchical) MCMC sampler.

Usage

log_posterior_hyper(
  phi_j_mat,
  theta_j_mat,
  log_prior_hyper_fun,
  log_prior_lower_fun,
  temperatures,
  suppress_warnings = TRUE
)

log_posterior_lower(
  thetas_one_subj_mat,
  all_phis_mat,
  model_subj,
  log_prior_lower_funs,
  temperatures,
  suppress_warnings = TRUE
)

Arguments

Details

log_posterior_hyper() computes the conditional log-posterior for a group-level hyperparameter matrix phi_j_mat, given the individual-level parameters across subjects theta_j_mat (for one type of model parameter).

log_posterior_lower() computes the conditional log-posterior for an individual participant’s parameter matrix thetas_one_subj_mat, given prior distributions. In the hierarchical setting, the prior distributions are conditioned on the group-level parameters.

Value

A list with two elements:

• posterior_vals, the total log-posterior values (log-likelihood + log-prior) per chain.

• log_like_vals, the log-likelihood components only, per chain.

Get the maximum number from an internal entry or flex_prms_obj (internal docu)

Description

The entries of the internal list are either digits (0-x) or expressions.

Usage

max_number_one_internal_entry(one_internal_entry)

get_number_prms(flex_prms_obj)

Arguments

Value

the largest digit in the entry or of the linear_list in the supplied flex_prms_obj (0 if there are only expressions). The largest number of the linear_list corresponds to the number of model parameters.

Convert MCMC Chain Array to a coda::mcmc.list Object

Description

Converts a 3D MCMC chain array (parameters × chains × iterations) into a coda::mcmc.list object for compatibility with diagnostic and summary functions from the coda package.

Usage

mcmc_dm_to_coda_mcmc(chains)

Arguments

Value

An object of class mcmc.list containing one mcmc object per chain.

See Also

coda::mcmc(), coda::mcmc.list()

Perform a Migration Step Between Chains

Description

In a migration step, a random subset of chains is selected. Each selected chain k receives a proposal based on the next chain (k + 1) in the sequence (cycling back to the first chain after the last). The proposed parameters are slightly perturbed by uniform noise controlled by b. All proposals are then evaluated simultaneously using Metropolis acceptance probabilities via call_log_posterior_m(), and accepted proposals replace the current values.

Usage

migration_crossover(
  prms_across_chains,
  pis_across_chains,
  log_likes_across_chains,
  b = 0.001,
  ...
)

Arguments

Value

A list with the following components:

• new_prms_across_chains: The updated parameter matrix of shape ⁠p × n⁠.

• new_pis_across_chains: The updated vector of log-posterior values.

• new_log_likes_across_chains: The updated vector of log-likelihood values.

See Also

full_crossover, call_log_posterior_m()

Set Instructions to a flex_prms object

Description

Functions to carry out the "instructions" on how to modify a flex_prms object, specified as a string.

Usage

modify_flex_prms(object, instr, ...)

## S3 method for class 'drift_dm'
modify_flex_prms(object, instr, ..., eval_model = FALSE)

## S3 method for class 'flex_prms'
modify_flex_prms(object, instr, ..., messaging = NULL)

Arguments

Details

modify_flex_prms is a generic function. The default methods pass forward a set of "instructions" to modify the (underlying) flex_prms object.

These instructions are inspired by the model syntax of the lavaan package. Note that specifying multiple instructions is possible, but each instruction has to be defined in its own line. Comments with '#' are possible, also line continuations are possible, if the last symbol is a "+","-", "*", "/", "(", or "[". The following instructions are implemented:

The "vary" instruction:

• Looks something like "a ~ foo + bar"

• This means that the parameter 'a' is allowed to vary independently for the conditions 'foo' and 'bar'

• Thus, when estimating the model, the user will have independent values for 'a' in conditions 'foo' and 'bar'

The "restrain" instruction:

• Looks something like "a ~! foo + bar "

• This means that the parameter 'a' is assumed to be identical for the conditions 'foo' and 'bar'

• Thus, when estimating the model, the user will have only a single value for 'a' in conditions 'foo' and 'bar'

The "set" instruction:

• Users may not always estimate a model directly but rather explore the model behavior. In this case setting the value of a parameter is necessary.

• The corresponding instruction looks something like "a ~ foo => 0.3"

• This will set the value for 'a' in condition 'foo' to the value of 0.3

The "fix" instruction:

• Oftentimes, certain parameters of a model are considered "fixed", so that they don't vary while the remaining parameters are estimated. An example would be the shape parameter 'a' of DMC (see dmc_dm).

• The corresponding instruction looks something like "a <!> foo + bar"

• Usually, users want to call the "set" instruction prior or after the "fix" instruction, to set the corresponding parameter to a certain value.

The "special dependency" instruction:

• Sometimes, users want to allow one parameter to depend on another. For instance, in DMC (see dmc_dm), the parameter A is positive in the congruent condition, but negative in the incongruent condition. Thus, parameters may have a 'special depencency' which can be expressed as an equation.

• To define a special dependency, users can use the operation "==". The parameter that should have the dependency is on the left-hand side, while the mathematical relationship to other parameters is defined on the right-hand side.

• This then looks something like "a ~ foo == -(a ~ bar)".

• This means that the parameter a in condition foo will always be -1 * the parameter a in condition bar. Thus, if a in condition bar has the value 5, then a in condition foo will be -5.

• The expression on the right-side can refer to any arbitrary mathematical relation.

• Important: Make sure that each 'parameter ~ condition' combination on the right-hand side of the equation are set in brackets.

• Another example: Parameter a in condition foo should be the mean of the parameter b in conditions bar and baz; this would be the instruction "a ~ foo == 0.5*(b ~ bar) + 0.5*(b ~ baz)"

The "additional/custom parameter combination" instruction:

• Sometimes, users may want to combine multiple parameters to summarize a certain property of the model. For example, in DMC (see dmc_dm), the shape and rate parameter jointly determine the peak latency.

• To avoid having to calculate this manually, users can define "custom" parameter combinations using the ":=" operation:

• An exemplary instruction might look like this: "peak_l := (a - 1) * tau"

• Expressions and values that provide calculations for those parameters are stored in a separate list cust_prms.

Value

For drift_dm objects, the updated drift_dm object.

For flex_prms, the updated flex_prms object.

See Also

flex_prms()

Examples

# Example 1: Modify a flex_prms object  directly ---------------------------
# create an auxiliary flex_prms object
a_flex_prms_obj <- flex_prms(
  c(muc = 3, b = 0.5, non_dec = 0.3),
  conds = c("foo", "bar")
)

# then carry out some "instructions". Here (arbitrary operations):
# 1.) Consider b as fixed
# 2.) Let muc vary independently for the conditions foo and bar
# 3.) Set non_dec in condition bar to be half as large as non_dec in
#     condition bar
instr <-
  "b <!>
 muc ~
 non_dec ~ bar == (non_dec ~ foo) / 2
"
modify_flex_prms(object = a_flex_prms_obj, instr = instr)


# Example 2: Modify a flex_prms object stored inside a drift_dm object -----
a_model <- ratcliff_dm() # get a model for demonstration purpose
modify_flex_prms(object = a_model, instr = "muc ~ => 4")

Constant Drift Rate

Description

Constant Drift Rate

Usage

mu_constant(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

a vector of the same length as t_vec with the drift rate for each element of the vector.

Drift Rate for DMC

Description

Provides the drift rate of the superimposed decision process. That is the derivative of the rescaled gamma function plus a constant drift rate for the controlled process.

Usage

mu_dmc(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

provides the first derivative of the superimposed process with respect to t_vec.

Integral of Constant Drift Rate

Description

Integral of Constant Drift Rate

Usage

mu_int_constant(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

a vector calculated as t_vec*muc

Integral of DMC's Drift Rate

Description

Provides the integral of the drift rate of the superimposed decision process. This is the sum of the rescaled gamma function and the linear function.

Usage

mu_int_dmc(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

provides the scaled gamma distribution function of the superimposed process for each time step in t_vec.

Drift Rate for SSP

Description

Provides the drift rate for the SSP model. That is, the sum of attention attributed to the flankers and central target, scaled by the perceptual input.

Usage

mu_ssp(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

provides the drift rate for SSP with respect to t_vec

Create a DDM model — internal

Description

This function assembles all components to create a drift_dm object.

Usage

new_drift_dm(
  flex_prms_obj,
  sigma,
  t_max,
  dt,
  dx,
  solver,
  comp_funs,
  cost_function,
  subclass,
  b_coding = NULL,
  obs_data = NULL
)

Arguments

Details

We do not perform input checks here; we just assemble the object. Any pre-wrangling is done in drift_dm(). Checks are performed by validate_drift_dm(), which is called indirectly via the setters (e.g., prms_solve() and obs_data()).

Value

A list with elements flex_prms_obj, prms_solve, solver, comp_funs, and cost_function. The object has class attributes c(subclass, "drift_dm") and an attribute "b_coding" containing the boundary coding. If obs_data is not NULL, the observed data are attached via obs_data().

See Also

drift_dm(), validate_drift_dm(), obs_data(), drift_dm_default_b_coding(), prms_solve().

Create a New stats_dm Object

Description

new_stats_dm initializes a stats_dm object to label statistic types and store necessary attributes for the custom methods (such as ⁠plot_*⁠)

Usage

new_stats_dm(stat_df, type, ...)

Arguments

Details

new_stats_dm sets up the stat_df object by assigning it the class stats_dm, along with additional classes based on the specified type. For "basic_stats", "cafs", "quantiles", "delta_funs", this will be c(type, "sum_dist", "stats_dm", "data.frame")". For "fit_stats", this will be c(type, "stats_dm", "data.frame")".

For basic stats, Conditional Accuracy Functions (CAFs), Quantiles, and Delta Functions, the function requires a b_coding argument, which specifies boundary coding details and is set as an attribute.

The function performs validation through validate_stats_dm to ensure that the stats_dm object is well formatted.

Value

An object of class stats_dm, with additional classes and attributes depending on type.

Get the Number of Observations for a drift_dm Object

Description

This method retrieves the total number of observations in the obs_data list of a drift_dm object.

Usage

## S3 method for class 'drift_dm'
nobs(object, ...)

Arguments

Details

The function iterates over each element in object$obs_data, counts the entries in each nested component, and returns the cumulative sum as the total observation count.

It was written to provide an nobs method for calculating the log-likelihood (logLik), AIC (stats::AIC), and BIC (stats::BIC) statistics for objects of type drift_dm.

Value

An integer representing the total number of observations across all conditions in object$obs_data. If obs_data doesn't exist, the function returns 0

Examples

# get a pre-built model and data set for demonstration purpose
a_model <- dmc_dm()
obs_data(a_model) <- dmc_synth_data

# then get the number of observations by accessing the model
nobs(a_model)

# same number of observations as in the original data set
nrow(dmc_synth_data)

Constant Non-Decision time

Description

A dirac delta on "non_dec", to provide a constant non-decision time.

Usage

nt_constant(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

a vector of the same length as t_vec with zeros, except for the element matching with "non_dec" with respect to "t_vec"

Truncated Normally-Distributed Non-Decision Time

Description

Truncated Normally-Distributed Non-Decision Time

Usage

nt_truncated_normal(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

returns the PDF of a truncated normal distribution for t_vec, with mean "non_dec" and standard deviation "sd_non_dec". Lower truncation is 0. Upper truncation is max(t_vec)

Uniform Non-Decision Time

Description

Uniform Non-Decision Time

Usage

nt_uniform(prms_model, prms_solve, t_vec, one_cond, ddm_opts)

Arguments

Value

returns the PDF of a uniform distribution for t_vec, centered around "non_dec" and with a range of "range_non_dec".

The Observed Data

Description

Functions to get or set the "observed data" of an object.

Usage

obs_data(object, ...) <- value

## S3 replacement method for class 'drift_dm'
obs_data(object, ..., eval_model = FALSE) <- value

obs_data(object, ...)

## S3 method for class 'drift_dm'
obs_data(object, ..., messaging = TRUE)

## S3 method for class 'fits_ids_dm'
obs_data(object, ...)

## S3 method for class 'fits_agg_dm'
obs_data(object, ...)

Arguments

Details

obs_data() is a generic accessor function, and ⁠obs_data<-()⁠ is a generic replacement function. The default methods get and set the "observed data". Their behavior, however, may be a bit unexpected.

In drift_dm objects, the observed data are not stored as a data.frame. Instead, any supplied observed data set is disassembled into RTs for the upper and lower boundary and with respect to the different conditions (ensures more speed and easier programming in the depths of the package). Yet, obs_data() returns a data.frame for drift_dm objects. This implies that obs_data() does not merely access the observed data, but re-assembles it. Consequently, a returned data.frame for the observed data is likely sorted differently than the data.frame that was originally set to the model via ⁠obs_data<-()⁠. Also, when the originally supplied data set provided more conditions than the model, the unused conditions will not be part of the returned data.frame.

For fits_ids_dm (see load_fits_ids), the observed data are stored as a data.frame in the general fit procedure info. This is the data.frame that obs_data() will return. Thus, the returned data.frame will match with the data.frame that was initially supplied to estimate_model_ids, although with unused conditions being dropped.

In theory, it is possible to update parts of the "observed data". However, because obs_data() returns a re-assembled data.frame for drift_dm objects, great care has to be taken with respect to the ordering of the argument value. A message is displayed to remind the user that the returned data.frame may be sorted differently than expected.

Value

For obs_data() a data.frame of the observed data. The method obs_data.drift_dm() per default displays a message to remind the user that the returned data.frame is likely sorted differently than expected.

For ⁠obs_data<-()⁠ the updated drift_dm object.

Note

There is only a replacement function for drift_dm objects. This is because replacing the observed data after the model has been fitted (i.e., for a fits_ids_dm object) doesn't make sense.

See Also

drift_dm()

Examples

# Set some data to a model -------------------------------------------------
my_model <- dmc_dm() # DMC is pre-built and directly available
# synthetic data suitable for DMC; comes with dRiftDM
some_data <- dmc_synth_data
obs_data(my_model) <- some_data

# Extract data from a model ------------------------------------------------
head(obs_data(my_model))

# Important: ---------------------------------------------------------------
# The returned data.frame may be sorted differently than the one initially
# supplied.
some_data <- some_data[sample(1:nrow(some_data)), ] #' # shuffle the data set
obs_data(my_model) <- some_data
all.equal(obs_data(my_model), some_data)
# so don't do obs_data(my_model)["Cond"] <- ...

# Addition: ----------------------------------------------------------------
# accessor method also available for fits_ids_dm objects
# (see estimate_model_ids)
# get an exemplary fits_ids_dm object
fits <- get_example_fits("fits_ids_dm")
head(obs_data(fits))

Disassemble an Observed Data set

Description

Takes a data.frame with columns RT, Cond, and column matching with b_coding, and disassembles it into a list of rts

Usage

obs_data_to_rt_lists(obs_data, b_coding = NULL)

Arguments

Details

performs checks on b_coding (check_b_coding) and checks/reduces obs_data (check_reduce_raw_data) before disassembling the data set.

Value

a list of rts with entries

• rts_u -> containing a list of numeric vectors, with names according to the values in Cond

• rts_l -> containing a list of numeric vectors, with names according to the values in Cond

Access the Probability Density Functions of a Model

Description

Functions to obtain the probability density functions (PDFs) of a model. These PDFs represent the convolution of the first-passage-time (decision time) with the non-decision time.

Usage

pdfs(object, ...)

## S3 method for class 'drift_dm'
pdfs(object, ...)

## S3 method for class 'fits_agg_dm'
pdfs(object, ...)

Arguments

Details

If the model has not been evaluated, re_evaluate_model() is called before returning the PDFs.

Value

A list with the entries:

• pdfs, contains another named list with entries corresponding to the conditions of the model (see conds()). Each of these elements is another named list, containing the entries pdf_u and pdf_l, which are numeric vectors for the PDFs of the upper and lower boundary, respectively.

• t_vec, containing a numeric vector of the time domain.

See Also

drift_dm(), re_evaluate_model(), conds()

Examples

# get a pre-built model for demonstration purpose
a_model <- dmc_dm()
str(pdfs(a_model))

Plot Conditional Accuracy Functions (CAFs)

Description

Visualizes conditional accuracy functions (CAFs) for observed and/or predicted data. This is useful for assessing model fit or exploring response patterns across conditions or participants.

Usage

## S3 method for class 'cafs'
plot(
  x,
  ...,
  id = NULL,
  conds = NULL,
  col = NULL,
  xlim = NULL,
  ylim = c(0, 1),
  xlab = "Bins",
  ylab = NULL,
  interval_obs = TRUE,
  interval_pred = TRUE
)

Arguments

Details

If x contains multiple IDs and no specific id is provided, the function aggregates across participants before plotting.

Observed CAFs are shown as points, and predicted CAFs as lines. When interval = TRUE and the input includes interval estimates (i.e., the column Estimate exists), the plot includes error bars for observed data and shaded contours for model predictions.

Colors, symbols, and line styles can be customized via ....

Value

Returns NULL invisibly. The function is called for its side effect of generating a plot.

Examples

# Example 1: Model predictions only ---------------------------------------
a_model <- dmc_dm()
cafs <- calc_stats(a_model, type = "cafs")
plot(cafs)
plot(cafs, col = c("green", "red"), ylim = c(0.5, 1))

# Example 2: Observed and predicted data ----------------------------------
obs_data(a_model) <- dmc_synth_data
cafs <- calc_stats(a_model, type = "cafs")
plot(cafs)

# Example 3: Observed data only -------------------------------------------
cafs <- calc_stats(dmc_synth_data, type = "cafs")
plot(cafs)

# Example 4: Observed data with interval ----------------------------------
cafs <- calc_stats(dmc_synth_data, type = "cafs", resample = TRUE)
plot(cafs)

Plot Delta Functions

Description

Visualizes delta functions for observed and/or predicted data. This is useful for assessing model fit or exploring the model behavior

Usage

## S3 method for class 'delta_funs'
plot(
  x,
  ...,
  id = NULL,
  conds = NULL,
  dv = NULL,
  col = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = "RT [s]",
  ylab = expression(Delta),
  interval_obs = TRUE,
  interval_pred = TRUE
)

Arguments

Details

If x contains multiple IDs and no specific id is provided, the function aggregates across participants before plotting.

Observed delta functions are shown as points, and predicted delta functions as lines. When interval_obs = TRUE or interval_pred = TRUE and the input includes interval estimates (i.e., the column Estimate exists), the plot includes error bars for observed data and shaded contours for model predictions.

Colors, symbols, and line styles can be customized via ....

Value

Returns NULL invisibly. The function is called for its side effect of generating a plot.

Examples

# Example 1: Model predictions only ---------------------------------------
a_model <- dmc_dm()
deltas <- calc_stats(
  a_model,
  type = "delta_funs",
  minuends = "incomp",
  subtrahends = "comp"
)
plot(deltas)
plot(deltas, col = "black", lty = 2, xlim = c(0.2, 0.65))

# Example 2: Observed and predicted data ----------------------------------
obs_data(a_model) <- dmc_synth_data
deltas <- calc_stats(
  a_model,
  type = "delta_funs",
  minuends = "incomp",
  subtrahends = "comp"
)
plot(deltas)

# Example 3: Observed data only -------------------------------------------
deltas <- calc_stats(
  dmc_synth_data,
  type = "delta_funs",
  minuends = "incomp",
  subtrahends = "comp"
)
plot(deltas)

# Example 4: Observed data with intervals ---------------------------------
deltas <- calc_stats(
  dmc_synth_data,
  type = "delta_funs",
  minuends = "incomp",
  subtrahends = "comp",
  resample = TRUE
)
plot(deltas)

Plot Distributions of Predicted and Observed Data

Description

Visualizes observed and/or predicted response time distributions. Useful for assessing model fit or exploring model behavior.

Usage

## S3 method for class 'densities'
plot(
  x,
  ...,
  id = NULL,
  conds = NULL,
  col = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = "RT [s]",
  ylab = "Density",
  obs_stats = "hist",
  interval_obs = FALSE,
  interval_pred = TRUE
)

Arguments

Details

If x contains multiple IDs and no specific id is provided, the function aggregates across participants before plotting. You can provide a vector of ids to produce separate plots for each participant.

Observed densities are shown as histograms (default: gray shaded areas), or KDE lines (default: black, dotted). Predicted densities are shown as lines (default: colorized). Distributions associated with the upper boundary are shown with values > 0 (i.e., the upper part of the plot), distributions associated with the lower boundary are shown with values < 0 (i.e., the lower part of the plot).

Axis limits, colors, and styling options can be customized via .... If interval information is provided (i.e., the column Estimate exists in x), error bars or shading will be added, depending on the type of statistic.

A legend is only displayed if there is predicted data.

Value

Returns NULL invisibly. The function is called for its side effect of generating a plot.

Examples

# Example 1: Predicted densities only -------------------------------------
a_model <- dmc_dm()
dens <- calc_stats(a_model, type = "densities")
plot(dens, xlim = c(0, 1))
plot(dens, xlim = c(0, 1), conds = "comp")


# Example 2: Observed and predicted densities -----------------------------
obs_data(a_model) <- dmc_synth_data
dens <- calc_stats(a_model, type = "densities")
plot(dens, xlim = c(0, 1), conds = "comp", col = "blue")

# Example 3: Observed densities only --------------------------------------
dens <- calc_stats(dmc_synth_data, type = "densities")
plot(dens, conds = "comp", obs.hist.col = "green", alpha = 1)

# Example 4: With interval estimates --------------------------------------
dens <- calc_stats(dmc_synth_data, type = "densities", resample = TRUE)
plot(dens, interval_obs = TRUE, conds = "comp")

Plot Components of a Drift Diffusion Model

Description

This function generates plots for all components of a drift diffusion model (DDM), such as drift rate, boundary, and starting condition. Each component is plotted against the time or evidence space, allowing for visual inspection of the model's behavior across different conditions.

Usage

## S3 method for class 'drift_dm'
plot(x, ..., conds = NULL, col = NULL, xlim = NULL, bundle_plots = TRUE)

Arguments

Details

The plot.drift_dm function provides an overview of key DDM components, which include:

• mu_fun: Drift rate over time.

• mu_int_fun: Integrated drift rate over time (if required by the specified solver of the model).

• x_fun: Starting condition as a density across evidence values.

• b_fun: Boundary values over time.

• dt_b_fun: Derivative of the boundary function over time.

• nt_fun: Non-decision time as a density over time.

Value

Nothing (NULL; invisibly)

Examples

# plot the component functions of the Ratcliff DDM
plot(ratcliff_dm())
plot(ratcliff_dm(var_non_dec = TRUE))
# Note: the variability in the drift rate for the Ratcliff DDM
# is not plotted! This is because it is not actually stored as a component
# function.

# plot the component functions of the DMC model
plot(dmc_dm(), col = c("green", "red"))

Plot MCMC Results and Diagnostics for mcmc_dm Objects

Description

Visualize MCMC results and diagnostics for mcmc_dm objects. The function plot.mcmc() is typically called when users supply an mcmc_dm object returned by estimate_dm() to the generic base::plot() function.

Usage

## S3 method for class 'mcmc_dm'
plot(x, ..., id = NULL, what = "trace", bundle_plots = TRUE)

Arguments

Details

This function provides diagnostic and summary visualizations of MCMC samples. It handles results from both hierarchical and non-hierarchical MCMC runs:

• If id is provided, the plot refers to the requested participant, with MCMC results extracted at the individual level.

• If id is omitted, plots refer to group-level parameters (i.e., the hyperparameters)

The following plot types are supported:

• Trace plots (what = "trace"): These plots show sampled parameter values across MCMC iterations for each chain. They are primarily used to inspect convergence and mixing behavior. Ideally, all chains should appear well-mixed (i.e., they should overlap and sample in a similar range). Lack of convergence is indicated by chains that remain in separate regions or exhibit trends over time.

• Density plots (what = "density"): These plots display smoothed marginal posterior distributions for each parameter, collapsed over chains and iterations. They are useful for understanding the central tendency, variance, and shape of the posterior distributions.

• Autocorrelation plots (what = "auto"): These plots display the autocorrelation at different lags, averaged across chains. They are useful to judge how quickly the chains produced independent samples.

Value

Returns NULL invisibly.

See Also

plot_mcmc_trace(), plot_mcmc_marginal(), plot_mcmc_auto()

Examples

# get an examplary `mcmc_dm` object
chains_obj <- get_example_fits("mcmc")
plot(chains_obj)
plot(chains_obj, what = "density")
plot(chains_obj, what = "density", which_prm = "b", bundle_plots = FALSE)

Plot Response Time Quantiles

Description

Visualizes response time quantiles for observed and/or predicted data across experimental conditions. This is useful for assessing model fit or exploring response patterns across conditions or participants.

Usage

## S3 method for class 'quantiles'
plot(
  x,
  ...,
  id = NULL,
  conds = NULL,
  dv = NULL,
  col = NULL,
  xlim = NULL,
  ylim = c(0, 1),
  xlab = "RT [s]",
  ylab = "F(RT)",
  interval_obs = TRUE,
  interval_pred = TRUE
)

Arguments

Details

If x contains multiple IDs and no specific id is provided, the function aggregates across participants before plotting.

Observed quantiles are shown as points, and predicted quantiles as lines. When interval = TRUE and the input includes interval estimates (i.e., the column Estimate exists), the plot includes error bars for observed data and shaded contours for model predictions.

Colors, symbols, and line styles can be customized via ....

Value

Returns NULL invisibly. The function is called for its side effect of generating a plot.

Examples

# Example 1: Model predictions only ---------------------------------------
a_model <- dmc_dm()
quantiles <- calc_stats(a_model, type = "quantiles")
plot(quantiles)
plot(quantiles, col = c("green", "red"), xlim = c(0.2, 0.6))

# Example 2: Observed and predicted data ----------------------------------
obs_data(a_model) <- dmc_synth_data
quantiles <- calc_stats(a_model, type = "quantiles")
plot(quantiles)

# Example 3: Observed data only -------------------------------------------
quantiles <- calc_stats(dmc_synth_data, type = "quantiles")
plot(quantiles)

# Example 4: Observed data with interval ----------------------------------
cafs <- calc_stats(dmc_synth_data, type = "quantiles", resample = TRUE)
plot(cafs)

Plot Multiple Statistics

Description

This function iterates over a list of statistics data, resulting from a call to calc_stats(), and subsequently plots each statistic. It allows for a simple arrangement of multiple plots on a single graphics device.

Usage

## S3 method for class 'stats_dm_list'
plot(x, ..., mfrow = NULL)

Arguments

Details

The plot.stats_dm_list() function "merely" iterates over each entry of x and calls the respective plot() method. If dRiftDM doesn't provide a plot() method for an object stored in x, the respective entry is skipped and a message is displayed.

When users want more control over each plot, it is best to call the plot() function separately for each statistic in the list (e.g., plot(x$cafs); plot(x$quantiles))

Value

Nothing (NULL; invisibly)

See Also

plot.cafs(), plot.quantiles(), plot.delta_funs(), plot.densities()

Examples

# get a list of statistics for demonstration purpose
all_fits <- get_example_fits("fits_ids_dm")
stats <- calc_stats(all_fits, type = c("cafs", "quantiles"))

# then call the plot function.
plot(stats, mfrow = c(1, 2))

Plot Traces of a Drift Diffusion Model

Description

Creates a plot of simulated traces (i.e., simulated evidence accumulation processes) from a drift diffusion model. Such plots are useful for exploring and testing model behavior.

Usage

## S3 method for class 'traces_dm_list'
plot(
  x,
  ...,
  conds = NULL,
  col = NULL,
  col_b = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = "Time",
  ylab = "Evidence"
)

## S3 method for class 'traces_dm'
plot(x, ...)

Arguments

Details

plot.traces_dm_list() iterates over all conditions and plots the traces. It includes a legend with condition labels.

plot.traces_dm plots a single set of traces. Because simulate_traces() returns an object of type traces_dm_list per default, users will likely call plot.traces_dm_list() in most cases; and not plot.traces_dm. plot.traces_dm is only relevant if users explicitly extract and provide an object of type traces_dm.

The function automatically generates the upper and lower boundaries based on the information stored within x.

Value

NULL invisibly

See Also

simulate_traces

Examples

# get a couple of traces for demonstration purpose
a_model <- dmc_dm()
some_traces <- simulate_traces(a_model, k = 3)

# Plots for traces_dm_list objects ----------------------------------------
# basic plot
plot(some_traces)

# a slightly more beautiful plot :)
plot(some_traces,
  col = c("green", "red"),
  xlim = c(0, 0.35),
  xlab = "Time [s]",
  ylab = bquote(Realizations ~ of ~ X[t]),
  legend_pos = "bottomright"
)

# Plots for traces_dm objects ---------------------------------------------
# we can also extract a single set of traces and plot them
one_set_traces <- some_traces$comp
plot(one_set_traces)

# modifications to the plot work in the same way
plot(one_set_traces,
  col = "green",
  xlim = c(0, 0.35),
  xlab = "Time [s]",
  ylab = bquote(Realizations ~ of ~ X[t]),
  legend = "just comp"
)

Plot MCMC Chains for Drift Diffusion Model Parameters

Description

The functions provide visualizations of MCMC results obtained for drift_dm() objects. Users won't call the functions directly. Instead, they are called via plot.mcmc_dm(), and the following documentation helps to define optional arguments that can be passed on via ....

Usage

plot_mcmc_trace(
  chains,
  col_palette = grDevices::rainbow,
  col_chains = NULL,
  which_prms = NULL,
  xlab = NULL,
  ylab = NULL,
  xlim = NULL,
  ylim = NULL
)

plot_mcmc_marginal(
  chains,
  col_line = NULL,
  col_shade = NULL,
  which_prms = NULL,
  xlab = NULL,
  ylab = NULL,
  xlim = NULL,
  ylim = NULL
)

plot_mcmc_auto(
  chains,
  lags = 1:30,
  col_line = NULL,
  which_prms = NULL,
  xlab = NULL,
  ylab = NULL,
  xlim = NULL,
  ylim = NULL,
  type = NULL,
  main = NULL
)

Arguments

Details

• plot_mcmc_trace() plots one panel per parameter, with lines for each chain showing how values evolve over iterations. This is useful for diagnosing convergence and mixing.

• plot_mcmc_marginal() plots smoothed marginal posterior densities collapsed over chains and iterations for each parameter, useful for inspecting posterior distributions.

• plot_mcmc_auto() plots the autocorrelation among samples as function of lags.

Value

These functions are called for their side effects (producing plots). They return NULL invisibly.

See Also

plot.mcmc_dm()

Print Functions for Stats Objects

Description

when calling calc_stats(), each returned statistic will be a subclass of stats_dm and data.frame. The following print() methods will call the more generic print.stats_dm() function.

Usage

## S3 method for class 'cafs'
print(x, ...)

## S3 method for class 'basic_stats'
print(x, ...)

## S3 method for class 'quantiles'
print(x, ...)

## S3 method for class 'delta_funs'
print(x, ...)

## S3 method for class 'fit_stats'
print(x, ...)

## S3 method for class 'sum_dist'
print(x, ...)

Arguments

Value

x (invisibly)

Summary and Printing for fits_agg_dm Objects

Description

Methods for summarizing and printing objects of the class fits_agg_dm, which contain model fits based on aggregated data across participants.

Usage

## S3 method for class 'summary.fits_agg_dm'
print(x, ..., just_header = FALSE, round_digits = drift_dm_default_rounding())

## S3 method for class 'fits_agg_dm'
summary(object, ..., select_unique = FALSE)

Arguments

Details

The summary.fits_agg_dm function creates a structured summary of a fits_agg_dm object, containing:

• summary_drift_dm_obj: A list with information about the underlying drift diffusion model (as returned by summary.drift_dm()).

• prms: Parameter estimates obtained from the model fit. This is equivalent to a call to coef.drift_dm() on the stored model object.

• obs_data: A list providing the number of individual participants and the average number of trials per condition across participants.

The print.summary.fits_agg_dm function formats and prints the above summary in a human-readable form.

Value

summary.fits_agg_dm() returns a list of class summary.fits_agg_dm (see Details for its structure).

print.summary.fits_agg_dm() returns the input object invisibly.

See Also

summary.drift_dm, coef.drift_dm

Examples

# Load example fit object
fits_agg <- get_example_fits("fits_agg")
sum_obj <- summary(fits_agg)
print(sum_obj, round_digits = 2)

Summary and Printing for fits_ids_dm Objects

Description

Methods for summarizing and printing objects of the class fits_ids_dm, which contain multiple fits across individuals.

Usage

## S3 method for class 'summary.fits_ids_dm'
print(x, ..., just_header = FALSE, round_digits = drift_dm_default_rounding())

## S3 method for class 'fits_ids_dm'
summary(object, ..., select_unique = FALSE)

Arguments

Details

The summary.fits_ids_dm function creates a summary object. The contents of this summary object depends on whether the user supplies a fits_ids_dm object that was created with estimate_dm() or the deprecated function load_fits_ids().

• In the first case, the object contains:

• summary_drift_dm_obj: A list with information about the underlying drift diffusion model (as returned by summary.drift_dm()).

• prms: All parameter values across all conditions (essentially a call to coef() with the argument select_unique = FALSE).

• stats: A named list of matrices for each condition, including mean and standard error for each parameter.

• obs_data: A list providing the number of individual participants and the average number of trials per condition across participants.

• optimizer: A string of the optimizer that was used

• conv_info: A list providing a summary of the convergance and messages for all IDs

• In the second case, the object contains:

• lower and upper: Lower and upper bounds of the search space.

• model_type: Description of the model type, based on class information.

• prms: All parameter values across all conditions (essentially a call to coef() with the argument select_unique = FALSE).

• stats: A named list of matrices for each condition, including mean and standard error for each parameter.

• N: The number of individuals.