Feature sequential sampling models

Project.toml

@@ -17,6 +17,7 @@ MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 MixedModelsSim = "d5ae56c5-23ca-4a1f-b505-9fc4796fc1fe"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SequentialSamplingModels = "0e71a2a6-2b30-4447-8742-d083a85e82d1"
 SignalAnalysis = "df1fea92-c066-49dd-8b36-eace3378ea47"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
@@ -36,6 +37,7 @@ MixedModels = "4,5"
 MixedModelsSim = "0.2"
 Parameters = "0.12"
 Random = "1"
+SequentialSamplingModels = "0.12"
 SignalAnalysis = "0.4, 0.5,0.6,0.7,0.8,0.9,0.10"
 Statistics = "1"
 StatsModels = "0.6,0.7"

docs/Project.toml

@@ -10,6 +10,7 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SequentialSamplingModels = "0e71a2a6-2b30-4447-8742-d083a85e82d1"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"

docs/literate/tutorials/simulateDriftOverlap.jl

@@ -0,0 +1,152 @@
+# # Simulate an Evidence Accumulation Overlap and Deconvolution
+
+# We want to use the SequenceDesign with three components to simulate an full evidence accumulation process with a Overlap between the evidence accumulation and the response signal. Afterwards we use the Deconvolution Method from Unfold to unfold the signals again.
+
+
+
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+using UnfoldSim
+using Unfold
+using StableRNGs
+using CairoMakie, UnfoldMakie
+using SequentialSamplingModels
+
+fs = 500
+Δt = 1 / fs; # time step
+tEnd = 1.0 # trial Duration
+time_vec = 0:Δt:tEnd; # time base - let's make it a typical stimulus duration
+max_length = tEnd / Δt
+
+c_stimulus = LinearModelComponent(;
+    basis = ones(Int(0.1 * fs)),
+    formula = @formula(0 ~ 1 + condition),
+    β = [1, 0],
+)
+c_response = LinearModelComponent(;
+    basis = ones(Int(0.2 * fs)),
+    formula = @formula(0 ~ 1 + condition),
+    β = [-0.5, 0],
+)
+
+# ```@raw html
+# </details >
+# ```
+# ## Design
+# Let's generate a single design with two drift_rates
+design_single =
+    SingleSubjectDesign(conditions = Dict(:drift_rate => [1, 2], :condition => [1]))
+design_seq = SequenceDesign(design_single, "SCR_")
+design_rep = RepeatDesign(design_seq, 100) # 100 SCR-repeats
+
+# # Simplified S-C-R drift-model with LBA
+# ## LBA
+# The following parameters are described in the [SequentialSamplingModels.LBA model documentation](https://itsdfish.github.io/SequentialSamplingModels.jl/dev/lba/)
+v = [0.8] ## Drift Rate, needs to be in [ ] for compatability with `SequentialSamplingModels.jl``
+A = 0.2 ## The starting point is sampled uniformly between `[0, A]`
+k = 0.4 ## Difference of A with the bound/threshold ($A + k = b$), where b is the decision threshold.
+
+lba_parameter = Dict(:ν => v, :A => A, :k => k, :τ => 0.0) # we didnt specify σ thus it is 1 by default,
+
+# !!! important
+# The non-decision time has to be specified here as "τ=>0." because it is not defined whether non-decision time applies to stimulus encoding or response execution. If you have a non-0 value here, the events following the DriftComponent will be shifted by τ.
+drift = DriftComponent(max_length, fs, LBA, lba_parameter)
+
+# Let's have a quick look at the simulation
+series(
+    hcat(simulate_component.(StableRNG.(1:10), Ref(drift), Ref(design_single))...)',
+    solid_color = (:black, 0.5),
+)
+# Super! Our start values are sampled between [0, 0.2] as indicated, and the boundary is at A+k = 0.2 + 0.4 = 0.6
+
+# ### Simulate EEG data 
+components = Dict('S' => [c_stimulus], 'C' => [drift], 'R' => [c_response])
+
+# Next we define a SequenceOnset ingredient, this allows us fine-grained control on when events should happen.
+
+seq_onset = SequenceOnset(
+    Dict(
+        'S' => UniformOnset(width = 0, offset = 0.2 * fs), # distance S<->C - we could choose 0.1*fs to immediately start after the stimulus component
+        'C' => DriftOnset(), ## automatically determines which is the respective DriftComponent
+        'R' => UniformOnset(width = 0, offset = 2 * fs), # distance R<->S
+    ),
+)
+
+# ## Simulate data with the created design, components and onsets
+data, evts = UnfoldSim.simulate(
+    StableRNG(12),
+    design_rep,
+    components,
+    seq_onset,
+    NoNoise(), # PinkNoise(noiselevel=1)
+)
+# ## Plot EEG to see the overlap as it appears
+f_data, _ax, _h = lines(data[1:5000])
+vlines!(evts.latency, color = Int.(evts.event), alpha = 0.3)
+xlims!(700, 3500)
+f_data
+
+# # Complex  S-C-R drift-model ala `Kelly et al.`
+
+# In addition to changing the model to a more realistic version, we will also make use of the `drift_rate` specification in our `design``
+
+v = "drift_rate" ## take the drift rate directly from the design
+
+# The KellyModel has _many_ default parameters, let's grab them and immediately change a few
+
+kelly_model =
+    UnfoldSim.KellyModel(drift_rate = v, motor_delay = 0.4, sensor_encoding_delay = 0.2)
+
+kelly_model_parameter = UnfoldSim.create_kelly_parameters_dict(kelly_model)
+drift = DriftComponent(max_length, fs, UnfoldSim.KellyModel, kelly_model_parameter)
+components = Dict('S' => [c_stimulus], 'C' => [drift], 'R' => [c_response])
+
+# We are ready to simulate!
+data, evts = UnfoldSim.simulate(
+    StableRNG(12),
+    design_rep,
+    components,
+    seq_onset,
+    NoNoise(), # PinkNoise(noiselevel=1)
+)
+# ## Plot EEG to see the overlap as it appears
+f_data, _ax, _h = lines(data[1:5000])
+vlines!(evts.latency, color = Int.(evts.event), alpha = 0.3)
+xlims!(700, 3500)
+f_data
+
+# !!! note
+# There are varying gaps between the stimulus offset <-> driftstart and driftend <-> response start. This is because the `KellyModel` by default has a `sensor_encoding_delay` and a  `motor_delay` (including a `motor_delay_variability` etc.)
+
+# # ERP & Deconvolution 
+# Let's calculate traditional ERPs and deconvolved ones!
+evts.event = string.(evts.event) # due to a Unfold-bug we have to convert :S
+data_epochs, times_epoch =
+    Unfold.epoch(data = data, tbl = evts, τ = (-1.3, 1.3), sfreq = fs);
+f = @formula(0 ~ 1 + drift_rate)
+m = fit(UnfoldModel, ["S" => (f, times_epoch), "R" => (f, times_epoch)], evts, data_epochs);
+plot_erp(
+    effects(Dict(:drift_rate => [2, 3]), m);
+    mapping = (;
+        color = :drift_rate => nonnumeric,
+        col = :eventname => UnfoldMakie.sorter(["S", "R"]),
+    ),
+)
+
+# ## Deconvolution of the overlap
+fir_stim = firbasis(τ = (-0.1, 1.3), sfreq = fs)
+fir_resp = firbasis(τ = (-1, 0.3), sfreq = fs)
+
+m = fit(UnfoldModel, ["S" => (f, fir_stim), "R" => (f, fir_resp)], evts, data);
+plot_erp(
+    effects(Dict(:drift_rate => [2, 3]), m);
+    mapping = (;
+        color = :drift_rate => nonnumeric,
+        col = :eventname => UnfoldMakie.sorter(["S", "R"]),
+    ),
+)
+# Deconvolution sucessfully removed the stimulus and response overlaps from our estimates, but it kept the drift-diffusion aspect nearly untouched. It is important to stress the "nearly", because depending on the `delay` parameters, this can potentially lead to problematic situation where it is unclear whether the drift activity should be assigned to S or R. In most realistic situations we found to be safe.
+# Further note that the drift-diffusion amplitude can be reduced because it is split up between `S` and `R` events.`

docs/make.jl

@@ -41,6 +41,7 @@ makedocs(;
             "Simulate event-related potentials (ERPs)" => "generated/tutorials/simulateERP.md",
             "Power analysis" => "generated/tutorials/poweranalysis.md",
             "Multi-subject simulation" => "generated/tutorials/multisubject.md",
+            "S-C-R & Evidence accumulation models" => "./generated/tutorials/simulateDriftOverlap.md",
         ],
         "Reference" => [
             "Overview of functionality" => "./generated/reference/overview.md",
@@ -58,6 +59,7 @@ makedocs(;
             "Use existing experimental designs & onsets in the simulation" => "./generated/HowTo/predefinedData.md",
             "Simulate ground truth marginal effects" => "./generated/HowTo/getGroundTruth.md",
             "Sequence of events (e.g. SCR)" => "./generated/HowTo/sequence.md",
+            #"Simulate LBA evidence accumulation models" => "./generated/HowTo/driftComponent.md",
         ],
         "Developer documentation" => "developer_docs.md",
         "API / Docstrings" => "api.md",

src/UnfoldSim.jl

@@ -17,6 +17,8 @@ using HDF5, Artifacts, FileIO
 using Automa # for sequence
 using LinearAlgebra # headmodel
 
+using SequentialSamplingModels # for SequentialSamplingModels
+
 import DSP.hanning
 import Base.length
 import Base.size
@@ -32,6 +34,7 @@ include("headmodel.jl")
 include("helper.jl")
 include("sequence.jl")
 include("bases.jl")
+include("sequentialSamplingModels.jl")
 
 export size, length
 export AbstractComponent, AbstractNoise, AbstractOnset, AbstractDesign
@@ -79,4 +82,7 @@ export AbstractHeadmodel, Hartmut, headmodel, leadfield, orientation, magnitude
 
 # multichannel
 export MultichannelComponent
+
+# evidence accumulation
+export DriftComponent, DriftOnset, SequenceOnset, KellyModel, LBA, DDM
 end

src/component.jl

@@ -737,3 +737,128 @@ function simulate_and_add!(
     @views epoch_data[:, (1+off):(length(component)+off), :] .+=
         simulate_component(rng, component, simulation)
 end
+
+
+
+"""
+    DriftComponent <: AbstractComponent
+A component that adds an evidence accumulation trace, following a specified sequential sampling model.
+
+Works best with [`SequenceDesign`](@ref).
+
+# Fields
+- `max_length::Int`: maximum length of the simulated trace.
+- `sfreq::Real`: sample frequency used to simulate the trace, used to evaluate the parameters in the model_type.
+- `model_type`: Model struct which defines the model to use to generate the trace, e.g. `KellyModel`
+- `model_parameters`: Dict. Containing the parameters for the sequential sampling simulation model specified in `model_type`.
+
+# Examples
+```julia-repl
+# use the KellyModel and its default parameters to simulate traces from 0:1/500:1.0
+sfreq=500
+tEnd = 1.0
+max_length = tEnd*sfreq;
+model_parameter = Dict()
+DriftComponent(
+    max_length,
+    sfreq,
+    KellyModel,
+    model_parameter)
+```
+
+"""
+struct DriftComponent <: AbstractComponent
+    max_length::Int
+    sfreq::Real
+    model_type::Any
+    model_parameters::Dict
+end
+"""
+    simulate_component(rng, c::DriftComponent, design::AbstractDesign)
+
+Generate evidence accumulation traces based on the model specified in `c::DriftComponent` and (optionally) parameters from `design::AbtractDesign`
+
+# Returns
+- `Matrix{Float64}`: Simulated component for each event in the events data frame. The output dimensions are `length(c.time_vec) x size(events, 1)`.
+
+# Examples
+```julia-repl
+# use the KellyModel and its default parameters to simulate traces from 0:1/500:1.0
+julia> model_parameter = Dict();
+
+julia> c = DriftComponent(500, 500, KellyModel, model_parameter);
+
+julia> design_single = SingleSubjectDesign(conditions = Dict(:drift_rate => [0.5, 0.8], :condition => [1]));
+
+julia> design_seq = SequenceDesign(design_single,"SCR_");
+
+julia> simulate_component(StableRNG(1),c,design_seq)
+501x6 Matrix{Float64}:
+0.0  0.0  0.0  0.0  0.0  0.0
+0.0  0.0  0.0  0.0  0.0  0.0
+⋮                        ⋮
+```
+"""
+function UnfoldSim.simulate_component(rng, c::DriftComponent, design::AbstractDesign)
+    _, traces = simulate_drift_component(deepcopy(rng), c, design)
+    return traces
+end
+"""
+    calculate_response_times_for_ssm(rng, component::DriftComponent, design::AbstractDesign)
+
+Generate response times of the evidence accumulation by using the model and its parameters specified in the component.
+
+Important: specify the same rng as in `simulate_component(rng, component::DriftComponent, design::AbstractDesign)` to ensure that the response times match the generated traces. 
+
+# Arguments
+- `rng::AbstractRNG`: Random seed to ensure the same traces are created as in the use of the [`simulate_component`](@ref) function.
+- `component::DriftComponent`: Component to specify the model and its parameters to simulate the evidence accumulation.
+- `design::UnfoldSim.SubselectDesign`: Subselection of the Sequence Design to ensure we only generate rt for the drift component events.
+
+# Returns
+- `Vector{Float64}`: Simulated response times for each event in the events data frame. The output dimension is `size(events, 1)`.
+
+```julia-repl
+# use the KellyModel and its default parameters to simulate traces from 0:1/500:1.0
+julia> model_parameter = Dict();
+
+julia> c = DriftComponent(500, 500, KellyModel, model_parameter);
+
+julia> design_single = SingleSubjectDesign(conditions = Dict(:drift_rate => [0.5, 0.8], :condition => [1]));
+
+julia> design_seq = SequenceDesign(design_single,"SCR_");
+
+julia> calculate_response_times_for_ssm(StableRNG(1),c,design_seq)
+2-element Vector{Float64}:
+ 260.70134768436486
+ 360.1329203034039
+```
+"""
+function calculate_response_times_for_ssm(rng, component::DriftComponent, design::UnfoldSim.SubselectDesign)
+    rts, _ = simulate_drift_component(deepcopy(rng), component, design)
+    return rts
+end
+Base.length(c::DriftComponent) = c.max_length
+"""
+    get_model_parameter(rng, evt, d::Dict)
+
+Construct a parameter dictionary of parameter names and values.
+"""
+function get_model_parameter(rng, evt, d::Dict)
+    result_parameter = Dict()
+    for key in keys(d)
+        result_parameter[key] = get_model_parameter(rng, evt, d[key])
+    end
+    return result_parameter
+end
+"""
+    get_model_parameter(rng, evt, val)
+
+Recursive default call return the model parameter as it is.
+"""
+get_model_parameter(rng, evt, val) = val
+"""
+    get_model_parameter(rng, evt, val::String)
+Recursive call if the model parameter specified as string from the conditions in the events data frame.
+"""
+get_model_parameter(rng, evt, val::String) = evt[val]