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

docs/literate/HowTo/driftComponent.jl

@@ -0,0 +1,100 @@
+# # 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
+# ```@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 Drift_Component which implies the evidence accumulation. For the Drift_Component 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 = Drift_Component(
+    simulate_component,
+    time_vec,
+    Δt,
+    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 Drift_Component 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,97 @@
+# # 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
+
+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 = Drift_Component(
+    simulate_component,
+    time_vec,
+    Δt,
+    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

@@ -18,6 +18,8 @@ using Automa # for sequence
 
 using LinearAlgebra # headmodel
 
+using SequentialSamplingModels # for SequentialSamplingModels
+
 import DSP.hanning
 import Base.length
 import Base.size
@@ -33,6 +35,7 @@ include("headmodel.jl")
 include("helper.jl")
 include("sequence.jl")
 include("bases.jl")
+include("sequentialSamplingModelSimulation.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 Drift_Component, DriftOnset, SequenceOnset, KellyModel
 end

src/component.jl

@@ -513,3 +513,133 @@ function simulate_and_add!(
     @views epoch_data[:, 1+off:length(component)+off, :] .+=
         simulate_component(rng, component, simulation)
 end
+
+
+
+"""
+    Drift_Component <: AbstractComponent
+A component that adds an evidence accumulation signal according to a sequential sampling model selected in the field model_type.
+
+All fields are mandatory. Works best with [`SequenceDesign`](@ref).
+
+# Fields
+- `basis`: an object, if accessed, provides a 'basis-function', e.g. `hanning(40)`, this defines the response at a single event. It will be weighted by the model-prediction. Can also be a tuple `(fun::Function,maxlength::Int)` with a function `fun` that either generates a matrix `size = (maxlength,size(design,1))` or a vector of vectors. If a larger matrix is generated, it is automatically cutoff at `maxlength`
+- `time_vec`: Vector which defines the length of the signal
+- `Δt`: Float that indicates the timesteps used to generate the time_vec
+- `model_type`: Model struct which defines the model to use to generate the traces, e.g. `KellyModel`
+- `model_parameters`: Dict. Containing the parameters for the 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
+fs=500
+Δt = 1/fs;
+tEnd = 1.0
+time_vec = 0:Δt:tEnd;
+model_parameter = create_kelly_parameters_dict(KellyModel())
+Drift_Component(
+    simulate_component,
+    time_vec,
+    Δt,
+    KellyModel,
+    model_parameter)
+```
+
+See also [`create_kelly_parameters_dict`](@ref), [`KellyModel`](@ref).
+"""
+struct Drift_Component <: AbstractComponent
+    basis::Any
+    time_vec::Any
+    Δt::Any
+    model_type::Any
+    model_parameters::Dict
+end
+"""
+    simulate_component(rng, c::Drift_Component, design::AbstractDesign)
+
+Generate evidence accumulation traces by using the model and its parameters specified in the component c.
+
+# Returns
+- `Matrix{Float64}`: Simulated component for each event in the events data frame. The output dimensions are `length(c) 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 = create_kelly_parameters_dict(KellyModel());
+
+julia> c = Drift_Component(simulate_component, 0:1/500:1.0, 1/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::Drift_Component, design::AbstractDesign)
+    _, traces = trace_sequential_sampling_model(deepcopy(rng), c, design)
+    return traces
+end
+"""
+    calculate_response_times_for_ssm(rng, component::Drift_Component, design::AbstractDesign)
+
+Generate response times of the evidence accumulation by using the model and its parameters specified in the component.
+
+Using same rng as in simulate_component(rng, component::Drift_Component, design::AbstractDesign) to ensure that the response times match the generated traces. 
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure the same traces are created as in the use of the [`simulate_component`](@ref) function.
+- `component::Drift_Component`: 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 time 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 = create_kelly_parameters_dict(KellyModel());
+
+julia> c = Drift_Component(simulate_component, 0:1/500:1.0, 1/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::Drift_Component, design::UnfoldSim.SubselectDesign)
+    rts, _ = trace_sequential_sampling_model(deepcopy(rng), component, design)
+    return rts
+end
+Base.length(c::Drift_Component) = length(c.time_vec)
+"""
+    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]