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"
@@ -35,7 +36,7 @@ MixedModels = "4"
 MixedModelsSim = "0.2"
 Parameters = "0.12"
 Random = "1"
-SignalAnalysis = "0.4, 0.5,0.6,0.7,0.8,0.9,0.10"
+SequentialSamplingModels = "0.12"
 Statistics = "1"
 StatsModels = "0.6,0.7"
 ToeplitzMatrices = "0.7, 0.8"

docs/literate/HowTo/driftComponent.jl

@@ -0,0 +1,96 @@
+# # Simulate an Evidence Accumulation EEG
+
+# We want to use the SequenceDesign with a stimulus component, a evidence accumulation component (DriftComponent) and a response component to simulate an full evidence accumulation process.
+
+
+
+# ### 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
+# ```@raw html
+# </details >
+# ```
+
+# ## Design
+# Let's generate a single design
+design_single = SingleSubjectDesign(conditions = Dict(:condition => [1]))
+# Now we convert the SingleSubjectDesign to an SequenceDesign with a Sequence of S: stimulus, C: component, R: response, _: gap
+design_seq = SequenceDesign(design_single, "SCR_")
+# On top we can now repeat the Sequence to have multiple trials. For this experiment we use 100 replicas.
+design_rep = RepeatDesign(design_seq, 100)
+
+# ## Create the 3 components we want to use
+# First the stimulus component as simple LinearModelComponent
+p3 = LinearModelComponent(;
+    basis = UnfoldSim.hanning(Int(0.5 * fs)),
+    formula = @formula(0 ~ 1 + condition),
+    β = [0.4, 0],
+)
+# Second a response component
+resp = LinearModelComponent(;
+    basis = UnfoldSim.hanning(Int(0.5 * fs)),
+    formula = @formula(0 ~ 1 + condition),
+    β = [0.6, 0],
+)
+# Third we create our DriftComponent which implies the evidence accumulation. For the DriftComponent we can choose between different Models which are used for the Simulation.
+v = [0.8] # Drift Rate
+A = 0.1 # The maximum start point, an indicator of an an initial bias towards a decision.
+k = 0.4 # A + k = b, where b is the decision threshold.
+t = 0.2 # The duration for a non-decisional processes (encoding and response execution).
+lba_parameter = Dict(:ν => v, :A => A, :k => k, :τ => t)
+drift = DriftComponent(max_length, fs, LBA, lba_parameter)
+# As last step we have to specify the components as a Dict connection the components with the events of the design.
+components = Dict('S' => [p3], 'C' => [drift], 'R' => [resp])
+
+# ## Create the Onset Component for the design
+# For the stimulus and response we use a simple Uniform onset with a distance of 250 and 300. The DriftComponent uses the special DriftOnset in combination with an UniformOnset. The DriftOnset returns the response time after the model reaches the decision threshold.
+# Important is to know that, the onset for each component defines the onset for the next. So in this case: S->C, C->R, R->S.
+seq_onset = SequenceOnset(
+    Dict(
+        'S' => UniformOnset(width = 0, offset = 250),
+        'C' => (DriftOnset(), UniformOnset(width = 0, offset = 150)),
+        'R' => UniformOnset(width = 0, offset = 300),
+    ),
+)
+
+# ## 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 ERP of simulated EEG
+# First the simulated EEG
+lines(data)
+vlines!(evts.latency[evts.event.=='S'], color = (:green, 0.5))
+vlines!(evts.latency[evts.event.=='C'], color = (:blue, 0.5))
+vlines!(evts.latency[evts.event.=='R'], color = (:orange, 0.5))
+CairoMakie.xlims!(0, 2000)
+current_figure()
+# Second the ERP
+evts.event = string.(evts.event)
+data_epochs, times_epoch = Unfold.epoch(data = data, tbl = evts, τ = (0, 1.0), sfreq = fs);
+f = @formula(0 ~ 1)
+m = fit(
+    UnfoldModel,
+    ["S" => (f, times_epoch), "R" => (f, times_epoch), "C" => (f, times_epoch)],
+    evts,
+    data_epochs,
+);
+results = coeftable(m)
+plot_erp(results; mapping = (; color = :eventname))

docs/literate/HowTo/simulateDriftOverlap.jl

@@ -0,0 +1,104 @@
+# # 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
+
+p3 = LinearModelComponent(;
+    basis = UnfoldSim.hanning(Int(0.7 * fs)),
+    formula = @formula(0 ~ 1 + condition),
+    β = [0.8, 0],
+)
+resp = LinearModelComponent(;
+    basis = UnfoldSim.hanning(Int(0.7 * fs)),
+    formula = @formula(0 ~ 1 + condition),
+    β = [0.5, 0],
+)
+
+# ```@raw html
+# </details >
+# ```
+# ## Design
+# Let's generate a single design with two different drift_rates as condition, so we can use them later for the Simulation
+design_single =
+    SingleSubjectDesign(conditions = Dict(:drift_rate => [3, 5], :condition => [1]))
+design_seq = SequenceDesign(design_single, "SCR_")
+design_rep = RepeatDesign(design_seq, 500)
+
+# ## Create the Drift Component using the special KellyModel for the evidence accumulation
+# As parameter for the models drift_rate we use the drift_rate specified in the design by directly naming the condition as string.
+v = "drift_rate" # Drift Rate from the design
+# Now we retrieve the models default parameters, but we change the drift_rate and
+kelly_model = KellyModel(drift_rate = v, motor_onset = 0.4, event_onset = 0.2)
+kelly_model_parameter = UnfoldSim.create_kelly_parameters_dict(kelly_model)
+drift = DriftComponent(max_length, fs, KellyModel, kelly_model_parameter)
+components = Dict('S' => [p3], 'C' => [drift], 'R' => [resp])
+# ## Create the Onset Component to simulate the overlap
+# To simulate an overlap between the drift_component in 'C' and the response component in 'R'. We have to specify the UniformOnset of the 'C' Component therefore with an negative offset to produce the overlap.
+seq_onset = SequenceOnset(
+    Dict(
+        'S' => UniformOnset(width = 0, offset = 1.2 * fs),
+        'C' => (DriftOnset(), UniformOnset(width = 0, offset = -100)),
+        'R' => UniformOnset(width = 0, offset = 2 * fs),
+    ),
+)
+
+# ## 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
+lines(data)
+vlines!(evts.latency[evts.event.=='S'], color = (:green, 0.5))
+vlines!(evts.latency[evts.event.=='C'], color = (:blue, 0.5))
+vlines!(evts.latency[evts.event.=='R'], color = (:orange, 0.5))
+CairoMakie.xlims!(0, 2000)
+current_figure()
+
+# ## Plot ERP with the overlap in the 'R' Component
+evts.event = string.(evts.event)
+data_epochs, times_epoch =
+    Unfold.epoch(data = data, tbl = evts, τ = (-0.1, 0.9), sfreq = fs);
+f = @formula(0 ~ 1)
+m = fit(
+    UnfoldModel,
+    ["S" => (f, times_epoch), "R" => (f, times_epoch), "C" => (f, times_epoch)],
+    evts,
+    data_epochs,
+);
+results = coeftable(m)
+plot_erp(results; mapping = (; color = :eventname))
+
+# ## Deconvolution of the overlap
+fir = firbasis(τ = (-0.1, 0.9), sfreq = fs)
+evts2 = deepcopy(evts)
+evts2.event = string.(evts2.event)
+m = fit(
+    UnfoldModel,
+    ["S" => (f, fir), "R" => (f, deepcopy(fir)), "C" => (f, deepcopy(fir))],
+    evts2,
+    data,
+);
+results = coeftable(m)
+plot_erp(results; mapping = (; color = :eventname))

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("sequentialSamplingModelSimulation.jl")
 
 export size, length
 export AbstractComponent, AbstractNoise, AbstractOnset, AbstractDesign
@@ -78,4 +81,7 @@ export AbstractHeadmodel, Hartmut, headmodel, leadfield, orientation, magnitude
 
 # multichannel
 export MultichannelComponent
+
+# evidence accumulation
+export DriftComponent, DriftOnset, SequenceOnset
 end

src/component.jl

@@ -706,3 +706,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]