inst/doc/model_recovery.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)

## -----------------------------------------------------------------------------
# 1.1) create the Ratcliff model and set new conditions
ratcliff_model <- ratcliff_dm()
conds(ratcliff_model) <- c("comp", "incomp")

# allow the drift rate to vary
instructions <- "muc ~ "

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



# 1.2) now create DMC (without trial-by-trial variability in the starting point
# and non-decision time; matching with the Ratcliff model)
dmc_model <- dmc_dm(var_non_dec = FALSE, var_start = FALSE)


# 1.3) create the parameter space for generating synthetic data
lower_dmc <- c(muc = 2, b = 0.3, non_dec = 0.2, tau = 0.02, A = 0.02)
upper_dmc <- c(muc = 6, b = 0.8, non_dec = 0.4, tau = 0.15, A = 0.20)

lower_ratcliff <- c(muc = 2, b = 0.3, non_dec = 0.2)
upper_ratcliff <- c(muc = 6, b = 0.8, non_dec = 0.4)

## -----------------------------------------------------------------------------
coef(ratcliff_model)
coef(dmc_model)

## -----------------------------------------------------------------------------
set.seed(1)
# Data generated under the Ratcliff model
data_ratcliff <- simulate_data(
  ratcliff_model, #  the model
  n = 200, # choose this value reasonably
  k = 100, # more data to filter afterwards
  lower = lower_ratcliff, # the lower parameter space for simulation
  upper = upper_ratcliff # the upper parameter space for simulation
)
# choose data were the drift rate in incompatible conditions is
# reasonably larger than in compatible conditions
prms_ratcliff <- data_ratcliff$prms
diffs <- prms_ratcliff$muc.comp - prms_ratcliff$muc.incomp
filter <- which(diffs < 1.5 & diffs > 0.1)
ids <- prms_ratcliff$ID[filter[1:10]]
data_ratcliff <- data_ratcliff$synth_data # extract the data.frame of raw RTs
data_ratcliff <- data_ratcliff[data_ratcliff$ID %in% ids, ]


# Data generated under the DMC model
data_dmc <- simulate_data(
  dmc_model, #  the model
  n = 200, # choose this value reasonably
  k = 10, # in practice, simulate more data sets!
  lower = lower_dmc, # the lower parameter space for simulation
  upper = upper_dmc # the upper parameter space for simulation
)
data_dmc <- data_dmc$synth_data # extract the data.frame of raw RTs

## -----------------------------------------------------------------------------
calc_stats(data_ratcliff, type = "basic_stats", level = "group")
calc_stats(data_dmc, type = "basic_stats", level = "group")

## ----eval = FALSE-------------------------------------------------------------
# # get the parameter ranges for optimization
# # (adapt this if you have a custom model)
# l_u_ratcliff <- get_lower_upper(ratcliff_model)
# l_u_dmc <- get_lower_upper(dmc_model)
# 
# 
# # Fit the Ratcliff model to its data
# fits_r_r <- estimate_dm(
#   drift_dm_obj = ratcliff_model,
#   obs_data = data_ratcliff,
#   optimizer = "DEoptim",
#   lower = l_u_ratcliff$lower,
#   upper = l_u_ratcliff$upper,
#   n_cores = 1, # you might want to increase this
#   verbose = 0,
#   progress = 0,
#   messaging = 0 # suppresses all status infos
# )
# 
# # Fit the Ratcliff model to DMC's data
# fits_r_d <- estimate_dm(
#   drift_dm_obj = ratcliff_model,
#   obs_data = data_dmc,
#   optimizer = "DEoptim",
#   lower = l_u_ratcliff$lower,
#   upper = l_u_ratcliff$upper,
#   n_cores = 1, # you might want to increase this
#   verbose = 0,
#   progress = 0,
#   messaging = 0 # suppresses all status infos
# )
# 
# # Fit the DMC model to its data
# fits_d_d <- estimate_dm(
#   drift_dm_obj = dmc_model,
#   obs_data = data_dmc,
#   optimizer = "DEoptim",
#   lower = l_u_dmc$lower,
#   upper = l_u_dmc$upper,
#   n_cores = 1, # you might want to increase this
#   verbose = 0,
#   progress = 0,
#   messaging = 0 # suppresses all status infos
# )
# 
# # Fit the DMC model to the Ratcliff model's data
# fits_d_r <- estimate_dm(
#   drift_dm_obj = dmc_model,
#   obs_data = data_ratcliff,
#   optimizer = "DEoptim",
#   lower = l_u_dmc$lower,
#   upper = l_u_dmc$upper,
#   n_cores = 1, # you might want to increase this
#   verbose = 0,
#   progress = 0,
#   messaging = 0 # suppresses all status infos
# )

## ----eval = FALSE-------------------------------------------------------------
# # Ratcliff model -> Ratcliff data
# stats_r_r <- BIC(fits_r_r, type = "fit_stats")
# # DMC model -> Ratcliff data
# stats_d_r <- BIC(fits_d_r, type = "fit_stats")
# # Ratcliff model -> DMC data
# stats_r_d <- BIC(fits_r_d, type = "fit_stats")
# # DMC model -> DMC data
# stats_d_d <- BIC(fits_d_d, type = "fit_stats")
# 
# # ensure the order of fits matches
# stopifnot(stats_r_r$ID == stats_d_r$ID)
# stopifnot(stats_r_d$ID == stats_d_d$ID)
# 
# # now count how often for each data set a model
# winner_ratcliff_data <- ifelse(stats_r_r$BIC < stats_d_r$BIC, 1, 2)
# winner_ratcliff_data <- factor(winner_ratcliff_data, levels = c(1, 2))
# winner_dmc_data <- ifelse(stats_r_d$BIC < stats_d_d$BIC, 1, 2)
# winner_dmc_data <- factor(winner_dmc_data, levels = c(1, 2))

## ----eval = FALSE-------------------------------------------------------------
# # convert to proportions
# table_ratcliff <- prop.table(table(winner_ratcliff_data))
# table_dmc <- prop.table(table(winner_dmc_data))
# 
# # assemble matrix
# mat <- rbind(table_ratcliff, table_dmc)
# rownames(mat) <- c("sim_r", "sim_d")
# colnames(mat) <- c("fit_r", "fit_d")

## ----echo = FALSE, eval = FALSE-----------------------------------------------
# use_directory("inst")
# saveRDS(
#   object = mat,
#   file = file.path("inst", "model_recovery_vignette_mat.rds")
# )

## ----echo = FALSE-------------------------------------------------------------
sys_path <- system.file(package = "dRiftDM")
mat <- readRDS(
  file = file.path(sys_path, "model_recovery_vignette_mat.rds")
)

## -----------------------------------------------------------------------------
mat

## -----------------------------------------------------------------------------
inversion_mat <- apply(mat, 2, function(col) col / sum(col))
inversion_mat