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)
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)
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)
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
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)
(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))
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 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.)
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)
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 ) )
(fit_stats <- itchmodel::get_fit_stats(optim_output, model = "DDM", parameterization = parameterization, n_data_points = nrow(stimuli)) )
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.