inst/doc/customize_ddms.R

In dRiftDM: Estimating (Time-Dependent) Drift Diffusion Models

## ----echo = FALSE, message=FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(dRiftDM)
set.seed(1014)

## -----------------------------------------------------------------------------
a_model <- ratcliff_dm() # some model
names(a_model) # the entries of the underlying list

## ----echo=FALSE, out.width = '90%'--------------------------------------------
knitr::include_graphics("evaluation_chain.png")

## -----------------------------------------------------------------------------
# Create a model (here the Diffusion Model for Conflict Tasks, DMC)
ddm <- dmc_dm()
flex_prms(ddm)

## -----------------------------------------------------------------------------
ddm_free_muc <- modify_flex_prms(
  ddm, # the model
  instr = "muc ~" # an instructions in a formula-like format
)
flex_prms(ddm_free_muc)

## -----------------------------------------------------------------------------
coef(ddm) # in this model muc is the same for all conditions

coef(ddm_free_muc) # here muc can be different for the conditions

coef(ddm_free_muc)[1] <- 5
coef(ddm_free_muc)

## -----------------------------------------------------------------------------
conds(ddm)

## -----------------------------------------------------------------------------
conds(ddm) <- c("comp", "neutral", "incomp")

## -----------------------------------------------------------------------------
flex_prms(ddm)

## -----------------------------------------------------------------------------
instructions <- "
a <!>                      # 'a' is fixed across all conditions
A ~ incomp == -(A ~ comp)  # A in incomp is -A in comp
A <!> neutral              # A is fixed for the neutral condition
A ~ neutral => 0           # A is zero for the neutral condition
"

ddm <- modify_flex_prms(
  object = ddm,
  instr = instructions
)
print(ddm)

## -----------------------------------------------------------------------------
ddm <- dmc_dm() # some pre-built model
names(comp_funs(ddm))

## -----------------------------------------------------------------------------
all_funs <- component_shelf()
names(all_funs)

## ----eval = F-----------------------------------------------------------------
# ... <- function(prms_model, prms_solve, t_vec, one_cond, ddm_opts) {
#   ...
# }

## -----------------------------------------------------------------------------
flex_prms(ddm)

## -----------------------------------------------------------------------------
prms_solve(ddm)

## -----------------------------------------------------------------------------
t_max <- prms_solve(ddm)["t_max"]
dt <- prms_solve(ddm)["dt"]
t_vec <- seq(0, t_max, dt)
head(t_vec)
tail(t_vec)

## ----eval = F-----------------------------------------------------------------
# ... <- function(prms_model, prms_solve, x_vec, one_cond, ddm_opts) {
#   ...
# }

## -----------------------------------------------------------------------------
dx <- prms_solve(ddm)["dx"]
x_vec <- seq(-1, 1, dx)

## -----------------------------------------------------------------------------
cust_mu <- function(prms_model, prms_solve, t_vec, one_cond, ddm_opts) {
  # extract all parameters (one row of the parameter matrix)
  muc <- prms_model[["muc"]]
  mua <- prms_model[["mua"]]
  sign <- prms_model[["sign"]]

  # and return the drift rate at each time step
  mu <- rep(muc + sign * mua, length(t_vec))
  return(mu)
}

## -----------------------------------------------------------------------------
cust_model <- function() {
  # define all parameters and conditions
  prms_model <- c(
    muc = 3, mua = 1, sign = 1, # parameters for the custom drift rate function
    b = .6, # parameter for a time-independent boundary "b"
    non_dec = .2 # parameter for a non-decision time "non_dec"
  )
  conds <- c("comp", "incomp")

  # get access to pre-built component functions
  comps <- component_shelf()

  # call the drift_dm function which is the backbone of dRiftDM
  ddm <- drift_dm(
    prms_model = prms_model,
    conds = conds,
    subclass = "my_custom_model",
    mu_fun = cust_mu, # your custom drift rate function
    mu_int_fun = comps$dummy_t, # a dummy function, because no integral
    # of the drift rate is required per default
    x_fun = comps$x_dirac_0, # pre-built dirac delta on zero for the starting point
    b_fun = comps$b_constant, # pre-built time-independent boundary with parameter b
    dt_b_fun = comps$dt_b_constant, # pre-built derivative of the boundary
    nt_fun = comps$nt_constant # pre-built non-decision time with parameter non_dec
  )

  # modify the flex_prms object to achieve the desired behavior of 'sign'
  # -> don't consider 'sign' a free parameter to estimate and set it to -1
  # for incompatible conditions
  ddm <- modify_flex_prms(ddm, instr = "sign <!>
                                        sign ~ incomp => -1")

  return(ddm)
}

ddm <- cust_model()

## -----------------------------------------------------------------------------
cust_mu_int <- function(prms_model, prms_solve, t_vec, one_cond, ddm_opts) {
  # extract all parameters (one row of the parameter matrix)
  muc <- prms_model[["muc"]]
  mua <- prms_model[["mua"]]
  sign <- prms_model[["sign"]]

  # and return the integral of the drift rate
  mu <- (muc + sign * mua) * t_vec
  return(mu)
}

## ----eval = F, echo = F-------------------------------------------------------
# cust_model <- function() {
#   # define all parameters and conditions
#   prms_model <- c(
#     muc = 3, mua = 1, sign = 1, # prms for the custom drift rate function
#     b = .6, # parameter for a time-independent boundary b
#     non_dec = .2
#   ) # parameter for a non-decision time
#   conds <- c("comp", "incomp")
# 
#   # get access to pre-build component functions
#   comps <- component_shelf()
# 
#   # call the drift_dm function which is the backbone of dRiftDM
#   ddm <- drift_dm(
#     prms_model = prms_model,
#     conds = conds,
#     subclass = "my_custom_model",
#     mu_fun = cust_mu, # custom defined
#     mu_int_fun = cust_mu_int, # to test :)
#     x_fun = comps$x_dirac_0, # dirac delta on zero for the starting point
#     b_fun = comps$b_constant, # time-independent boundary with parameter b
#     dt_b_fun = comps$dt_b_constant, # respective derivative of the boundary
#     nt_fun = comps$nt_constant # non-decision time with parameter non_dec
#   )
# 
#   # modify the flex_prms object to achieve the desired behavior of 'sign'
#   # -> don't consider 'sign' a free parameter to estimate and set it to -1
#   # for incompatible conditions
#   ddm <- modify_flex_prms(
#     ddm,
#     instr = "sign <!>
#              sign ~ incomp => -1"
#   )
# 
#   return(ddm)
# }
# 
# ddm <- cust_model()
# plot(re_evaluate_model(ddm)$pdfs$comp$pdf_u)
# solver(ddm) <- "im_zero"
# points(re_evaluate_model(ddm)$pdfs$comp$pdf_u, col = "red", ty = "l")

## -----------------------------------------------------------------------------
cust_x <- function(prms_model, prms_solve, x_vec, one_cond, ddm_opts) {
  dx <- prms_solve[["dx"]]
  z <- prms_model[["z"]]
  stopifnot(z > 0, z < 1)

  # create a dirac delta for the starting point
  x <- numeric(length = length(x_vec))
  index <- round(z * (length(x) - 1)) + 1
  x[index] <- 1 / dx # make sure it integrates to 1
  return(x)
}

## -----------------------------------------------------------------------------
cust_model <- function() {
  # define all parameters and conditions
  prms_model <- c(
    z = .75, # parameter for the custom starting point
    muc = 4, # parameter for a time-independent drift rate "muc"
    b = .6, # parameter for a time-independent boundary "b"
    non_dec = .2 # parameter for a non-decision time "non_dec"
  )

  # each model must have a condition (in this example it is not relevant, so we just call it "foo")
  conds <- c("foo")

  # get access to pre-built component functions
  comps <- component_shelf()

  # call the drift_dm function which is the backbone of dRiftDM
  ddm <- drift_dm(
    prms_model = prms_model,
    conds = conds,
    subclass = "my_custom_model",
    mu_fun = comps$mu_constant, # time-independent drift rate with parameter muc
    mu_int_fun = comps$mu_int_constant, # respective integral of the drift rate
    x_fun = cust_x, # custom starting point function with parameter z
    b_fun = comps$b_constant, # time-independent boundary with parameter b
    dt_b_fun = comps$dt_b_constant, # respective derivative of the boundary
    nt_fun = comps$nt_constant # non-decision time with parameter non_dec
  )
  return(ddm)
}

ddm <- cust_model()

## -----------------------------------------------------------------------------
# the boundary function
cust_b <- function(prms_model, prms_solve, t_vec, one_cond, ddm_opts) {
  b0 <- prms_model[["b0"]]
  kappa <- prms_model[["kappa"]]
  t05 <- prms_model[["t05"]]

  return(b0 * (1 - kappa * t_vec / (t_vec + t05)))
}

## -----------------------------------------------------------------------------
# the derivative of the boundary function
cust_dt_b <- function(prms_model, prms_solve, t_vec, one_cond, ddm_opts) {
  b0 <- prms_model[["b0"]]
  kappa <- prms_model[["kappa"]]
  t05 <- prms_model[["t05"]]

  return(-(b0 * kappa * t05) / (t_vec + t05)^2)
}

## -----------------------------------------------------------------------------
cust_model <- function() {
  # define all parameters and conditions
  prms_model <- c(
    b0 = .6, # parameters for the custom boundary
    kappa = .5,
    t05 = .15,
    muc = 4, # parameter for a time-independent drift rate "muc"
    non_dec = .2 # parameter for a non-decision time "non_dec"
  )

  # each model must have a condition (in this example it is not relevant, so we just call it "foo")
  conds <- c("foo")

  # get access to pre-built component functions
  comps <- component_shelf()

  # call the drift_dm function which is the backbone of dRiftDM
  ddm <- drift_dm(
    prms_model = prms_model,
    conds = conds,
    subclass = "my_custom_model",
    mu_fun = comps$mu_constant, # time-independent drift rate with parameter muc
    mu_int_fun = comps$mu_int_constant, # respective integral of the drift rate
    x_fun = comps$x_dirac_0, # dirac delta on zero
    b_fun = cust_b, # custom time-dependent boundary
    dt_b_fun = cust_dt_b, # respective derivative of the boundary
    nt_fun = comps$nt_constant # non-decision time with parameter non_dec
  )
  return(ddm)
}
ddm <- cust_model()

## ----eval = F, echo = F-------------------------------------------------------
# plot(
#   simulate_traces(ddm, 1, sigma = 0)
# )

## -----------------------------------------------------------------------------
cust_nt <- function(prms_model, prms_solve, t_vec, one_cond, ddm_opts) {
  # get the relevant parameters
  non_dec <- prms_model[["non_dec"]]
  range_non_dec <- prms_model[["range_non_dec"]]

  # get the settings for the time space discretization
  t_max <- prms_solve[["t_max"]]
  dt <- prms_solve[["dt"]]

  # calculate the density
  d_nt <- dunif(x = t_vec, min = non_dec - range_non_dec / 2, max = non_dec + range_non_dec / 2)
  d_nt <- d_nt / (sum(d_nt) * dt) # ensure it integrates to 1
  return(d_nt)
}

## -----------------------------------------------------------------------------
cust_model <- function() {
  # define all parameters and conditions
  prms_model <- c(
    non_dec = .2, # parameters for the custom non-decision time
    range_non_dec = .05,
    muc = 4, # parameter for a time-independent drift rate
    b = .6 # parameter for a time-independent boundary
  )

  # each model must have a condition (in this example it is not relevant, so we just call it "foo")
  conds <- c("foo")

  # get access to pre-built component functions
  comps <- component_shelf()

  # call the drift_dm function which is the backbone of dRiftDM
  ddm <- drift_dm(
    prms_model = prms_model,
    conds = conds,
    subclass = "my_custom_model",
    mu_fun = comps$mu_constant, # time-independent drift rate with parameter muc
    mu_int_fun = comps$mu_int_constant, # respective integral of the drift rate
    x_fun = comps$x_dirac_0, # dirac delta on zero
    b_fun = comps$b_constant, # time-independent boundary b
    dt_b_fun = comps$dt_b_constant, # respective derivative of the boundary
    nt_fun = cust_nt # custom non-decision time
  )
  return(ddm)
}
ddm <- cust_model()

## ----eval = F, echo = F-------------------------------------------------------
# plot(ddm)

## ----eval = F-----------------------------------------------------------------
# ?component_shelf

## ----eval = F-----------------------------------------------------------------
# all_funs <- component_shelf() # all functions
# View(all_funs$x_uniform) # see the code for x_uniform

## -----------------------------------------------------------------------------
ddm <- dmc_dm()
test <- simulate_traces(ddm, k = 1, sigma = 0, add_x = FALSE)
plot(test)

## ----eval = F, echo = T-------------------------------------------------------
# plot(ddm, col = c("green", "red"))

## ----eval = T, echo = F-------------------------------------------------------
plot(ddm, mfrow = c(1, 1), col = c("green", "red")) # otherwise the vignette won't build

## -----------------------------------------------------------------------------
cust_mu <- function(prms_model, prms_solve, t_vec, one_cond, ddm_opts) {
  print(ddm_opts) # print out the values of ddm_opts
  muc <- rep(prms_model[["muc"]], length(t_vec))
  return(muc)
}

a_model <- ratcliff_dm() # a dummy model for demonstration purpose
ddm_opts(a_model) <- "Hello World" # attach "Hello World" to the model
comp_funs(a_model)[["mu_fun"]] <- cust_mu # swap in the custom component function

# evaluate the model (which will evaluate the custom drift rate function with the user-defined R
# object for ddm_opts)
a_model <- re_evaluate_model(a_model)