In bramzandbelt/itchmodel: Cognitive model of time framing effects in intertemporal choice

Overview

The aim of this notebook is to determine whether the optimization procedure can recover model parameters with which data were simulated.

Parameter recovery will be performed for the four core models ("parameterizations") that will be tested: - scaling of the value function - scaling of the time function - sensitivity of the value function (optional) - sensitivity of the time function (optional)

Preliminaries

Clear data, values, and functions

rm(list = ls())

Make magrittr's pipe accessible from notebook.

"%>%" <- magrittr::`%>%`

Set state for random number generation

set.seed(19821101)

Model settings

Parameterization, parameter bounds, and parameter names

parameterization <- "defer_speedup_value_scaling"

# Parameter names
names <- get_par_names(model = "DDM", 
                       parameterization = parameterization)

# Lower and upper bounds
lowers <- get_par_bounds(model = "DDM", 
                         parameterization = parameterization, 
                         bound = "lower")
uppers <- get_par_bounds(model = "DDM", 
                         parameterization = parameterization, 
                         bound = "upper")

# Show
tibble::tibble(names = names, 
               lower_bound = lowers, 
               upper_bound = uppers)

Number of free parameters

n_free_param <- sum(!(lowers == uppers))

Default values, group-average parameters from Dai & Busemeyer, JEP General, 2014, Table 10

# def_alpha <- 0.85 # Sensitivity of the value function
# def_mu <- 1 # Scaling of the value function
# def_beta <- 0.70 # Sensitivity of the time function
# def_kappa <- 1 # Scaling of the time function
# def_w <- 0.5 # Attentional weight
# def_a <- 0.5 # Threshold separation
# def_t0 <- 0.3 # Non-decision time
# 
def_parameter_list <- 
  list('alpha' = .94, # Sensitivity of the value function, varies between frames
       'mu' = 1, # Scaling of the value function
       'beta' = .6, # Sensitivity of the time function
       'kappa' = 1, # Scaling of the time function
       'w' = 0.48, # Attentional weight
       "a" = 1.94/2, # Threshold separation
       "t0" = 1.23 # Non-decision time,
       )

# Parameters
parameters <- def_parameter_list

if (parameterization == "defer_speedup_value_scaling") {
  # Scaling of the value function varies between frames: mu_neutral < mu_defer = mu_speedup (see Scholten & Read, JEPHPP, 2013 for more details)
  parameters[['mu']] <- c(1, 1.1, 1.1)
} else if (parameterization == "defer_speedup_time_scaling") {
  # Scaling of the time function varies between frames: kappa_defer > kappa_neutral > kappa_speedup (see Scholten & Read, JEPHPP, 2013 for more details)
  kappa_coef <- 1.1

  parameters[['kappa']] <- c(1, kappa_coef, 1/kappa_coef)
} else if (parameterization == "date_delay_value_scaling") {
  # Sensitivity of the value function varies between frames: mu_delay < mu_date
  parameters[['mu']] <- c(1, 1.1)
} else if (parameterization == "date_delay_time_scaling") {
  # Scaling of the time function varies between frames: kappa_delay > kappa_date
  parameters[['kappa']] <- c(1.5, 1)
}

sampled_par_vals <- unlist(parameters)

Task settings

Define task parameters

m_ll <- 43.52 # 0.17 * 2^8 = 43.52 -> this means we can vary amounts in steps of 17 or 34 cents
t_ss <- 0 # small amount delivered immediately
t_lls <- c(2,4,8,16,32,64,128) # Exponentially increasing delays, recommended by Frye et al., JOVE, 2016
n_staircase_trials <- 7 # For each delay, participants perform 8 choices in staircase procedure
n_reps <- 5 # Number of repetitions of each trial

Simulate indifference point procedure

1-up-1-down procedure - small number of decreasing steps

(responses_1U1D_decreasing_step <- 
   itchmodel::staircase_procedure(m_ll = m_ll, 
                                  t_ss = t_ss, 
                                  t_lls = t_lls,
                                  n_rep = 1, 
                                  procedure = '1up1down',
                                  parameters = def_parameter_list, 
                                  step_size_decreasing = TRUE,
                                  decrease_every_trial = TRUE,
                                  n_trial = n_staircase_trials
                                  )
 )

itchmodel::plot_my_indifference_point_procedure(responses_1U1D_decreasing_step)

Determine data points to be used in experiment

(experimental_stimuli <- 
  responses_1U1D_decreasing_step %>%
  dplyr::filter(trial == max(trial)) %>%
  dplyr::select(t_ll, m_ss, true_i_p) %>%
  dplyr::rename(around = m_ss) %>%
  dplyr::mutate(below = around - m_ll * 2^-(n_staircase_trials-1),
                above = around + m_ll * 2^-(n_staircase_trials-1)
                ) %>%
  tidyr::gather(key = "m_ss_type", value = "m_ss", below, around, above) %>%
  dplyr::arrange(t_ll, m_ss))

Make task stimuli

stimuli <- 
  experimental_stimuli %>%
  dplyr::rename(t_l = t_ll,
                m_s = m_ss) %>%
  tidyr::crossing(frame = factor(itchmodel::get_frames("defer_speedup_value_scaling"), 
                               levels = itchmodel::get_frames("defer_speedup_value_scaling")),
                  rep = 1:n_reps
                  ) %>%
  dplyr::mutate(m_l = m_ll,
                t_s = t_ss) %>%
  dplyr::select(frame, m_s, t_s, m_l, t_l, m_ss_type)

Optimization settings

Optimization parameters

max_iter <- 500
reltol <- 1e-6
np_per_free_param <- 20

(Additional constraints, other than bounds, may be set inside the cost function, ll_diffusion.)

Simulate synthetic data

This needs to be implemented later:

Sample parameter values from a uniform distribution with bounds

# sampled_par_vals <- 
#   purrr::pmap_dbl(.l = list(min = lowers, 
#                             max = uppers),
#                   .f = runif, # Sample from uniform distribution
#                   n = 1)
# 
# names(sampled_par_vals) <- 
#   get_par_names(model = "DDM",
#                 parameterization = parameterization)
# 
# tibble::tibble(names = names, 
#                values = sampled_par_vals)

# # Put parameters in a format that sim_data can deal with
# sampled_par_vals_per_frame <- 
#   sampled_par_vals %>%
#   itchmodel::get_par_values(model = "DDM", 
#                             parameterization = parameterization)

# parameters <- 
#     list('alpha' = purrr::map_dbl(.x = sampled_par_vals_per_frame, .f = function(x) {x['alpha']}),
#          'mu' = purrr::map_dbl(.x = sampled_par_vals_per_frame, .f = function(x) {x['mu']}),
#        'beta' = purrr::map_dbl(.x = sampled_par_vals_per_frame, .f = function(x) {x['beta']}),
#        'kappa' = purrr::map_dbl(.x = sampled_par_vals_per_frame, .f = function(x) {x['kappa']}),
#        'w' = purrr::map_dbl(.x = sampled_par_vals_per_frame, .f = function(x) {x['w']}),
#        'a' = purrr::map_dbl(.x = sampled_par_vals_per_frame, .f = function(x) {x['a']}),
#        't0' = purrr::map_dbl(.x = sampled_par_vals_per_frame, .f = function(x) {x['t0']})
# )

synthetic_data <- 
  itchmodel::sim_data(stimuli = stimuli,
                      parameters = parameters,
                      parameterization = parameterization) %>% 
  dplyr::select(-parameters)

Fit the model to the data by optimizig parameters

optim_output <- 
  itchmodel::fit_model(data = synthetic_data,
                       model = "DDM", 
                       parameterization = parameterization, 
                       lowers = itchmodel::get_par_bounds(model = "DDM",
                                                          parameterization = parameterization,
                                                          bound = 'lower'),
                       uppers = itchmodel::get_par_bounds(model = "DDM",
                                                          parameterization = parameterization,
                                                           bound = 'upper'),
                       control = list(itermax = max_iter,
                                      reltol = reltol,
                                      NP = np_per_free_param * n_free_param
                                      )
                       )

Evaluate output

Fitting stats

(fit_stats <- 
   itchmodel::get_fit_stats(optim_output, 
                            model = "DDM", 
                            parameterization = parameterization,
                            n_data_points = nrow(stimuli))
 )

Best fitting parameters

The following table shows for each parameter, the best-fitting parameter values (best), the parameter values with which the data were simulated (generating), the lower and upper bounds on the parameters (l_bound and u_bound, respectively), as well as the difference between the best-fitting parameter values and the parameter values with which the data were simulated, and the two bounds (d_generating, d_lower, and d_upper, respectively)

tibble::tibble(names = names, 
               l_bound = lowers,
               generating = sampled_par_vals,
               best = optim_output$optim$bestmem,
               u_bound = uppers,
               d_generating = optim_output$optim$bestmem - sampled_par_vals,
               d_lower = optim_output$optim$bestmem - lowers,
               d_upper = uppers - optim_output$optim$bestmem
               )

Plots

itchmodel::plot_parameter_recovery_results(model = "DDM",
                                           parameterization = parameterization,
                                           orig_params = sampled_par_vals,
                                           optim_output = optim_output
                                           )

# Plot negative maximum likelihood evolved across iterations
itchmodel::plot_optim_fun_evol(optim_output = optim_output)

bramzandbelt/itchmodel documentation built on May 7, 2019, 8:42 a.m.