## compute_drift_rate ##########################################################
#' Compute drift rate
#'
#' compute_drift_rate uses approaches described by:
#' - Dai & Busemeyer, J Exp Psychol Gen, 2014
#' - Scholten & Read, J Exp Psychol Learn Mem Cogn, 2013
#'
#' In general, drift rate is computed according to Dai & Busemeyer (2014):
#' v = x['w'] * (um_l - um_s) - (1 - x['w']) * (pt_l - pt_s)
#'
#' In case of the defer-speedup effect under outcome framing, we follow Scholten & Read (2013). See their p. 1197.
#'
#' @param parameters parameters
#' @param stimuli stimuli
#' @param frame the way in which the choice was framed: "defer", "speedup", "date", "delay
#' @inheritParams itchmodel::fit_model
#' @export
compute_drift_rate <- function(parameters, stimuli, parameterization = "", frame = "") {
# 1. Unlist parameters =======================================================
x <- unlist(parameters)
# 2. Compute drift rate ======================================================
# Drift rate, v, is computed as:
#
# v = w * du - (1 - w) * dp, where
#
# w is a weighting parameter
# du is the difference in utility between the large and small option
# dp is the difference in weighted/perceived time between the late and soon option
unname(x['w'] *
compute_transformation_diffs(parameters = parameters,
stimuli = stimuli,
parameterization = parameterization,
frame = frame,
variable = 'du') -
(1 - x['w']) *
compute_transformation_diffs(parameters = parameters,
stimuli = stimuli,
parameterization = parameterization,
frame = frame,
variable = 'dp'))
}
#'
#'
#' @export
compute_drift_rate_transformation_diffs <- function(parameters, stimuli, parameterization = "", frame = "") {
x <- unlist(parameters)
tibble::tibble(variable = character(),
value = double(),
weighted = logical()) %>%
tibble::add_row(variable = 'du',
value = compute_transformation_diffs(parameters = parameters,
stimuli = stimuli,
parameterization = parameterization,
frame = frame,
variable = 'du'),
weighted = FALSE) %>%
tibble::add_row(variable = 'dp',
value = compute_transformation_diffs(parameters = parameters,
stimuli = stimuli,
parameterization = parameterization,
frame = frame,
variable = 'dp'),
weighted = FALSE) %>%
tibble::add_row(variable = 'du',
value = unname(x['w']) *
compute_transformation_diffs(parameters = parameters,
stimuli = stimuli,
parameterization = parameterization,
frame = frame,
variable = 'du'),
weighted = TRUE) %>%
tibble::add_row(variable = 'dp',
value = (1 - unname(x['w'])) *
compute_transformation_diffs(parameters = parameters,
stimuli = stimuli,
parameterization = parameterization,
frame = frame,
variable = 'dp'),
weighted = TRUE) %>%
dplyr::mutate(variable = factor(variable, levels = c('du', 'dp')),
weighted = factor(weighted, levels = c('FALSE', 'TRUE'))
)
}
## compute_transformation_diffs ################################################
#' Computes the differences in utilities (du) and weighted times (dp) for the small-but-soon and large-but-later options, given a set of parameters and stimuli.
compute_transformation_diffs <- function(parameters, stimuli, parameterization = "", frame = "", variable = 'du') {
if (variable == 'du') {
# Differences in utility ===================================================
# Utility of the large option
compute_transformation(q = stimuli$m_l,
parameters = parameters,
variable = 'u_ll',
frame = frame) -
# Utility of the small option
compute_transformation(q = stimuli$m_s,
parameters = parameters,
variable = 'u_ss',
frame = frame)
} else if (variable == 'dp') {
# Differences in weighted/perceived time ===================================
# Weighted/perceived time of the late option
compute_transformation(q = stimuli$t_l,
parameters = parameters,
variable = 'p_ll',
frame = frame) -
# Weighted/perceived time of the soon option
compute_transformation(q = stimuli$t_s,
parameters = parameters,
variable = 'p_ss',
frame = frame)
} else {
NA
}
}
## compute_transformation ######################################################
#' Given a set of parameters, compute the money-to-utility and clock-time-to-perceived-time transformation functions
#'
#' @param parameters Named list of parameters, including alpha, mu, beta, and kappa
#' @param x optional double vector of values to be transformed (i.e. money/time)
#' @export
compute_transformation <- function(q, parameters, variable, frame = "") {
# 1. Unlist parameters =======================================================
x <- unlist(parameters)
# 2. Transform money/utility into utility/weighted time ======================
if (variable == 'u') {
# 2.1. Utility (used for date and delay frames) ----------------------------
unname(convert(q, type = 'power', param = c(x['alpha'], x['mu'])))
} else if (variable == 'u_ss') {
# 2.2. Utility of small-but-soon option (neutral, defer, speedup frames) ---
if (frame %in% c('neutral','speedup')) {
unname(convert(q, type = 'power', param = c(x['alpha'], 1)))
} else if (frame %in% c('defer', 'delay', 'date')) {
unname(convert(q, type = 'power', param = c(x['alpha'], x['mu'])))
} else {
NA
}
} else if (variable == 'u_ll') {
# 2.3. Utility of large-but-later option (neutral, defer, speedup frames) --
if (frame %in% c('neutral','defer')) {
unname(convert(q, type = 'power', param = c(x['alpha'], 1)))
} else if (frame %in% c('speedup', 'delay', 'date')) {
unname(convert(q, type = 'power', param = c(x['alpha'], x['mu'])))
} else {
NA
}
} else if (variable %in% c('p', 'p_ss', 'p_ll')) {
# 2.3. Weighted time (any frame) -------------------------------------------
unname(convert(q, type = 'power', param = c(x['beta'], x['kappa'])))
} else {
NA
}
}
## convert #####################################################################
#' Convert a quantity into another according to a transformation function. Transformation functions are 'power' (based on Dai & Busemeyer, J Exp Gen Psychol, 2014) and 'log' (based on Scholten & Read, J Exp Psychol Learn Mem Cogn, 2013).
#'
#' @param q quantity (either money in Euros or time in days)
#' @param type transformation function
#' @param params transformation function parameters
#' @export
convert <- function(q, type = 'power', param) {
# 1. Assertions ==============================================================
assertthat::assert_that(all(q >= 0), msg = 'Make sure that money and time should be ≥ 0')
assertthat::assert_that(is.character(type))
assertthat::assert_that(is.vector(param))
assertthat::assert_that(length(param) == 2)
assertthat::assert_that(is.numeric(param))
# 2. Transform data ==========================================================
switch(type,
power = param[2] * q^param[1],
log = param[2] * (1 / param[1] * log(1 + param[1] * q))
)
}
## db_bin_choice_prob ##########################################################
#' Dai & Busemeyer equation for binary choice probability of a specific option
#'
#' Note that this gives same results as rtdists::rdiffusion with identical parameters (but note that the latter has bound separation parameter a, which equals 2 * threshold in Dai & Busemeyer equation)
#'
#' The following should get the same output:
#' a <- runif(n = 1, min = 0.5, max = 2) # threshold separation
#' v <- runif(n = 1, min = -5, max = 5) # drift rate
#' s <- 1 # diffusion constant
#' t0 <- runif(n = 1, min = 0, max = 2) # nondecision time
#' itchmodel::db_bin_choice_prob(d = v, s = s, a = a, z = 0)
#'
#' rtdists::pdiffusion(rt = Inf, response = "upper", a = a, v = v, t0 = t0, s = s, z = 0.5*a)
#'
#' @export
db_bin_choice_prob <- function(d, s, a, z) {
theta <- a / 2
if (any(dplyr::near(d,0, tol = 1e-9))) {
d <- d + 1e-9
}
# Dividing infinite numbers results in NaNs; make negative drift rates smaller to prevent this
d <- ifelse((-2 * d * (theta + z) / s^2 > log(.Machine$double.xmax) |
-4 * d * theta / s^2 > log(.Machine$double.xmax)),
max(c((log(.Machine$double.xmax) * s^2) / (-2 * (theta + z)),
(log(.Machine$double.xmax) * s^2) / (-4 * theta))),
d)
# if (-2 * d * (theta + z) / s^2 > log(.Machine$double.xmax) |
# -4 * d * theta / s^2 > log(.Machine$double.xmax)) {
# d <-
# max(c((log(.Machine$double.xmax) * s^2) / (-2 * (theta + z)),
# (log(.Machine$double.xmax) * s^2) / (-4 * theta)))
# }
numerator <- (1 - exp(-2 * d * (theta + z) / s^2))
denominator <- (1 - exp(-4 * d * theta / s^2))
numerator / denominator
}
## dft_cp ######################################################################
#' Compute choice probability density, according to decision field theory
#'
#' Code is based on the function CP.m by Junyi Dai
#'
#' @export
#'
dft_cp <- function(d, s, theta, z, response) {
# Replace theta's with unlikely low values
theta[theta < 1e-4] <- 1e-4
# Invert d and z when response is lower
d[response == "lower"] <- -d[response == "lower"]
zvec = numeric(length(response))
zvec[response == "lower"] <- -z
zvec[response == "upper"] <- z
# Make vector for storing probabilities
cp = numeric(length(response))
# Set probability density of unlikely responses to 0.001
unlikely <- -4 * d * theta / s^2 > 500
cp[unlikely] <- 0.001
d[d == 0] <- 1e-100
# Compute probability densities for likely responses
cp[!unlikely] <-
expm1(-2 * d[!unlikely] * (theta[!unlikely] + zvec[!unlikely]) / s[!unlikely]^2) /
expm1(-4 * d[!unlikely] * theta[!unlikely] / s[!unlikely]^2)
# Make sure that cp cannot be 1 or 0
cp[cp > 0.9999] = 0.9999
cp[cp < 0.0001] = 0.0001
# if (-4 * d * theta / s^2 > 500) {
# cp <- .0001
# } else {
#
# if (d == 0) {
# d = 1e-100
# }
#
# cp <- expm1(-2 * d * (theta + z) / s^2) / expm1(-4 * d * theta / s^2)
# cp <- max(min(cp, .9999), .0001)
# }
return(cp)
}
## dft_dpd #####################################################################
#' Compute defective probability density of a specific response time , according
#' to to decision field theory
#'
#' Code is based on the function DPD.m by Junyi Dai
#'
#' @export
#'
dft_dpd <- function(d, s, theta, z, response, rtd) {
if (response == "lower") {
d <- -d
z <- -z
}
tolerance <- 1e-200
j <- 1
dpd <- 0
step_size <-
pi * (2 * theta / s)^(-2) * j *
exp(d * (theta - z) / s^2 -
((pi * j * s / (2 * theta))^2 + (d / s)^2) * rtd / 2)
while (step_size > tolerance) {
dpd <- dpd + step_size * sin(pi * (theta - z) * j / (2 * theta))
j <- j + 1
step_size <-
pi * (2 * theta / s)^(-2) * j *
exp(d * (theta - z) / s^2 -
((pi * j * s / (2 * theta))^2 + (d / s)^2) * rtd / 2);
}
if (dpd > 0) {
pd <- dpd
} else {
pd <- 1e-200
}
return(pd)
}
## eq_list_elements ############################################################
#' Ensure that list elemens are equal across conditions
#'
#' Parameter values not varying between conditions are repeated
#'
eq_list_elements <- function(lst, n_cond) {
expand_list <- function(x, n) {
if (length(x) == 1) {
x = rep(x, length.out = n)
} else if (length(x) == n) {
x = x
}
}
if (base::missing(n_cond)) {
n_cond <- max(unlist(purrr::map(lst, length)))
}
purrr::map(lst, expand_list, n_cond)
}
## fit_model ###################################################################
#' Function for fitting intertemporal choice diffusion model
#'
#' @param data Tibble with nested data
#' @param model Name of the model. Currently, three models are implemented: decision field theory - choices only (DFT_C), decision field theory - choices and response times (DFT_CRT), and the drift diffusion model (DDM)
#' @param parameterization Name of the parameterization. Currently, there are nine possibilities: (1) date_delay_time_scaling; (2) date_delay_time_scaling_t0; (3) date_delay_value_scaling; (4) date_delay_value_scaling_t0; (5) defer_speedup_value_scaling; (6) defer_speedup_value_scaling_t0; (7) defer_speedup_time_scaling; (8) defer_speedup_time_scaling_t0, and (9) one_condition. The parameterizations in which t0 varies between conditions is only available for DFT_CRT and DDM.
#' @param control A list of control parameters for the Differential Evolution algorithm
fit_model <- function(data,
model = "",
parameterization = "",
lowers = get_par_bounds(model = model,
parameterization = parameterization,
bound = 'lower'),
uppers = get_par_bounds(model = model,
parameterization = parameterization,
bound = 'upper'),
control = list(itermax = 250, reltol = 1e-5)
) {
# 1. Assert that inputs meet certain conditions ==============================
assertthat::assert_that(model %in% c("DDM", "DFT_C", "DFT_CRT"),
msg = "model is ill-specified. Currently, only the following models are supported: \n'DDM', 'DFT_C', and 'DFT_CRT'")
assertthat::assert_that(parameterization %in% c('one_condition',
'date_delay_time_scaling',
'date_delay_time_scaling_t0',
'date_delay_time_and_value_scaling',
'date_delay_time_and_value_scaling_t0',
'date_delay_value_scaling',
'date_delay_value_scaling_t0',
'defer_speedup_time_scaling',
'defer_speedup_time_scaling_t0',
'defer_speedup_time_and_value_scaling',
'defer_speedup_time_and_value_scaling_t0',
'defer_speedup_value_scaling',
'defer_speedup_value_scaling_t0'),
msg = "parameterization is ill-specified. It should be any of the following: \n'date_delay_time_scaling', \n'date_delay_time_scaling_t0', \n'date_delay_time_and_value_scaling', \n'date_delay_time_and_value_scaling_t0', \n'date_delay_value_scaling', \n'date_delay_value_scaling_t0', \n'defer_speedup_time_scaling', \n'defer_speedup_time_scaling_t0', \n'defer_speedup_time_and_value_scaling', \n'defer_speedup_time_and_value_scaling_t0', \n'defer_speedup_value_scaling', \n'defer_speedup_value_scaling_t0', or \n'one_condition', .")
assertthat::assert_that(is.double(lowers), msg = "lowers should be a double-precision vector")
assertthat::assert_that(is.double(uppers), msg = "uppers should be a double-precision vector")
# 2. Use Differential Evolution to optimize model parameters =================
optim_out <-
DEoptim::DEoptim(fn = itchmodel::get_log_likelihood, # optimization function
lower = lowers,
upper = uppers,
data = data,
model = model,
parameterization = parameterization,
control = control
)
# Nonlinear Constrained and Unconstrained Optimization via Differential Evolution
# optim_out <-
# DEoptimR::JDEoptim(fn = itchmodel::get_log_likelihood, # optimization function
# lower = lowers,
# upper = uppers,
# constr = get_nonlinear_constraints(data = data,
# model = model,
# parameterization = parameterization),
# data = data,
# model = model,
# parameterization = parameterization,
# control = control
# )
# Add model and parameterization to output
attributes(optim_out)$model = model
attributes(optim_out)$parameterisation = parameterization
# 3. Return output ===========================================================
optim_out
}
## get_default_parameter_values ################################################
#' Returns default values for model parameters
#'
#' @export
get_default_parameter_values <- function() {
list('alpha' = 0.85, # Sensitivity of the value function
'mu' = 1, # Scaling of the value function
'beta' = 0.70, # Sensitivity of the time function
'kappa' = 1, # Scaling of the time function
'w' = 0.5, # Attentional weight
"a" = 1, # Bound separation (threshold equals a * 2)
"t0" = 0.3 # Non-decision time,
)
}
## get_fit_stats ###############################################################
#' Returns model fit statistics
#'
#' @inheritParams fit_model
#' @param optim_output output from DEoptim::DEoptim
#' @param n_data_points number of data points (trials)
#' DEoptimR minimizes negative log likelihood. get_fit_stats computes log likelihood (LL), Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) based on that.
#'
#' @export
get_fit_stats <- function(optim_output, algorithm = "DEoptimR", model = "", parameterization = "", n_data_points, max_iter=0) {
# Number of free parameters
n_free_param <- get_n_free_param(model = model,
parameterization = parameterization)
if (algorithm == "DEoptimR") {
# Log-likelihood
LL <- -optim_output$value
n_iter <- optim_output$iter
converged <- ifelse(optim_output$convergence == 0, "TRUE", "FALSE")
} else if (algorithm == "DEoptim") {
LL <- -optim_output$optim$bestval
n_iter <- optim_output$optim$iter
converged <- ifelse(n_iter < max_iter,
TRUE,
FALSE)
}
# Return fit statistics
tibble::tibble(model = model,
parameterization = parameterization,
n_iter = n_iter,
converged = converged,
LL = LL,
AIC = -2 * LL + 2 * n_free_param,
BIC = -2 * LL + log(n_data_points) * n_free_param
)
}
## get_frames ##################################################################
#' Returns the names of the framing conditions, given a parameterization
#'
#' @export
get_frames <- function(parameterization) {
# 1. Assert that inputs meet expectations ====================================
assertthat::assert_that(parameterization %in% c("one_condition",
"date_delay_time_scaling",
"date_delay_time_scaling_t0",
"date_delay_time_and_value_scaling",
"date_delay_time_and_value_scaling_t0",
"date_delay_value_scaling",
"date_delay_value_scaling_t0",
"defer_speedup_time_scaling",
"defer_speedup_time_scaling_t0",
"defer_speedup_time_and_value_scaling",
"defer_speedup_time_and_value_scaling_t0",
"defer_speedup_value_scaling",
"defer_speedup_value_scaling_t0"
),
msg = "Specify parameterization as \n\"one-condition\" (i.e., estimate parameters from a single experimental condition), \n\"date_delay_value_scaling\" (i.e., estimate free value function scaling parameter across two experimental conditions), \n\"date_delay_value_scaling_t0\" (i.e., estimate free value function scaling parameter and t0 across two experimental conditions), \n\"date_delay_time_scaling\" (i.e., estimate free time function sensitivity parameter across two experimental conditions), \n\"date_delay_time_scaling_t0\" (i.e., estimate free time function sensitivity parameter and t0 across two experimental conditions), \n\"date_delay_time_and_value_scaling\", \n\"date_delay_time_and_value_scaling_t0\", \n\"defer_speedup_value_scaling\" (i.e., estimate free value function scaling parameter across two experimental conditions), \n\"defer_speedup_value_scaling_t0\" (i.e., estimate free value function scaling parameter and t0 across two experimental conditions), \n\"defer_speedup_time_scaling\" (i.e., estimate free time function scaling parameter across two experimental conditions), \n\"defer_speedup_time_scaling_t0\" (i.e., estimate free time function scaling parameter and t0 across two experimental conditions), \n\"defer_speedup_time_and_value_scaling\", or \n\"defer_speedup_time_and_value_scaling_t0\"."
)
# 2. Return frames ===========================================================
if (startsWith(parameterization, "one_condition")) {
'neutral'
} else if (startsWith(parameterization, "date_delay")) {
c('delay', 'date')
} else if (startsWith(parameterization, "defer_speedup")) {
c('neutral', 'defer', 'speedup')
} else {
NA
}
}
## get_log_likelihood #########################################################
#' Get the log-likelihood of the parameters, given the data
#'
#' This function is a wrapper around ll_diffusion.
#'
#' @param x parameters
#' @param pcrit critical probabilities, only used by get_nonlinear_constraints
#' @inheritParams fit_model
#' @export
get_log_likelihood = function(x, data, model = "DFT_C", parameterization = "", pcrit = c(.99, .01)) {
# 1. Get parameter values ==================================================
params <- get_par_values(x, model = model, parameterization = parameterization)
# 2. Return negative log-likelihood for function minimization ================
ll <- 0
if (model == "DFT_C") {
if (parameterization == "one_condition") {
ll <-
ll +
ll_dft(x = params,
stimuli = data$stimuli[[1]],
frame = data$frame[[1]],
observations = data$observations[[1]],
rt = FALSE)
} else {
for (i_cond in 1:length(params)) {
ll <-
ll +
ll_dft(x = params[[i_cond]],
stimuli = data$stimuli[[i_cond]],
frame = data$frame[[i_cond]],
observations = data$observations[[i_cond]],
rt = FALSE)
}
}
} else if (model == "DFT_CRT") {
if (parameterization == "one_condition") {
ll <-
ll +
ll_dft(x = params,
stimuli = data$stimuli[[1]],
frame = data$frame[[1]],
observations = data$observations[[1]],
rt = TRUE)
} else {
for (i_cond in 1:length(params)) {
ll <-
ll +
ll_dft(x = params[[i_cond]],
stimuli = data$stimuli[[i_cond]],
frame = data$frame[[i_cond]],
observations = data$observations[[i_cond]],
rt = TRUE)
}
}
} else if (model == "DDM") {
if (parameterization == "one_condition") {
ll <-
ll +
ll_diffusion(x = params,
stimuli = data$stimuli[[1]],
frame = data$frame[[1]],
observations = data$observations[[1]])
} else {
for (i_cond in 1:length(params)) {
ll <-
ll +
ll_diffusion(x = params[[i_cond]],
stimuli = data$stimuli[[i_cond]],
frame = data$frame[[i_cond]],
observations = data$observations[[i_cond]])
}
}
}
ll
}
## get_par_names ###############################################################
#' Get parameter names, given a model and parameterization
#'
#' @inheritParams fit_model
#' @export
get_par_names = function(model = "DFT_C", parameterization = "") {
# 1. Get model parameter names ===============================================
list(DDM =
list(
# 1.1. One condition ------------------------------------------------
one_condition = c("alpha", "mu", "beta", "kappa", "w", "a", "t0"),
# 1.2.1.1. Date/delay effect - changes in time scaling (kappa) ----------
# - delay framing: kappa = 1
# - date framing: kappa = kappa_loss < 1
date_delay_time_scaling = c("alpha", "mu", "beta", "kappa1", "kappa2", "w", "a", "t0"),
# 1.2.1.2. Date/delay effect - changes in time scaling (kappa) & t0----
# - delay framing: kappa = 1, t0 = t01
# - date framing: kappa = kappa_loss < 1, , t0 = t02
date_delay_time_scaling_t0 = c("alpha", "mu", "beta", "kappa1", "kappa2", "w", "a", "t01", "t02"),
# 1.2.2.1. Date/delay effect - changes in value scaling (mu) ----------
# - delay framing: mu = 1
# - date framing: mu = mu_gain > 1
date_delay_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "a", "t0"),
# 1.2.2.2. Date/delay effect - changes in value scaling (mu) & t0 -----
# - delay framing: mu = 1, t0 = t01
# - date framing: mu = mu_gain > 1, t0 = t02
date_delay_value_scaling_t0 = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "a", "t01", "t02"),
# 1.2.3.1. Date/delay effect - changes in time scaling (kappa) & value scaling (mu) ----------
# - delay framing: kappa = 1, mu = 1
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1
date_delay_time_and_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "w", "a", "t0"),
# 1.2.3.2. Date/delay effect - changes in time scaling (kappa), value scaling (mu), & t0 ----
# - delay framing: kappa = 1, mu = 1, t0 = t01
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1, t0 = t02
date_delay_time_and_value_scaling_t0 = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "w", "a", "t01", "t02"),
# 1.3.1.1. Defer/speedup effect - changes in time scaling (kappa) -----
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals)
# - speedup framing: kappa < 1 (under-responsive to speedups)
defer_speedup_time_scaling = c("alpha", "mu", "beta", "kappa1", "kappa2", "kappa3", "w", "a", "t0"),
# 1.3.1.2. Defer/speedup effect - changes in time scaling (kappa) & t0
# - neutral framing: kappa = 1, t0 = t01
# - defer framing: kappa > 1 (over-responsive to deferrals), t0 = t02
# - speedup framing: kappa < 1 (under-responsive to speedups), t0 = t03
defer_speedup_time_scaling_t0 = c("alpha", "mu", "beta", "kappa1", "kappa2", "kappa3", "w", "a", "t01", "t02", "t03"),
# 1.3.2.1. Defer/speedup effect - changes in value scaling (mu) -------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "a", "t0"),
# 1.3.2.2. Defer/speedup effect - changes in value scaling (mu) & t0 --
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling_t0 = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "a", "t01", "t02", "t03"),
# 1.3.3.1. Defer/speedup effect - changes in time scaling (kappa) & value scaling (mu) -----
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du)
# - speedup framing: kappa < 1 (under-responsive to speedups), mu > 1 (larger du)
defer_speedup_time_and_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "kappa3", "w", "a", "t0"),
# 1.3.3.2. Defer/speedup effect - changes in time scaling (kappa), value scaling (mu), & t0
# - neutral framing: kappa = 1, t0 = t01
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du), t0 = t02
# - speedup framing: kappa < 1 (under-responsive to speedups), mu > 1 (larger du), t0 = t03
defer_speedup_time_and_value_scaling_t0 = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "kappa3", "w", "a", "t01", "t02", "t03")
),
DFT_C =
list(
# 1.1. One condition ------------------------------------------------
one_condition = c("alpha", "mu", "beta", "kappa", "w", "theta_star"),
# 1.2.1. Date/delay effect - changes in time scaling (kappa) ----------
# - delay framing: kappa = 1
# - date framing: kappa = kappa_loss < 1
date_delay_time_scaling = c("alpha", "mu", "beta", "kappa1", "kappa2", "w", "theta_star"),
# 1.2.2. Date/delay effect - changes in value scaling (mu) ----------
# - delay framing: mu = 1
# - date framing: mu = mu_gain > 1
date_delay_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "theta_star"),
# 1.2.3. Date/delay effect - changes in time scaling (kappa) and value scaling ----------
# - delay framing: kappa = 1, mu = 1
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1
date_delay_time_scaling = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "w", "theta_star"),
# 1.3.1. Defer/speedup effect - changes in time scaling (kappa) -----
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals)
# - speedup framing: kappa < 1 (under-responsive to speedups)
defer_speedup_time_scaling = c("alpha", "mu", "beta", "kappa1", "kappa2", "kappa3", "w", "theta_star"),
# 1.3.2. Defer/speedup effect - changes in value scaling (mu) -------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "theta_star"),
# 1.3.3. Defer/speedup effect - changes in time scaling (kappa) and value scaling (mu) -----
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du)
# - speedup framing: kappa < 1 (under-responsive to speedups), , mu > 1 (larger du)
defer_speedup_time_and_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "kappa3", "w", "theta_star")
),
DFT_CRT =
list(
# 1.1. One condition ------------------------------------------------
one_condition = c("alpha", "mu", "beta", "kappa", "w", "theta_star", "t0"),
# 1.2.1.1. Date/delay effect - changes in time scaling (kappa) ----------
# - delay framing: kappa = 1
# - date framing: kappa = kappa_loss < 1
date_delay_time_scaling = c("alpha", "mu", "beta", "kappa1", "kappa2", "w", "theta_star", "t0"),
# 1.2.1.2. Date/delay effect - changes in time scaling (kappa) & t0----
# - delay framing: kappa = 1, t0 = t01
# - date framing: kappa = kappa_loss < 1, , t0 = t02
date_delay_time_scaling_t0 = c("alpha", "mu", "beta", "kappa1", "kappa2", "w", "theta_star", "t01", "t02"),
# 1.2.2.1. Date/delay effect - changes in value scaling (mu) ----------
# - delay framing: mu = 1
# - date framing: mu = mu_gain > 1
date_delay_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "theta_star", "t0"),
# 1.2.2.2. Date/delay effect - changes in value scaling (mu) & t0 -----
# - delay framing: mu = 1, t0 = t01
# - date framing: mu = mu_gain > 1, t0 = t02
date_delay_value_scaling_t0 = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "theta_star", "t01", "t02"),
# 1.2.3.1. Date/delay effect - changes in time scaling (kappa) & value scaling (mu) ----------
# - delay framing: kappa = 1, mu = 1
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1
date_delay_time_and_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "w", "theta_star", "t0"),
# 1.2.3.2. Date/delay effect - changes in time scaling (kappa), value scaling (mu), & t0 ----
# - delay framing: kappa = 1, mu = 1, t0 = t01
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1, t0 = t02
date_delay_time_and_value_scaling_t0 = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "w", "theta_star", "t01", "t02"),
# 1.3.1.1. Defer/speedup effect - changes in time scaling (kappa) -----
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals)
# - speedup framing: kappa < 1 (under-responsive to speedups)
defer_speedup_time_scaling = c("alpha", "mu", "beta", "kappa1", "kappa2", "kappa3", "w", "theta_star", "t0"),
# 1.3.1.2. Defer/speedup effect - changes in time scaling (kappa) & t0
# - neutral framing: kappa = 1, t0 = t01
# - defer framing: kappa > 1 (over-responsive to deferrals), t0 = t02
# - speedup framing: kappa < 1 (under-responsive to speedups), t0 = t03
defer_speedup_time_scaling_t0 = c("alpha", "mu", "beta", "kappa1", "kappa2", "kappa3", "w", "theta_star", "t01", "t02", "t03"),
# 1.3.2.1. Defer/speedup effect - changes in value scaling (mu) -------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "theta_star", "t0"),
# 1.3.2.2. Defer/speedup effect - changes in value scaling (mu) & t0 --
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling_t0 = c("alpha", "mu1", "mu2", "beta", "kappa", "w", "theta_star", "t01", "t02", "t03"),
# 1.3.3.1. Defer/speedup effect - changes in time scaling (kappa) & value scaling (mu) -----
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du)
# - speedup framing: kappa < 1 (under-responsive to speedups), mu > 1 (larger du)
defer_speedup_time_and_value_scaling = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "kappa3", "w", "theta_star", "t0"),
# 1.3.3.2. Defer/speedup effect - changes in time scaling (kappa), value scaling (mu), & t0
# - neutral framing: kappa = 1, t0 = t01
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du), t0 = t02
# - speedup framing: kappa < 1 (under-responsive to speedups), mu > 1 (larger du), t0 = t03
defer_speedup_time_and_value_scaling_t0 = c("alpha", "mu1", "mu2", "beta", "kappa1", "kappa2", "kappa3", "w", "theta_star", "t01", "t02", "t03")
)
)[[model]][[parameterization]]
}
## get_n_free_param ##############################################################
#' Get number of free parameters, given parameterization
#'
#' @inheritParams get_par_bounds
#' @export
get_n_free_param = function(model = "DFT_C", parameterization = "") {
# 1. Get lower and upper bounds for all parameters ===========================
l <- get_par_bounds(model = model,
parameterization = parameterization,
bound = 'lower')
u <- get_par_bounds(model = model,
parameterization = parameterization,
bound = 'upper')
# 2. Determine number of free parameters based on bounds =====================
sum(!(l == u))
}
## get_nonlinear_constraints ###################################################
#' Get nonlinear constraints
#' @param x vector of parameters
#' @param model
#' @param parameterizationn model parameterization
#' @export
get_nonlinear_constraints <- function(x, data, model = "DFT_C", parameterization = "", pcrit = c(.99, .01)) {
# Formulate subfunction for computing probabilities of choosing LL option
comp_p_ll <- function(resp, rt, parameters, model, du, dp, d, errorval) {
if (model == "DDM") {
p_ll <-
tryCatch(rtdists::pdiffusion(rt = rep(Inf, length(d)), # Note this should be Inf rather than actual RT for this purpose, because we're interested in cumulative p_ll across all possible RTs
response = rep("upper", length(d)),
a = parameters["a"],
v = d,
t0 = parameters["t0"],
s = 1,
z = 0.5 * parameters["a"]),
error = function(e) errorval)
} else if (model %in% c("DFT_C", "DFT_CRT")) {
# Probability of sampling from or attending to a specific attribute at any
# time is assumed to equal the corresponding attention weight. See Eq. 15 in
# Dai & Busemeyer
s = unname(sqrt(parameters["w"] * du^2 + (1 - parameters["w"]) * dp^2 - d^2))
p_ll <-
tryCatch(dft_cp(d = d,
s = s,
theta = unname(parameters["theta_star"] * s),
z = 0,
response = rep("upper", length(d))
),
error = function(e) errorval)
}
# Replace NaNs, if any; otherwise DEoptimR crashes
p_ll[is.na(p_ll)] <- errorval
return(p_ll)
}
# Get parameter values =======================================================
x_named <- get_par_values(x, model = model, parameterization = parameterization)
# Get frames =================================================================
frames <- as.character(data$frame)
# Formulate linear inequalities ==============================================
lineq <- double()
for (i_frame in seq(length(frames))) {
frame <- frames[i_frame]
parameters <- unlist(x_named[i_frame])
stimuli <- data$stimuli[data$frame == frame][[1]]
# Data need to be sorted for rtdists::pdiffusion
obs <-
data$observations[[i_frame]] %>%
dplyr::arrange(rt)
resp <- obs$response
rt <- obs$rt
# Constraint 1: When SS amount is 0, P(choose LL) >= 0.99 ------------------
du_1 <-
compute_transformation(q = stimuli$m_l,
parameters = parameters,
variable = 'u_ll',
frame = frame) -
compute_transformation(q = rep(0,length(stimuli$m_l)),
parameters = parameters,
variable = 'u_ss',
frame = frame)
dp_1 <-
compute_transformation(q = stimuli$t_l,
parameters = parameters,
variable = 'p_ll',
frame = frame) -
compute_transformation(q = stimuli$t_s,
parameters = parameters,
variable = 'p_ss',
frame = frame)
d_1 <- unname(parameters["w"] * du_1 - (1 - parameters["w"]) * dp_1)
p_ll_1 <-
comp_p_ll(resp, rt, parameters, model, du = du_1, dp = dp_1, d = d_1, errorval = 0)
# P(LL) >= p_crit_1 for all trials
lineq <- c(lineq, pcrit[1] - p_ll_1)
# Constraint 2: When SS and LL amounts are equal, P(choose LL) <= 0.01 -----
du_2 <-
compute_transformation(q = stimuli$m_l,
parameters = parameters,
variable = 'u_ll',
frame = frame) -
compute_transformation(q = stimuli$m_l,
parameters = parameters,
variable = 'u_ss',
frame = frame)
dp_2 <-
compute_transformation(q = stimuli$t_l,
parameters = parameters,
variable = 'p_ll',
frame = frame) -
compute_transformation(q = stimuli$t_s,
parameters = parameters,
variable = 'p_ss',
frame = frame)
d_2 <- unname(parameters["w"] * du_2 - (1 - parameters["w"]) * dp_2)
p_ll_2 <-
comp_p_ll(resp, rt, parameters, model, du = du_2, dp = dp_2, d = d_2, errorval = 1)
# P(LL) >= p_crit_1 for all trials
lineq <- c(lineq, p_ll_2 - pcrit[2])
}
# Output
unname(lineq)
}
## get_par_bounds ##############################################################
#' Get parameter lower and upper bounds, given a model and parameterization
#'
#' These minimum and maximum values are based on reported values in the literature or best guesses, see code for details.
#'
#' @inheritParams fit_model
#' @param bound bound for which to obtain values, either "lower" or "upper"
#' @export
get_par_bounds = function(model = "DFT_C", parameterization = "", bound = "lower", bound_setting = "standard") {
# 1. Define lower and upper bounds for all parameters ========================
# Values are based on the following sources:
# DB_JEPG_2014 - Dai & Busemeyer, J Exp Psychol Gen, 2014 - p. 1512
# SR_JEPLMC_2013 - Scholten & Read, J Exp Psychol Learn Mem Cogn, 2013, p. 1197
# rtdists - documentation accompanying the R package rtdists
if (bound_setting == "standard") {
lowers = list('alpha' = .01, # DB_JEPG_2014
'mu' = 1, # SR_JEPLMC_2013
'mu_gain' = 1,
'mu_loss' = 1, # SR_JEPLMC_2013
'beta' = .01, # DB_JEPG_2014
'kappa' = 1, # SR_JEPLMC_2013
'kappa_gain' = 0,
'kappa_loss' = 1,
'w' = 0.001, # DB_JEPG_2014 & SR_JEPLMC_2013; note that w is primarily a scaling parameter
"a" = 0.1, # Adjusted by BBZ; rtdists: 0.5
"t0" = 0.05, # rtdists; N.B. lowers identical across different parameterizations
"theta_star" = 0.01 # DB_JEPG_2014
)
uppers = list('alpha' = 2, # DB_JEPG_2014
'mu' = 1, #
'mu_gain' = 3,
'mu_loss' = 3, # SR_JEPLMC_2013
'beta' = 2, # DB_JEPG_2014
'kappa' = 1, # SR_JEPLMC_2013
'kappa_gain' = 1,
'kappa_loss' = 10, # guess
'w' = 1 - 0.001, # DB_JEPG_2014 & SR_JEPLMC_2013; note that w is primarily a scaling parameter
"a" = 10, # Adjusted by BBZ; rtdists: 2
"t0" = 3, # DB_JEPG_2014 (data in Table 10); N.B. uppers identical across different parameterizations
"theta_star" = 100 # DB_JEPG_2014
)
} else if (bound_setting == "wide") {
# Wider bounds for alpha, beta, w
lowers = list('alpha' = .01, # DB_JEPG_2014
'mu' = 1, # SR_JEPLMC_2013
'mu_gain' = 0, # This allows fitting when people choose LL more often on defer than speedup frames
'mu_loss' = 0, # This allows fitting when people choose LL more often on defer than speedup frames
'beta' = .01, # DB_JEPG_2014
'kappa' = 1, # SR_JEPLMC_2013
'kappa_gain' = 0,
'kappa_loss' = 0,
'w' = 0.00001, # DB_JEPG_2014 & SR_JEPLMC_2013; note that w is primarily a scaling parameter
"a" = 0.001, # Adjusted by BBZ; rtdists: 0.5
"t0" = 0.05, # rtdists; N.B. lowers identical across different parameterizations
"theta_star" = 0.01 # DB_JEPG_2014
)
uppers = list('alpha' = 3, # DB_JEPG_2014
'mu' = 1, #
'mu_gain' = 3,
'mu_loss' = 3, # SR_JEPLMC_2013
'beta' = 3, # DB_JEPG_2014
'kappa' = 1, # SR_JEPLMC_2013
'kappa_gain' = 10,
'kappa_loss' = 10, # guess
'w' = 1 - 0.00001, # DB_JEPG_2014 & SR_JEPLMC_2013; note that w is primarily a scaling parameter
"a" = 10, # Adjusted by BBZ; rtdists: 2
"t0" = 3, # DB_JEPG_2014 (data in Table 10); N.B. uppers identical across different parameterizations
"theta_star" = 100 # DB_JEPG_2014
)
} else if (bound_setting == "wide_a") {
# Wider bounds for alpha, beta, w
lowers = list('alpha' = .01, # DB_JEPG_2014
'mu' = 1, # SR_JEPLMC_2013
'mu_gain' = 0, # This allows fitting when people choose LL more often on defer than speedup frames
'mu_loss' = 0, # This allows fitting when people choose LL more often on defer than speedup frames
'beta' = .01, # DB_JEPG_2014
'kappa' = 1, # SR_JEPLMC_2013
'kappa_gain' = 0,
'kappa_loss' = 0,
'w' = 0.00001, # DB_JEPG_2014 & SR_JEPLMC_2013; note that w is primarily a scaling parameter
"a" = 0.001, # Adjusted by BBZ; rtdists: 0.5
"t0" = 0.05, # rtdists; N.B. lowers identical across different parameterizations
"theta_star" = 0.01 # DB_JEPG_2014
)
uppers = list('alpha' = 3, # DB_JEPG_2014
'mu' = 1, #
'mu_gain' = 3,
'mu_loss' = 3, # SR_JEPLMC_2013
'beta' = 3, # DB_JEPG_2014
'kappa' = 1, # SR_JEPLMC_2013
'kappa_gain' = 10,
'kappa_loss' = 10, # guess
'w' = 1 - 0.00001, # DB_JEPG_2014 & SR_JEPLMC_2013; note that w is primarily a scaling parameter
"a" = 100, # Adjusted by BBZ; rtdists: 2
"t0" = 3, # DB_JEPG_2014 (data in Table 10); N.B. uppers identical across different parameterizations
"theta_star" = 100 # DB_JEPG_2014
)
}
# Put in a vector
l <- unlist(lowers)
u <- unlist(uppers)
# 2. Return correct values, based on model, parameterization, and bound ======
unname(
list(DDM =
list(one_condition =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'], l['w'], l['a'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'], u['w'], u['a'], u['t0'])
),
# 1. Date-delay effect =========================================
# 1.1.1. Changes in time scaling (kappa) -------------------------
# - delay framing: kappa = 1
# - date framing: kappa = kappa_loss < 1
date_delay_time_scaling =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'], l['kappa_gain'], l['w'], l['a'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'], u['kappa_gain'], u['w'], u['a'], u['t0'])
),
# 1.1.2. Changes in time scaling (kappa) & t0 ------------------
# - delay framing: kappa = 1, t0 = t01
# - date framing: kappa = kappa_loss < 1, t0 = t02
date_delay_time_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['kappa_gain'], l['w'], l['a'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['kappa_gain'], u['w'], u['a'], u['t0'], u['t0'])
),
# 1.2.1. Changes in value scaling (mu) ---------------------------
# - delay framing: mu = 1
# - date framing: mu = mu_gain > 1
date_delay_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'],
l['kappa'], l['w'], l['a'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'],
u['kappa'], u['w'], u['a'], u['t0'])
),
# 1.2.2. Changes in value scaling (mu) ---------------------------
# - delay framing: mu = 1, t0 = t01
# - date framing: mu = mu_gain > 1, t0 = t02
date_delay_value_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'],
l['kappa'], l['w'], l['a'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'],
u['kappa'], u['w'], u['a'], u['t0'], u['t0'])
),
# 1.3.1. Changes in time scaling (kappa) -------------------------
# - delay framing: kappa = 1, mu = 1
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1
date_delay_time_and_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'], l['kappa'], l['kappa_gain'], l['w'], l['a'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'], u['kappa'], u['kappa_gain'], u['w'], u['a'], u['t0'])
),
# 1.3.2. Changes in time scaling (kappa) & t0 ------------------
# - delay framing: kappa = 1, mu = 1, t0 = t01
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1, t0 = t02
date_delay_time_and_value_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'], l['kappa'],
l['kappa_gain'], l['w'], l['a'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'], u['kappa'],
u['kappa_gain'], u['w'], u['a'], u['t0'], u['t0'])
),
# 2. Defer-speedup effect ======================================
# 2.1.1. Changes in time scaling (kappa) -----------------------
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals)
# - speedup framing: kappa < 1 (under-responsive to speedups)
defer_speedup_time_scaling =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'], l['a'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'], u['a'], u['t0'])
),
# 2.1.2. Changes in time scaling (kappa) -----------------------
# - neutral framing: kappa = 1, t0 = t01
# - defer framing: kappa > 1 (over-responsive to deferrals), t0 = t02
# - speedup framing: kappa < 1 (under-responsive to speedups), t0 = t03
defer_speedup_time_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'],
l['a'], l['t0'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'],
u['a'], u['t0'], u['t0'], u['t0'])
),
# 2.2.1. Changes in value scaling (mu) ---------------------------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'],
l['kappa'], l['w'], l['a'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'],
u['kappa'], u['w'], u['a'], u['t0'])
),
# 2.2.2. Changes in value scaling (mu) ---------------------------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'],
l['kappa'], l['w'], l['a'], l['t0'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'],
u['kappa'], u['w'], u['a'], u['t0'], u['t0'], u['t0'])
),
# 2.3.1. Changes in time scaling (kappa) and value scaling (mu) -----------------------
# - neutral framing: kappa = 1, mu = 1
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du)
# - speedup framing: kappa < 1 (under-responsive to speedups), mu > 1 (larger du)
defer_speedup_time_and_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'], l['a'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'], u['a'], u['t0'])
),
# 2.3.2. Changes in time scaling (kappa), value scaling (mu), & t0 -----------------------
# - neutral framing: kappa = 1, mu = 1, t0 = t01
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du), t0 = t02
# - speedup framing: kappa < 1 (under-responsive to speedups), mu > 1 (larger du), t0 = t03
defer_speedup_time_and_value_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'],
l['a'], l['t0'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'],
u['a'], u['t0'], u['t0'], u['t0'])
)
),
DFT_C = list(one_condition =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['w'], l['theta_star']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['w'], u['theta_star'])
),
# 1. Date-delay effect =========================================
# 1.1.1. Changes in time scaling (kappa) -------------------------
# - delay framing: kappa = 1
# - date framing: kappa = kappa_loss < 1
date_delay_time_scaling =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['kappa_gain'], l['w'], l['theta_star']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['kappa_gain'], u['w'], u['theta_star'])
),
# 1.2.1. Changes in value scaling (mu) ---------------------------
# - delay framing: mu = 1
# - date framing: mu = mu_gain > 1
date_delay_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'],
l['kappa'], l['w'], l['theta_star']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'],
u['kappa'], u['w'], u['theta_star'])
),
# 1.3.1. Changes in time scaling (kappa) and value scaling (mu)-------------------------
# - delay framing: kappa = 1, mu = 1
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1
date_delay_time_and_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'], l['kappa'],
l['kappa_gain'], l['w'], l['theta_star']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'], u['kappa'],
u['kappa_gain'], u['w'], u['theta_star'])
),
# 2. Defer-speedup effect ======================================
# 2.1.1. Changes in time scaling (kappa) -----------------------
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals)
# - speedup framing: kappa < 1 (under-responsive to speedups)
defer_speedup_time_scaling =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'], l['theta_star']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'], u['theta_star'])
),
# 2.2.1. Changes in value scaling (mu) ---------------------------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'],
l['kappa'], l['w'], l['theta_star']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'],
u['kappa'], u['w'], u['theta_star'])
),
# 2.3.1. Changes in time scaling (kappa) & value scaling (mu) -----------------------
# - neutral framing: kappa = 1, mu = 1
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du)
# - speedup framing: kappa < 1 (under-responsive to speedups), mu > 1 (larger du)
defer_speedup_time_and_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'], l['theta_star']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'], u['theta_star'])
)
),
DFT_CRT = list(one_condition =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['w'], l['theta_star'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['w'], u['theta_star'], u['t0'])
),
# 1. Date-delay effect =========================================
# 1.1.1. Changes in time scaling (kappa) -------------------------
# - delay framing: kappa = 1
# - date framing: kappa = kappa_loss < 1
date_delay_time_scaling =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['kappa_gain'], l['w'], l['theta_star'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['kappa_gain'], u['w'], u['theta_star'], u['t0'])
),
# 1.1.2. Changes in time scaling (kappa) & t0 ------------------
# - delay framing: kappa = 1, t0 = t01
# - date framing: kappa = kappa_loss < 1, t0 = t02
date_delay_time_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['kappa_gain'], l['w'], l['theta_star'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['kappa_gain'], u['w'], u['theta_star'], u['t0'], u['t0'])
),
# 1.2.1. Changes in value scaling (mu) ---------------------------
# - delay framing: mu = 1
# - date framing: mu = mu_gain > 1
date_delay_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'],
l['kappa'], l['w'], l['theta_star'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'],
u['kappa'], u['w'], u['theta_star'], u['t0'])
),
# 1.2.2. Changes in value scaling (mu) ---------------------------
# - delay framing: mu = 1, t0 = t01
# - date framing: mu = mu_gain > 1, t0 = t02
date_delay_value_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'],
l['kappa'], l['w'], l['theta_star'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'],
u['kappa'], u['w'], u['theta_star'], u['t0'], u['t0'])
),
# 1.3.1. Changes in time scaling (kappa) & value scaling (mu) -------------------------
# - delay framing: kappa = 1, mu = 1
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1
date_delay_time_and_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'], l['kappa'],
l['kappa_gain'], l['w'], l['theta_star'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'], u['kappa'],
u['kappa_gain'], u['w'], u['theta_star'], u['t0'])
),
# 1.3.2. Changes in time scaling (kappa), value scaling (mu) & t0 ------------------
# - delay framing: kappa = 1, mu = 1, t0 = t01
# - date framing: kappa = kappa_loss < 1, mu = mu_gain > 1, t0 = t02
date_delay_time_and_value_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['mu_gain'], l['beta'], l['kappa'],
l['kappa_gain'], l['w'], l['theta_star'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_gain'], u['beta'], u['kappa'],
u['kappa_gain'], u['w'], u['theta_star'], u['t0'], u['t0'])
),
# 2. Defer-speedup effect ======================================
# 2.1.1. Changes in time scaling (kappa) -----------------------
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals)
# - speedup framing: kappa < 1 (under-responsive to speedups)
defer_speedup_time_scaling =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'], l['theta_star'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'], u['theta_star'], u['t0'])
),
# 2.1.2. Changes in time scaling (kappa) -----------------------
# - neutral framing: kappa = 1, t0 = t01
# - defer framing: kappa > 1 (over-responsive to deferrals), t0 = t02
# - speedup framing: kappa < 1 (under-responsive to speedups), t0 = t03
defer_speedup_time_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'],
l['theta_star'], l['t0'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'],
u['theta_star'], u['t0'], u['t0'], u['t0'])
),
# 2.2.1. Changes in value scaling (mu) ---------------------------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'],
l['kappa'], l['w'], l['theta_star'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'],
u['kappa'], u['w'], u['theta_star'], u['t0'])
),
# 2.2.2. Changes in value scaling (mu) ---------------------------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'],
l['kappa'], l['w'], l['theta_star'], l['t0'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'],
u['kappa'], u['w'], u['theta_star'], u['t0'], u['t0'], u['t0'])
),
# 2.3.1. Changes in time scaling (kappa) & value scaling (mu) -----------------------
# - neutral framing: kappa = 1, mu = 1
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du)
# - speedup framing: kappa < 1 (under-responsive to speedups), mu > 1 (larger du)
defer_speedup_time_and_value_scaling =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'], l['theta_star'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'], u['theta_star'], u['t0'])
),
# 2.3.2. Changes in time scaling (kappa), value scaling (mu), & t0 -----------------------
# - neutral framing: kappa = 1, mu = 1, t0 = t01
# - defer framing: kappa > 1 (over-responsive to deferrals), mu > 1 (smaller du), t0 = t02
# - speedup framing: kappa < 1 (under-responsive to speedups), mu > 1 (larger du), t0 = t03
defer_speedup_time_and_value_scaling_t0 =
list(lower = c(l['alpha'], l['mu'], l['mu_loss'], l['beta'], l['kappa'],
l['kappa_loss'], l['kappa_gain'], l['w'],
l['theta_star'], l['t0'], l['t0'], l['t0']),
upper = c(u['alpha'], u['mu'], u['mu_loss'], u['beta'], u['kappa'],
u['kappa_loss'], u['kappa_gain'], u['w'],
u['theta_star'], u['t0'], u['t0'], u['t0'])
)
)
)[[model]][[parameterization]][[bound]]
)
}
get_par_pop_stats <- function(model = "DDM", parameterization = "", descstat = "mean") {
means = list('alpha' = 0.94, # Source: DB_2014 (Table 10, Expt. 3)
'mu' = 1, # # Guesstimate
'mu_gain' = 1.2, # Guesstimate
'mu_loss' = 1.2, # Guesstimate
'beta' = 0.75, # DB_JEPG_2014
'kappa' = 1, # Guesstimate
'kappa_gain' = 1 / 1.2, # Guesstimate
'kappa_loss' = 1.2, # Guesstimate
'w' = 0.48, # Source: DB_2014 (Table 10, Expt. 3)
"a" = 1, # Source: guesstimate based on typical range according to rtdists
"t0" = 1.23, # # Source: DB_2014 (Table 10, Expt. 3)
"theta_star" = 1.94 # Source: DB_2014 (Table 10, Expt. 3)
)
sds = list('alpha' = 0.64, # Source: DB_2014 (Table 10, Expt. 3)
'mu' = 1e-6, # # Guesstimate; must be bigger than zero (as diagonals of covariance matrix need to be positive)
'mu_gain' = 0.1, # Guesstimate
'mu_loss' = 0.1, # Guesstimate
'beta' = 0.65, # DB_JEPG_2014
'kappa' = 1e-6, # Guesstimate; must be bigger than zero (as diagonals of covariance matrix need to be positive)
'kappa_gain' = 1 / 1.2, # Guesstimate
'kappa_loss' = 1.2, # Guesstimate
'w' = 0.18, # Source: DB_2014 (Table 10, Expt. 3)
"a" = 1, # Guesstimate
"t0" = 0.14, # # Source: DB_2014 (Table 10, Expt. 3)
"theta_star" = 0.61, # Source: DB_2014 (Table 10, Expt. 3, a = theta_star / 2)
)
M <- unlist(means)
S <- unlist(sds)
get_correlation_mat <- function(names) {
one <- 1
mkp <- 0.10 # Positive correlation of parameters with mu and kappa
mkn <- -0.10 # Negtaive correlation of parameters with mu and kappa
mkw <- 0.85 # Positive correlation within mu-kappa
alb <- 0.85 # Alpha-beta
alw <- 0.25 # Alpha-w
ala <- -0.15 # Alpha-a
alt <- 0.29 # Alpha-t0
b_w <- 0.26 # Beta-w
b_a <- -0.15 # Beta-a
b_t <- 0.23 # Beta-t0
w_a <- -0.15 # w-a
w_t <- 0.04 # w-t0
a_t <- 0.04 # a-t0
mat <-
tibble::frame_matrix(
~alpha, ~mu, ~mu_gain, ~mu_loss, ~beta, ~kappa, ~kappa_gain, ~kappa_loss, ~w, ~a, ~t0,
one, mkp, mkp, mkp, alb, mkp, mkp, mkp, alw, ala, alt,
mkp, one, mkw, mkw, mkp, mkw, mkw, mkw, mkp, mkn, mkp,
mkp, mkw, one, mkw, mkp, mkw, mkw, mkw, mkp, mkn, mkp,
mkp, mkw, mkw, one, mkp, mkw, mkw, mkw, mkp, mkn, mkp,
alb, mkp, mkp, mkp, one, mkp, mkp, mkp, b_w, b_a, b_t,
mkp, mkw, mkw, mkw, mkp, one, mkw, mkw, mkp, mkn, mkp,
mkp, mkw, mkw, mkw, mkp, mkw, one, mkw, mkp, mkn, mkp,
mkp, mkw, mkw, mkw, mkp, mkw, mkw, one, mkp, mkn, mkp,
alw, mkp, mkp, mkp, b_w, mkp, mkp, mkp, one, w_a, w_t,
ala, mkn, mkn, mkn, b_a, mkn, mkn, mkn, w_a, one, a_t,
alt, mkp, mkp, mkp, b_t, mkp, mkp, mkp, w_t, a_t, one
)
selection <- colnames(mat) %in% names
mat[selection, selection]
}
unname(
list(DDM =
list(one_condition =
list(mean = c(M['alpha'], M['mu'], M['beta'], M['kappa'], M['w'], M['a'], M['t0']),
sd = c(S['alpha'], S['mu'], S['beta'], S['kappa'], S['w'], S['a'], S['t0']),
correlation = get_correlation_mat(names = c('alpha', 'mu', 'beta', 'kappa', 'w', 'a', 't0'))
),
# 1. Date-delay effect =========================================
# 1.1.1. Changes in time scaling (kappa) -------------------------
# - delay framing: kappa = 1
# - date framing: kappa = kappa_loss < 1
date_delay_time_scaling =
list(mean = c(M['alpha'], M['mu'], M['beta'], M['kappa'], M['kappa_loss'], M['w'], M['a'], M['t0']),
sd = c(S['alpha'], S['mu'], S['beta'], S['kappa'], S['kappa_loss'], S['w'], S['a'], S['t0']),
correlation = get_correlation_mat(names = c('alpha', 'mu', 'beta', 'kappa', 'kappa_loss', 'w', 'a', 't0'))
),
# 1.1.2. Changes in time scaling (kappa) & t0 ------------------
# - delay framing: kappa = 1, t0 = t01
# - date framing: kappa = kappa_loss < 1, t0 = t01
date_delay_time_scaling_t0 =
list(mean = c(M['alpha'], M['mu'], M['beta'], M['kappa'], M['kappa_loss'], M['w'], M['a'], M['t0'], M['t0']),
sd = c(S['alpha'], S['mu'], S['beta'], S['kappa'], S['kappa_loss'], S['w'], S['a'], S['t0'], S['t0']),
correlation = get_correlation_mat(names = c('alpha', 'mu', 'beta', 'kappa', 'kappa_loss', 'w', 'a', 't0', 't0'))
),
# 1.2.1. Changes in value scaling (mu) ---------------------------
# - delay framing: mu = 1
# - date framing: mu = mu_gain > 1
date_delay_value_scaling =
list(mean = c(M['alpha'], M['mu'], M['mu_gain'], M['beta'], M['kappa'], M['w'], M['a'], M['t0']),
sd = c(S['alpha'], S['mu'], S['mu_gain'], S['beta'], S['kappa'], S['w'], S['a'], S['t0']),
correlation = get_correlation_mat(names = c('alpha', 'mu', 'mu_gain', 'beta', 'kappa', 'w', 'a', 't0'))
),
# 1.2.2. Changes in value scaling (mu) & t0 --------------------
# - delay framing: mu = 1, t0 = t01
# - date framing: mu = mu_gain > 1, t0 = t02
date_delay_value_scaling_t0 =
list(mean = c(M['alpha'], M['mu'], M['mu_gain'], M['beta'], M['kappa'], M['w'], M['a'], M['t0'], M['t0']),
sd = c(S['alpha'], S['mu'], S['mu_gain'], S['beta'], S['kappa'], S['w'], S['a'], S['t0'], S['t0']),
correlation = get_correlation_mat(names = c('alpha', 'mu', 'mu_gain', 'beta', 'kappa', 'w', 'a', 't0', 't0'))
),
# 2. Defer-speedup effect ======================================
# 2.1.1. Changes in time scaling (kappa) -----------------------
# - neutral framing: kappa = 1
# - defer framing: kappa > 1 (over-responsive to deferrals)
# - speedup framing: kappa < 1 (under-responsive to speedups)
defer_speedup_time_scaling =
list(mean = c(M['alpha'], M['mu'], M['beta'], M['kappa'], M['kappa_loss'], M['kappa_gain'], M['w'], M['a'], M['t0']),
sd = c(S['alpha'], S['mu'], S['beta'], S['kappa'], S['kappa_loss'], S['kappa_gain'], S['w'], S['a'], S['t0']),
correlation = get_correlation_mat(names = c('alpha', 'mu', 'beta', 'kappa', 'kappa_loss', 'kappa_gain', 'w', 'a', 't0'))
),
# 2.1.2. Changes in time scaling (kappa) & t0 ------------------
# - neutral framing: kappa = 1, t0 = t01
# - defer framing: kappa > 1 (over-responsive to deferrals), t0 = t02
# - speedup framing: kappa < 1 (under-responsive to speedups), t0 = t03
defer_speedup_time_scaling_t0 =
list(mean = c(M['alpha'], M['mu'], M['beta'], M['kappa'], M['kappa_loss'], M['kappa_gain'], M['w'], M['a'], M['t0'], M['t0'], M['t0']),
sd = c(S['alpha'], S['mu'], S['beta'], S['kappa'], S['kappa_loss'], S['kappa_gain'], S['w'], S['a'], S['t0'], S['t0'], S['t0']),
correlation = get_correlation_mat(names = c('alpha', 'mu', 'beta', 'kappa', 'kappa_loss', 'kappa_gain', 'w', 'a', 't0', 't0', 't0'))
),
# 2.2.1. Changes in value scaling (mu) ---------------------------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling =
list(mean = c(M['alpha'], M['mu'], M['mu_loss'], M['beta'], M['kappa'], M['w'], M['a'], M['t0']),
sd = c(S['alpha'], S['mu'], S['mu_loss'], S['beta'], S['kappa'], S['w'], S['a'], S['t0']),
correlation = get_correlation_mat(names = c('alpha', 'mu', 'mu_loss', 'beta', 'kappa', 'w', 'a', 't0'))
),
# 2.2.2. Changes in value scaling (mu) & t0---------------------
# - neutral framing:
# - defer framing:
# - speedup framing:
defer_speedup_value_scaling_t0 =
list(mean = c(M['alpha'], M['mu'], M['mu_loss'], M['beta'], M['kappa'], M['w'], M['a'], M['t0'], M['t0'], M['t0']),
sd = c(S['alpha'], S['mu'], S['mu_loss'], S['beta'], S['kappa'], S['w'], S['a'], S['t0'], S['t0'], S['t0']),
correlation = get_correlation_mat(names = c('alpha', 'mu', 'mu_loss', 'beta', 'kappa', 'w', 'a', 't0', 't0', 't0'))
)
)
)[[model]][[parameterization]][[descstat]]
)
}
## get_par_values #############################################################
#' Get parameter values
#'
#' @param x parameters
#' @inheritParams fit_model
#' @export
get_par_values = function(x, model = "DFT_C", parameterization = "") {
# 1. Define parameter names ==================================================
parameter_names <-
get_par_names(model = model,
parameterization = "one_condition")
params = list()
names(x) = get_par_names(model = model,
parameterization = parameterization)
# 2. Extract parameter values ================================================
if (model == "DDM") {
if (parameterization == "one_condition") {
params = c(x["alpha"], x["mu"], x["beta"], x["kappa"], x["w"], x["a"], x["t0"])
names(params) = parameter_names
} else if (parameterization == "date_delay_time_scaling") {
# 2.1.1.1. Date/delay effect - changes in time scaling (kappa) ---------------
# delay frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["a"], x["t0"])
# date frame: kappa = kappa_loss < 1
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["a"], x["t0"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_time_scaling_t0") {
# 2.1.1.2. Date/delay effect - changes in time scaling (kappa) & t0 ----------
# delay frame: kappa = 1, t0 = t01
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["a"], x["t01"])
# date frame: kappa = kappa_loss < 1, t0 = t02
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["a"], x["t02"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_value_scaling") {
# 2.1.2.1. Date/delay effect - changes in value scaling (mu) -----------------
# delay frame: mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["a"], x["t0"])
# date frame: mu = mu_gain > 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["a"], x["t0"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_value_scaling_t0") {
# 2.1.2.2. Date/delay effect - changes in value scaling (mu) & t0 ------------
# delay frame: mu = 1, t0 = t01
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["a"], x["t01"])
# date frame: mu = mu_gain > 1, t0 = t02
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["a"], x["t02"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_time_and_value_scaling") {
# 2.1.3.1. Date/delay effect - changes in time scaling (kappa) & value scaling (mu) ---------------
# delay frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["a"], x["t0"])
# date frame: kappa = kappa_loss < 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["a"], x["t0"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_time_and_value_scaling_t0") {
# 2.1.3.2. Date/delay effect - changes in time scaling (kappa), value scaling (mu), & t0 ----------
# delay frame: kappa = 1, t0 = t01
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["a"], x["t01"])
# date frame: kappa = kappa_loss < 1, t0 = t02
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["a"], x["t02"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "defer_speedup_time_scaling") {
# 2.2.1.1. Defer/speedup effect - changes in time scaling (kappa) --------------
# neutral frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["a"], x["t0"])
# defer frame: kappa = 1
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["a"], x["t0"])
# speedup frame: kappa = 1
params[[3]] = c(x["alpha"], x["mu"], x["beta"], x["kappa3"], x["w"], x["a"], x["t0"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_time_scaling_t0") {
# 2.2.1.2. Defer/speedup effect - changes in time scaling (kappa) & t0 -------
# neutral frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["a"], x["t01"])
# defer frame: kappa = 1
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["a"], x["t02"])
# speedup frame: kappa = 1
params[[3]] = c(x["alpha"], x["mu"], x["beta"], x["kappa3"], x["w"], x["a"], x["t03"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_value_scaling") {
# 2.2.2.1. Defer/speedup effect - changes in value scaling (mu) --------------
# neutral frame: mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["a"], x["t0"])
# defer frame: mu > 1 for
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["a"], x["t0"])
# speedup frame: mu > 1 for
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["a"], x["t0"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_value_scaling_t0") {
# 2.2.2.2. Defer/speedup effect - changes in value scaling (mu) & t0 ---------
# neutral frame: mu = 1, t0 = t01
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["a"], x["t01"])
# defer frame: mu > 1 for, t0 = t02
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["a"], x["t02"])
# speedup frame: mu > 1 for, t0 = t03
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["a"], x["t03"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_time_and_value_scaling") {
# 2.2.3.1. Defer/speedup effect - changes in time scaling (kappa) & value scaling (mu) --------------
# neutral frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["a"], x["t0"])
# defer frame: kappa = 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["a"], x["t0"])
# speedup frame: kappa = 1
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa3"], x["w"], x["a"], x["t0"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_time_and_value_scaling_t0") {
# 2.2.3.2. Defer/speedup effect - changes in time scaling (kappa), value scaling (mu) & t0 -------
# neutral frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["a"], x["t01"])
# defer frame: kappa = 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["a"], x["t02"])
# speedup frame: kappa = 1
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa3"], x["w"], x["a"], x["t03"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
}
} else if (model == "DFT_C") {
if (parameterization == "one_condition") {
params = c(x["alpha"], x["mu"], x["beta"], x["kappa"], x["w"], x["theta_star"])
names(params) = parameter_names
} else if (parameterization == "date_delay_time_scaling") {
# 2.1.1. Date/delay effect - changes in time scaling (kappa) ---------------
# delay frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["theta_star"])
# date frame: kappa = kappa_loss < 1
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["theta_star"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_value_scaling") {
# 2.1.2. Date/delay effect - changes in value scaling (mu) -----------------
# delay frame: mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["theta_star"])
# date frame: mu = mu_gain > 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["theta_star"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_time_and_value_scaling") {
# 2.1.3. Date/delay effect - changes in time scaling (kappa) & value scaling (mu) ---------------
# delay frame: kappa = 1, mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["theta_star"])
# date frame: kappa = kappa_loss < 1, mu = mu_gain > 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["theta_star"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "defer_speedup_time_scaling") {
# 2.2.1. Defer/speedup effect - changes in time scaling (kappa) --------------
# neutral frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["theta_star"])
# defer frame: kappa = 1
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["theta_star"])
# speedup frame: kappa = 1
params[[3]] = c(x["alpha"], x["mu"], x["beta"], x["kappa3"], x["w"], x["theta_star"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_value_scaling") {
# 2.2.2. Defer/speedup effect - changes in value scaling (mu) --------------
# neutral frame: mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["theta_star"])
# defer frame: mu > 1 for
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["theta_star"])
# speedup frame: mu > 1 for
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["theta_star"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_time_and_value_scaling") {
# 2.2.3. Defer/speedup effect - changes in time scaling (kappa) & value scaling (mu) --------------
# neutral frame: kappa = 1, mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["theta_star"])
# defer frame: kappa = 1, mu > 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["theta_star"])
# speedup frame: kappa = 1, mu > 1
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa3"], x["w"], x["theta_star"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
}
} else if (model == "DFT_CRT") {
if (parameterization == "one_condition") {
params = c(x["alpha"], x["mu"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t0"])
names(params) = parameter_names
} else if (parameterization == "date_delay_time_scaling") {
# 2.1.1.1. Date/delay effect - changes in time scaling (kappa) ---------------
# delay frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["theta_star"], x["t0"])
# date frame: kappa = kappa_loss < 1
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["theta_star"], x["t0"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_time_scaling_t0") {
# 2.1.1.2. Date/delay effect - changes in time scaling (kappa) & t0 ----------
# delay frame: kappa = 1, t0 = t01
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["theta_star"], x["t01"])
# date frame: kappa = kappa_loss < 1, t0 = t02
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["theta_star"], x["t02"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_value_scaling") {
# 2.1.2.1. Date/delay effect - changes in value scaling (mu) -----------------
# delay frame: mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t0"])
# date frame: mu = mu_gain > 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t0"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_value_scaling_t0") {
# 2.1.2.2. Date/delay effect - changes in value scaling (mu) & t0 ------------
# delay frame: mu = 1, t0 = t01
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t01"])
# date frame: mu = mu_gain > 1, t0 = t02
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t02"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_time_and_value_scaling") {
# 2.1.3.1. Date/delay effect - changes in time scaling (kappa) & value scaling (mu) ---------------
# delay frame: kappa = 1, mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["theta_star"], x["t0"])
# date frame: kappa = kappa_loss < 1, mu = mu_gain > 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["theta_star"], x["t0"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "date_delay_time_and_value_scaling_t0") {
# 2.1.3.2. Date/delay effect - changes in time scaling (kappa), value scaling (mu) & t0 ----------
# delay frame: kappa = 1, mu = 1, t0 = t01
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["theta_star"], x["t01"])
# date frame: kappa = kappa_loss < 1, mu = mu_gain > 1, t0 = t02
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["theta_star"], x["t02"])
names(params[[1]]) = names(params[[2]]) = parameter_names
} else if (parameterization == "defer_speedup_time_scaling") {
# 2.2.1.1. Defer/speedup effect - changes in time scaling (kappa) --------------
# neutral frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["theta_star"], x["t0"])
# defer frame: kappa = 1
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["theta_star"], x["t0"])
# speedup frame: kappa = 1
params[[3]] = c(x["alpha"], x["mu"], x["beta"], x["kappa3"], x["w"], x["theta_star"], x["t0"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_time_scaling_t0") {
# 2.2.1.2. Defer/speedup effect - changes in time scaling (kappa) & t0 -------
# neutral frame: kappa = 1
params[[1]] = c(x["alpha"], x["mu"], x["beta"], x["kappa1"], x["w"], x["theta_star"], x["t01"])
# defer frame: kappa = 1
params[[2]] = c(x["alpha"], x["mu"], x["beta"], x["kappa2"], x["w"], x["theta_star"], x["t02"])
# speedup frame: kappa = 1
params[[3]] = c(x["alpha"], x["mu"], x["beta"], x["kappa3"], x["w"], x["theta_star"], x["t03"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_value_scaling") {
# 2.2.2.1. Defer/speedup effect - changes in value scaling (mu) --------------
# neutral frame: mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t0"])
# defer frame: mu > 1 for
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t0"])
# speedup frame: mu > 1 for
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t0"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_value_scaling_t0") {
# 2.2.2.2. Defer/speedup effect - changes in value scaling (mu) & t0 ---------
# neutral frame: mu = 1, t0 = t01
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t01"])
# defer frame: mu > 1 for, t0 = t02
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t02"])
# speedup frame: mu > 1 for, t0 = t03
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa"], x["w"], x["theta_star"], x["t03"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_time_and_value_scaling") {
# 2.2.3.1. Defer/speedup effect - changes in time scaling (kappa) --------------
# neutral frame: kappa = 1, mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["theta_star"], x["t0"])
# defer frame: kappa = 1, mu > 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["theta_star"], x["t0"])
# speedup frame: kappa = 1, mu > 1
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa3"], x["w"], x["theta_star"], x["t0"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
} else if (parameterization == "defer_speedup_time_and_value_scaling_t0") {
# 2.2.3.2. Defer/speedup effect - changes in time scaling (kappa) & t0 -------
# neutral frame: kappa = 1, mu = 1
params[[1]] = c(x["alpha"], x["mu1"], x["beta"], x["kappa1"], x["w"], x["theta_star"], x["t01"])
# defer frame: kappa = 1, mu > 1
params[[2]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa2"], x["w"], x["theta_star"], x["t02"])
# speedup frame: kappa = 1, mu > 1
params[[3]] = c(x["alpha"], x["mu2"], x["beta"], x["kappa3"], x["w"], x["theta_star"], x["t03"])
names(params[[1]]) = names(params[[2]]) = names(params[[3]]) = parameter_names
}
}
# 3. Return parameter values =================================================
params
}
## itch_ddm ###################################################################
#' Intertemporal choice drift diffusion model
#'
#' @param stimuli bla # TODO: use @inher@inheritParams
#' @param parameters ble # TODO: use @inher@inheritParams
#' @export
#' This model is inspired by Dai & Busemeyer, J Exp Psychol Gen, 2014 and Scholten & Read, J Exp Psychol Learn Mem Cogn, 2013
itch_ddm <- function(stimuli, parameters, parameterization = "", frame = "", n = 1, sim_fun = 'rtdists_package') {
# 1. Unlist parameters =======================================================
x = unlist(parameters)
# 2. Compute drift rate ======================================================
v <- compute_drift_rate(parameters = parameters,
stimuli = stimuli,
parameterization = parameterization,
frame = frame)
# 3. Compute and return predicted responses ==================================
# Either using Dai & Busemeyer's method (Eq. 9 in their paper) or using the rtdists packages
if (sim_fun == 'dai_busemeyer') {
tibble::tibble(p_ll = itchmodel::db_bin_choice_prob(d = v,
s = 0.1,
a = unname(x['a']),
z = 0))
} else if (sim_fun == 'rtdists_package') {
purrr::map_df(.x = v,
.f = rtdists::rdiffusion,
n = n,
s = 1, # This is the default in the rtdists package, see documentation
a = x['a'],
t0 = x['t0']
) %>%
dplyr::mutate(v = v)
}
}
## ll_dft ################################################################
#' Log-likelihood for decision field theory
#'
#' This is based on Junyi Dai's intertemporal choice code
#'
#' @export
ll_dft <- function(x, stimuli, frame = "", observations, rt = TRUE) {
# 1. Unlist parameters =======================================================
x <- unlist(x)
# 2. Compute drift rate ======================================================
du <- compute_transformation_diffs(parameters = x,
stimuli = stimuli,
parameterization = parameterization,
frame = frame,
variable = 'du')
dp <- compute_transformation_diffs(parameters = x,
stimuli = stimuli,
parameterization = parameterization,
frame = frame,
variable = 'dp')
d <- unname(x["w"] * du - (1 - x["w"]) * dp)
# Probability of sampling from or attending to a specific attribute at any
# time is assumed to equal the corresponding attention weight. See Eq. 15 in
# Dai & Busemeyer
s = unname(sqrt(x["w"] * du^2 + (1 - x["w"]) * dp^2 - d^2))
# 3. Compute densities =======================================================
if (rt) {
densities <-
tryCatch(rtdists::ddiffusion(rt = observations$rt,
response = observations$response,
v = d,
s = s,
a = as.double(x["theta_star"]) * s,
t0 = as.double(x['t0'])),
error = function(e) 0)
} else {
densities <-
# dft_cp(d = d,
# s = s,
# theta = unname(x["theta_star"] * s),
# z = 0,
# response = observations$response)
tryCatch(dft_cp(d = d,
s = s,
theta = unname(x["theta_star"] * s),
z = 0,
response = observations$response),
error = function(e) 0)
}
# Replace any NA values to prevent that the next line returns an error
densities[is.na(densities)] <- 0
if (any(densities == 0)) return(1e6)
return(-sum(log(densities)))
}
## ll_diffusion ################################################################
#' Log-likelihood for drift diffusion model
#'
#' This is based on example in rtdists::ddiffusion package
#'
#' @param pars Vector of parameters
#' @inheritParams rtdists::ddiffusion
#' @export
ll_diffusion <- function(x, stimuli, frame = "", observations) {
# 1. Unlist parameters =======================================================
x <- unlist(x)
# 2. Compute drift rate ======================================================
v <- compute_drift_rate(x, stimuli,
frame = frame)
# 3. Compute densities =======================================================
# Compute densities
densities <-
tryCatch(rtdists::ddiffusion(rt = observations$rt,
response = observations$response,
v = v,
s = 1,
a = x['a'],
t0 = x['t0']),
error = function(e) 0)
# Replace any NA values to prevent that the next line returns an error
densities[is.na(densities)] <- 0
if (any(densities == 0)) return(1e6)
return(-sum(log(densities)))
}
## logistic ####################################################################
#' Logistic function
#'
#' @param x data
#' @export
logistic <- function(x) {
1 / (1 + exp(-x))
}
## make_stimulus_df ############################################################
#' Make stimulus list (i.e. trials)
#'
#' @param frames character vector of frame names, defaults to c('delay', 'date')
#' @param m_l double vector of monetary amounts of larger option (in Euros), defaults to c(20)
#' @param pct_m_l double vector of monetary amount of smaller option, experessed as percentage of larger option, defaults to c(0.2, 0.5, 0.7, 0.8, 0.88, 0.93, 0.97, 0.99)
#' @param t_s double vector of delay of reward delivery for sooner option (in days), defaults to c(0)
#' @param interval double vector of intervals between sooner and later option (in days), defaults to c(0)
#' @param n_reps number of repetitions of each combination of parameters
#' @return The output is a tibble with columns
make_stimulus_df <- function(frames = "neutral",
m_l = c(20),
pct_m_l = c(0.2, 0.5, 0.7, 0.8, 0.88, 0.93, 0.97, 0.99),
t_s = c(0),
interval = c(1,2,3,5,7,10,14),
n_reps = 1) {
if (all(c('delay', 'date') %in% frames)) {
factor_levels <- c('delay','date')
} else if (all(c('neutral', 'defer', 'speedup') %in% frames)) {
factor_levels <- c('neutral', 'defer', 'speedup')
} else { # e.g. for indifference point procedure
factor_levels <- 'neutral'
}
expand.grid(frame = factor(frames, levels = factor_levels),
pct_m_l = pct_m_l,
m_l = m_l,
t_s = t_s,
interval = interval,
rep = 1:n_reps
) %>%
dplyr::mutate(m_s = pct_m_l * m_l,
t_l = t_s + interval) %>%
dplyr::select(frame, m_s, t_s, m_l, t_l, dplyr::everything()) %>%
dplyr::select(-rep) %>%
tibble::as.tibble()
}
## nest_join_stimuli_parameters ################################################
#' Nest and join stimulu and parameters
#'
#' nest_join_stimuli_parameters
#'
nest_join_stimuli_parameters <- function(stimuli = make_stimulus_df(),
parameters, parameterization = "") {
# Assertions -----------------------------------------------------------------
assertthat::assert_that(parameterization %in% c("one_condition",
"date_delay_time_scaling",
"date_delay_time_scaling_t0",
"date_delay_time_and_value_scaling",
"date_delay_time_and_value_scaling_t0",
"date_delay_value_scaling",
"date_delay_value_scaling_t0",
"defer_speedup_time_scaling",
"defer_speedup_time_scaling_t0",
"defer_speedup_time_and_value_scaling",
"defer_speedup_time_and_value_scaling_t0",
"defer_speedup_value_scaling",
"defer_speedup_value_scaling_t0"
),
msg = "Specify parameterization as \n\"one-condition\" (i.e., estimate parameters from a single experimental condition), \n\"date_delay_value_scaling\" (i.e., estimate free value function scaling parameter across two experimental conditions), \n\"date_delay_value_scaling_t0\" (i.e., estimate free value function scaling parameter and t0 across two experimental conditions), \n\"date_delay_time_scaling\" (i.e., estimate free time function sensitivity parameter across two experimental conditions), \n\"date_delay_time_scaling_t0\" (i.e., estimate free time function sensitivity parameter and t0 across two experimental conditions), \n\"date_delay_time_and_value_scaling\", \n\"date_delay_time_and_value_scaling_t0\", \n\"defer_speedup_value_scaling\" (i.e., estimate free value function scaling parameter across two experimental conditions), \n\"defer_speedup_value_scaling_t0\" (i.e., estimate free value function scaling parameter and t0 across two experimental conditions), \n\"defer_speedup_time_scaling\" (i.e., estimate free time function scaling parameter across two experimental conditions), or \n\"defer_speedup_time_scaling_t0\" (i.e., estimate free time function scaling parameter and t0 across two experimental conditions), \n\"defer_speedup_time_and_value_scaling\", \n\"defer_speedup_time_and_value_scaling_t0\"."
)
# Nested stimuli -------------------------------------------------------------
nested_stimuli <-
stimuli %>%
dplyr::group_by(frame) %>%
tidyr::nest(.key = 'stimuli')
# Nested parameters ----------------------------------------------------------
nested_parameters <-
parameters %>%
get_par_values(model = "DDM", parameterization = parameterization) %>%
dplyr::bind_rows() %>%
dplyr::mutate(frame = factor(levels(stimuli$frame),
levels = levels(stimuli$frame))) %>%
dplyr::select(frame, dplyr::everything()) %>%
dplyr::group_by(frame) %>%
tidyr::nest(.key = 'parameters')
# This is the old code:
# nested_parameters <-
# parameters %>%
# tibble::as.tibble() %>%
# dplyr::mutate(frame = factor(levels(stimuli$frame),
# levels = levels(stimuli$frame))) %>%
# dplyr::select(frame, dplyr::everything()) %>%
# dplyr::group_by(frame) %>%
# tidyr::nest(.key = 'parameters')
# Join nested stimuli and parameters -----------------------------------------
nested_stimuli %>%
dplyr::left_join(nested_parameters, by = "frame")
}
## params_to_tibble ############################################################
#' Convert vector of parameters to a tibble with named columns
#'
#' @param x vector or parameters
#' @param model model name
#' @param parameterization model parameterization
params_to_tibble <- function(x, model = "", parameterization = "") {
par_names <-
itchmodel::get_par_names(model = model,
parameterization = parameterization)
params <-
matrix(x,
ncol = length(par_names)
)
colnames(params) <-
par_names
params %>%
tibble::as.tibble()
}
## sample_params_from_pop ######################################################
#' Sample model parameters from population
#'
#' @export
sample_params_from_pop <- function(n, model = "DDM", parameterization = "", stimuli) {
# Get parameter population means, standard deviations, and correlations ------
M <- get_par_pop_stats(model = model,
parameterization = parameterization,
descstat = "mean")
S <- get_par_pop_stats(model = model,
parameterization = parameterization,
descstat = "sd")
correlation_mat <- get_par_pop_stats(model = model,
parameterization = parameterization,
descstat = "correlation")
# Parameter covariance matrix ------------------------------------------------
stdev_ij <- S %*% t(S)
cov_mat <- correlation_mat * stdev_ij
# Sample parameters from truncated multivariate normal -----------------------
# constrained by lower and upper bounds + linear inequalities
iter <- 1
params <- double()
par_names <-
itchmodel::get_par_names(model = model,
parameterization = parameterization)
print(sprintf("Parameterization: %s", parameterization))
while (iter <= 50) {
# Sample parameters with bound constraints
params <-
rbind(params,
tmvtnorm::rtmvnorm2(n = n,
mean = M,
sigma = cov_mat,
lower = get_par_bounds(model = model,
parameterization = parameterization,
bound = "lower"),
upper = get_par_bounds(model = model,
parameterization = parameterization,
bound = "upper") + 1e-6,
)
)
# Check linear inequalities
params_list <-
purrr::array_branch(array = params,
margin = 1)
nested_stimuli <-
stimuli %>%
dplyr::group_by(frame) %>%
tidyr::nest() %>%
dplyr::rename(stimuli = data)
lineqs <-
purrr::pmap_dbl(.l = list(x = params_list),
.f = get_nonlinear_constraints,
data = nested_stimuli,
model = model,
parameterization = parameterization
)
# Only retain parameter samples that meet linear inequality constraints
params <- params[lineqs <= 0,]
# Add column names
params <- matrix(params, ncol = length(par_names))
colnames(params) <- par_names
if (nrow(params) >= n) {
params <- params[1:n,]
print(sprintf('Iteration %d: Sampled %d out of %d valid parameter sets', iter, nrow(params), n))
break
} else {
print(sprintf('Iteration %d: Sampled %d out of %d valid parameter sets', iter, nrow(params), n))
iter <- iter + 1
}
}
# Add column names
params <- matrix(params, ncol = length(par_names))
colnames(params) <- par_names
assertthat::assert_that(nrow(params) >= n,
msg = "Algorithm failed to sampled enough valid parameter sets. Consider adjusting settings.")
# Add index & parameterization column and return
params <-
params %>%
tibble::as.tibble() %>%
dplyr::mutate(ix = 1:n())
params %>%
dplyr::group_by(ix) %>%
tidyr::nest()
}
## sim_data ####################################################################
#' Simulate intertemporal choice task performance
#'
#' sim_data
#'
#' @param stimuli tibble with trial stimuli
#' @param parameters list with model parameters
#' @param parameterization char vectore specifying model parameterization, defaults to "date_delay_time_scaling"
sim_data <- function(stimuli = make_stimulus_df(),
parameters, parameterization = "date_delay_time_scaling", n = 1, sim_fun = 'rtdists_package') {
# Simulate observations
nest_join_stimuli_parameters(stimuli = stimuli,
parameters = parameters,
parameterization = parameterization) %>%
dplyr::mutate(observations = purrr::pmap(.l = list(stimuli = stimuli,
parameters = parameters,
parameterization = parameterization,
frame = as.character(frame)),
.f = itch_ddm,
n = n,
sim_fun = sim_fun),
model = parameterization
)
}
# sim_n_datasets(stimuli = make_stimulus_df(),
# parameters, parameterization = "date_delay_time_scaling", n_dataset = 1, n_trial_per_dataset = 1, sim_fun = 'rtdists_package')) {
#
#
#
#
# sim_data <- function(stimuli = make_stimulus_df(),
# parameters,
# parameterization = parameterization,
# n = 1,
# sim_fun = sim_fun)
#
#
#
#
#
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.