Getting started with ADDM.jl

In this tutorial we will introduce some of the core functionality of the toolbox. We will define an aDDM with specific parameters, simulate choice and response time using those parameters and recover the known parameters from the simulated data.

Load package

We begin with loading the packages we'll use in this tutorial.

julia> using ADDM, CSV, DataFrames, DataFramesMeta
julia> using Plots, StatsPlots, Random
julia> Random.seed!(38435)Random.TaskLocalRNG()

Parameter recovery on simulated data

Define model

The first component of the toolbox is a parameter container for the model class. This is a key-value pair mapping of parameter names to parameter values. It is a specific kind of dictionary of class ADDM.aDDM. We can create a model container using the ADDM.define_model function.

julia> MyModel = ADDM.define_model(d = 0.007, σ = 0.03, θ = .6, barrier = 1,
                       decay = 0, nonDecisionTime = 100, bias = 0.0)ADDM.aDDM(Dict{Symbol, Any}(:nonDecisionTime => 100, :σ => 0.03, :d => 0.007, :bias => 0.0, :barrier => 1, :decay => 0, :θ => 0.6, :η => 0.0))

Define stimuli

The second ingredient to simulating data is to define stimuli. For the pre-defined model in the module, stimuli consist of values associated with the options over which a decision is made. The stimuli should eventually be arranged into a NamedTuple with required field names (case sensitive): valueLeft and valueRight. There are several ways of defining the stimuli within this constraints. Below are a few examples.

Option 1: Load stimuli with corresponding fixation data

The toolbox comes with some datasets from published research. It also includes a function, ADDM.load_data_from_csv(), that can read in datasets stored in CSVs and wrangles it into a format expected by other functions.

ADDM.load_data_from_csv() expects columns parcode,trial, rt (optional), choice (optional), item_left, item_right in the CSVs and convert item_left anditem_right to valueLeft and valueRight. It organizes both the behavioral and the fixation data as a dictionary of ADDM.Trial objects indexed by subject.

Here, we are reading in empirical data that comes with the package but we will not be making use of the observed choices and response times. This is specified with the stimsOnly argument. The empirical data is only used to extract value difference information to index the fixation data correctly. The choices and response times will be simulated below based on the parameters we specified above.

julia> data = ADDM.load_data_from_csv(data_path * "stimdata.csv", data_path * "fixations.csv"; stimsOnly = true);

data is organized as a dictionary with keys corresponding to subject id's and values of arrays containing ADDM.Trials for each subject. Since we are only interested in the values associated with the stimuli, we specify the number of trials we want to simulate (nTrials) and extract that many trials' value information from the data as MyStims.

julia> nTrials = 1400;
julia> MyStims = (valueLeft = reduce(vcat, [[i.valueLeft for i in data[j]] for j in keys(data)])[1:nTrials], valueRight = reduce(vcat, [[i.valueRight for i in data[j]] for j in keys(data)])[1:nTrials]);

Option 2: Read in from CSV

Stimuli can also be read in directly from a file. This approach is simpler to define only the stimuli but would require additional steps to define and organize fixation data.

fn = data_path * "rawstims.csv"
tmp = DataFrame(CSV.File(fn, delim=","))
MyStims = (valueLeft = tmp.valueLeft, valueRight = tmp.valueRight)

Option 3: Create random stimuli

Alternatively, one can simulate stimuli with random values. Similar to reading in stimulus values directly, this approach would also require additional steps to defne and organize fixation data.

Random.seed!(38535)
MyStims = (valueLeft = randn(1000), valueRight = randn(1000))

Define fixation data

The last ingredient to simulate data is fixation patterns. These are necessary because the distinguishing feature of the aDDM is its ability to use eyetracking data to account for attentional biases in choice behavior. Importantly, aDDM treats fixation data as exogenous. This means, that the model does not describe how the fixation patterns arise, but only takes them as given to examine their effects on choice behavior.

The toolbox has a specific structure for fixation data organized in the FixationData type. This type organizes empirical fixations to distributions conditional on fixation type (first, second etc.) and value difference.

First, we extract value difference information from the dataset to use in processing the fixations.

julia> vDiffs = sort(unique([x.valueLeft - x.valueRight for x in data["1"]]));

Then we summarize the empricial data from all subjects as distributions from which the model samples from depending on the value difference and the fixation type (1st, 2nd etc.) using ADDM.process_fixations.

julia> MyFixationData = ADDM.process_fixations(data, fixDistType="fixation", valueDiffs = vDiffs);

Simulate data

Finally, we define additional arguments for the second component of the model, the trial simulator. Here, these are the fixation data, time step for simulations and the maximum length a trial is allowed (in ms). These arguments need to be specified as a NamedTuple, and must have at least two elements. Otherwise the function tries to apply iterate to the single element which would likely end with a MethodError. In this example, we specify timeStep and cutoff in addition to the only required argument without a default fixationData to avoid this.

julia> MyArgs = (timeStep = 10.0, cutOff = 20000, fixationData = MyFixationData);

Note that these are positional arguments for code efficiency.

julia> SimData = ADDM.simulate_data(MyModel, MyStims, ADDM.aDDM_simulate_trial, MyArgs);

Recover parameters using a grid search

Now that we have simulated data with known parameters we can use the third component of the module, the likelihood function to invert the model and recover those parameters from the data.

The built-in optimization algorithm function for this is ADDM.grid_search. It computes the sum of negative log likelihood across all trials in the data, for each parameter combination specified in param_grid. The parameter space defined in param_grid is organized as an array of NamedTuples indicating the parameter names and their specific values for each combination. For this toy example, the parameter space consists of 27 parameter combinations.

julia> fn = data_path * "addm_grid.csv";
julia> tmp = DataFrame(CSV.File(fn, delim=","));
julia> param_grid = NamedTuple.(eachrow(tmp))27-element Vector{NamedTuple{(:d, :sigma, :theta), Tuple{Float64, Float64, Float64}}}:
 (d = 0.003, sigma = 0.01, theta = 0.2)
 (d = 0.005, sigma = 0.01, theta = 0.2)
 (d = 0.007, sigma = 0.01, theta = 0.2)
 (d = 0.003, sigma = 0.03, theta = 0.2)
 (d = 0.005, sigma = 0.03, theta = 0.2)
 (d = 0.007, sigma = 0.03, theta = 0.2)
 (d = 0.003, sigma = 0.05, theta = 0.2)
 (d = 0.005, sigma = 0.05, theta = 0.2)
 (d = 0.007, sigma = 0.05, theta = 0.2)
 (d = 0.003, sigma = 0.01, theta = 0.4)
 ⋮
 (d = 0.003, sigma = 0.01, theta = 0.6)
 (d = 0.005, sigma = 0.01, theta = 0.6)
 (d = 0.007, sigma = 0.01, theta = 0.6)
 (d = 0.003, sigma = 0.03, theta = 0.6)
 (d = 0.005, sigma = 0.03, theta = 0.6)
 (d = 0.007, sigma = 0.03, theta = 0.6)
 (d = 0.003, sigma = 0.05, theta = 0.6)
 (d = 0.005, sigma = 0.05, theta = 0.6)
 (d = 0.007, sigma = 0.05, theta = 0.6)

Having defined the parameter space, we can compute the sum of negative log likelihoods (NLL) for each data point and select the parameter combination with the lowest NLL, which is equivalent to maximizing the likelihood. This combination of parameters is called the maximum likelihood estimate, or the MLE. ADDM.grid_search expects the data, the parameter space, the likelihood function and any fixed parameters for the model as its positional arguments. Additionaly, as a keyword argument, one can specify likelihood_args for values that should be passed onto the likelihood function. Here, we specify the timeStep and the stateStep for solving the Fokker-Planck Equation. The function returns an output dictionary with various components depending on additional arguments that are detailed in other tutorials and the API Reference.

julia> output = ADDM.grid_search(SimData, param_grid, ADDM.aDDM_get_trial_likelihood, Dict(:η=>0.0, :barrier=>1, :decay=>0, :nonDecisionTime=>100, :bias=>0.0); likelihood_args = (timeStep = 10.0, stateStep = 0.1), return_grid_nlls = true);
julia> mle = output[:mle];
julia> all_nll_df = output[:grid_nlls];

Below we display the sum of negative log likelihoods for each parameter combination.

julia> sort!(all_nll_df, [:nll])27×4 DataFrame
 Row │ d        sigma    theta    nll
     │ Float64  Float64  Float64  Float64
─────┼──────────────────────────────────────────
   1 │   0.007     0.05      0.6   8300.66
   2 │   0.007     0.05      0.4   8370.46
   3 │   0.005     0.05      0.6   8468.25
   4 │   0.005     0.05      0.4   8537.43
   5 │   0.007     0.05      0.2   8588.48
   6 │   0.005     0.05      0.2   8678.79
   7 │   0.003     0.05      0.6   8930.31
   8 │   0.003     0.05      0.4   8984.07
  ⋮  │    ⋮        ⋮        ⋮           ⋮
  21 │   0.005     0.01      0.2      2.06312e5
  22 │   0.007     0.01      0.4      2.06312e5
  23 │   0.007     0.01      0.6      2.06312e5
  24 │   0.007     0.01      0.2      2.06312e5
  25 │   0.003     0.01      0.4      2.06312e5
  26 │   0.005     0.01      0.6      2.06312e5
  27 │   0.005     0.01      0.4      2.06312e5
                                 12 rows omitted

output_path = '/home/jovyan/work/all_nll_df.csv'
CSV.write(output_path, all_nll_df)

You might have noticed that the grid search did not identify the true parameters (d = 0.007, σ = 0.03, θ = .6) as the ones with the highest likelihood. This highlights the importance of choosing good stepsizes for the temporal and spatial discretization. Let's reduce the spatial step size and see if we can recover the corect parameter combination.

julia> my_likelihood_args = (timeStep = 10.0, stateStep = 0.01)(timeStep = 10.0, stateStep = 0.01)
julia> output = ADDM.grid_search(SimData, param_grid, ADDM.aDDM_get_trial_likelihood, Dict(:η=>0.0, :barrier=>1, :decay=>0, :nonDecisionTime=>100, :bias=>0.0), likelihood_args=my_likelihood_args);
julia> mle = output[:mle];
julia> mleDict{Symbol, Real} with 13 entries:
  :barrier         => 1
  :decay           => 0
  :timeStep        => 10.0
  :nonDecisionTime => 100
  :σ               => 0.03
  :sigma           => 0.03
  :η               => 0.0
  :bias            => 0.0
  :d               => 0.007
  :nll             => 7824.37
  :θ               => 0.6
  :stateStep       => 0.01
  :theta           => 0.6

Success! Reducing the state step size was sufficient to recover the true parameter combination for the simulated dataset.

Ok, but what happened there? That's what we detail in the next tutorial.