OnsetFormulas, Sequences, and Component Functions - v0.4.0 -> actually v0.5.0

Project.toml

@@ -5,6 +5,7 @@ version = "0.4.1"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+Automa = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

docs/Project.toml

@@ -15,7 +15,5 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679"
 UnfoldMakie = "69a5ce3b-64fb-4f22-ae69-36dd4416af2a"
+UnfoldMixedModels = "019ae9e0-8363-565c-86e5-97a5a2fe84f4"
 UnfoldSim = "ed8ae6d2-84d3-44c6-ab46-0baf21700804"
-
-[compat]
-Unfold = "0.7"

docs/literate/HowTo/componentfunction.jl

@@ -0,0 +1,69 @@
+# # Component Basisfunctions
+# HowTo use functions that depend on the `design` and return  per-event basis-vectors, instead of the same basis vector for all events.
+
+
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+## Load required packages
+using UnfoldSim
+using Unfold
+using Random
+using DSP
+using CairoMakie, UnfoldMakie
+# ```@raw html
+# </details >
+# ```
+
+
+sfreq = 100;
+
+# ## Design
+# Let's generate a design with a categorical effect and a continuous duration effect
+design = UnfoldSim.SingleSubjectDesign(;
+    conditions = Dict(
+        :category => ["dog", "cat"],
+        :duration => Int.(round.(20 .+ rand(100) .* sfreq)),
+    ),
+);
+
+
+# Instead of defining a "boring" vector basis function e.g. `[0,0,1,2,3,3,2,1,0,0,0]`, let's use function - in our case a hanning window with the size depending on the experimental design's duration.
+# !!! important
+#     Two things have to be taken care of:
+#     1. in case a rng is required to e.g. generate the design, or your basisfunction depends on it, you have to specify a two-argument basis-function: `(rng,design)->...`
+#     2. a `maxlength` has to be specified via a tuple `(function,maxlength)``
+
+mybasisfun = design -> hanning.(generate_events(design).duration)
+signal = LinearModelComponent(;
+    basis = (mybasisfun, 100),
+    formula = @formula(0 ~ 1 + category),
+    β = [1, 0.5],
+);
+
+erp = UnfoldSim.simulate_component(MersenneTwister(1), signal, design);
+
+# After simulation, we are ready to plot it. We expect that the simulated responses are scaled by the design's duration. To show it more effectively, we sort by duration.
+##---
+f = Figure()
+df = UnfoldMakie.eeg_array_to_dataframe(erp')
+df.duration = repeat(generate_events(design).duration, inner = size(erp, 1))
+df.category = repeat(generate_events(design).category, inner = size(erp, 1))
+plot_erp!(
+    f[1, 1],
+    df,
+    mapping = (; group = :group => nonnumeric, col = :category), #  color = :duration, fails right nowUnfoldMakie#353
+    layout = (; legend_position = :left),
+    colorbar = (; label = "Duration"),
+)
+plot_erpimage!(
+    f[2, 1],
+    erp,
+    sortvalues = generate_events(design).duration,
+    layout = (; legend_position = :bottom),
+)
+f
+
+# The scaling by the two `condition` effect levels and the modified event duration by the `duration` are clearly visible

docs/literate/HowTo/getGroundTruth.jl

@@ -0,0 +1,116 @@
+# # Simulate ground truth marginal effects
+
+# Often when testing some algorithm, we want to compare our results to a known ground truth. In the case of marginal effects via the `Unfold.effects`/ `Effects.jl` interface, we can do this using an `EffectsDesign`.
+# You can find more on what marginalized effects are here in the [Unfold.jl documentation](https://unfoldtoolbox.github.io/Unfold.jl/dev/generated/HowTo/effects/).
+
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+## Load required packages
+using UnfoldSim
+using Unfold
+using CairoMakie
+using UnfoldMakie
+using Random
+# ```@raw html
+# </details >
+# ```
+# ## Simulation
+# First let's set up a SingleSubject simulation
+
+# !!! note
+#     An `EffectsDesign` for a `MultiSubjectDesign` is not implemented yet.
+
+design =
+    SingleSubjectDesign(;
+        conditions = Dict(
+            :condition => ["bike", "face"],
+            :continuous => range(0, 5, length = 10),
+        ),
+    ) |> x -> RepeatDesign(x, 5);
+
+# **n170** has a condition effect, faces are more negative than bikes
+n1 = LinearModelComponent(;
+    basis = n170(),
+    formula = @formula(0 ~ 1 + condition),
+    β = [5, 3],
+);
+# **p300** has a continuous effect, higher continuous values will result in larger P300's.
+# We include both a linear and a quadratic effect of the continuous variable.
+p3 = LinearModelComponent(;
+    basis = p300(),
+    formula = @formula(0 ~ 1 + continuous + continuous^2),
+    β = [5, 1, 0.2],
+);
+
+components = [n1, p3]
+data, evts = simulate(
+    MersenneTwister(1),
+    design,
+    components,
+    UniformOnset(; width = 0, offset = 1000),
+    PinkNoise(),
+);
+
+# ## Simulate marginal effects directly
+# To simulate marginal effects we first have to specify an effects dictionary and subsequently hand this dict plus the original design to `EffectsDesign()`
+
+effects_dict = Dict(:condition => ["bike", "face"])
+
+effects_design = EffectsDesign(design, effects_dict)
+
+# !!! note
+#     We only specified the condition levels here, by default every unspecified variable will be set to a "typical" (i.e. the mean) value.
+
+# And finally we can simulate our ground truth marginal effects
+
+gt_data, gt_events = simulate(
+    MersenneTwister(1),
+    effects_design,
+    components,
+    NoOnset(),
+    NoNoise(),
+    return_epoched = true,
+);
+@show gt_events
+
+# Additionally, we can get the simulated effects into a tidy dataframe using Unfold's `result_to_table`.
+# Note that the data has to be reshaped into a channel X times X predictor form. (In our one channel example `size(gt_data) = (45,2)`, missing the channel dimension.)
+
+g = reshape(gt_data, 1, size(gt_data)...)
+times = range(1, size(gt_data, 1));
+gt_effects = Unfold.result_to_table([g], [gt_events], [times], ["effects"])
+first(gt_effects, 5)
+
+
+# ## Compare with Unfold.jl results
+
+m = fit(
+    UnfoldModel,
+    [
+        Any => (
+            @formula(0 ~ 1 + condition + spl(continuous, 4)),
+            firbasis(τ = [-0.1, 1], sfreq = 100, name = "basis"),
+        ),
+    ],
+    evts,
+    data,
+);
+
+ef = effects(effects_dict, m);
+
+# !!! note
+#      The ground truth signal is shorter because the ground truth typically returns values between `[0 maxlength(components)]`, whereas in our Unfold model we included a baseline period of 0.1s.
+#      If you want to actually compare results with the ground truth, you could either us `UnfoldSim.pad_array()` or set the Unfold modelling window to `τ=[0,1]`.
+
+gt_effects.type .= "UnfoldSim effects"
+ef.type .= "Unfold effects"
+
+gt_effects.time = gt_effects.time ./ 100 .- 1 / 100
+ef.continuous .= 2.5 # needed to be able to easily merge the two dataframes
+comb = vcat(gt_effects, ef)
+plot_erp(comb; mapping = (; color = :type, col = :condition))
+
+# The simulated ground truth marginal effects, and the fitted marginal effects look similar as expected, but the fitted has some additional noise because of finite data (also as expected).

docs/literate/HowTo/newComponent.jl

@@ -1,13 +1,18 @@
 # # Define a new component (with variable duration and shift)
 
 # We want a new component that changes its duration and shift depending on a column in the event design. This is somewhat already implemented in the HRF + Pupil bases.
+# !!! hint
+#     if you are just interested to use duration-dependency in your simulation, check out the component-function tutorial
+
+
 
 # ### Setup
 # ```@raw html
 # <details>
 # <summary>Click to expand</summary>
 # ```
 using UnfoldSim
+import UnfoldSim.simulate_component
 using Unfold
 using Random
 using DSP
@@ -39,7 +44,7 @@ end
 Base.length(c::TimeVaryingComponent) = length(c.maxlength)
 
 # While we could have put the TimeVaryingComponent.basisfunction directly into the simulate function, I thought this is a bit more modular
-function UnfoldSim.simulate(rng, c::TimeVaryingComponent, design::AbstractDesign)
+function UnfoldSim.simulate_component(rng, c::TimeVaryingComponent, design::AbstractDesign)
     evts = generate_events(design)
     return c.basisfunction(evts, c.maxlength)
 end
@@ -62,7 +67,7 @@ function basis_shiftduration(evts, maxlength)
 end
 
 # ## Simulate data with the new component type
-erp = UnfoldSim.simulate(
+erp = UnfoldSim.simulate_component(
     MersenneTwister(1),
     TimeVaryingComponent(basis_shiftduration, 50),
     design,

docs/literate/HowTo/newDesign.jl

@@ -3,8 +3,10 @@
 # A design specifies how much data is generated, and how the event-table(s)
 # should be generated. Already implemented examples are `MultiSubjectDesign` and `SingleSubjectDesign`.
 
+
 # We need 3 things for a new design: a `struct<:AbstractDesign`, a `size` and a `generate_events` function.
 
+
 # ### Setup
 # ```@raw html
 # <details>
@@ -15,6 +17,7 @@ using StableRNGs
 using DataFrames
 using Parameters
 using Random
+
 # ```@raw html
 # </details>
 # <br />

docs/literate/HowTo/sequence.jl

@@ -0,0 +1,93 @@
+# # Sequence of events (e.g. SCR)
+
+# In this HowTo you will learn to simulate a "SR"-Sequence, a stimulus response, followed by a button press response. 
+
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+## Load required packages
+using UnfoldSim
+using CairoMakie
+using StableRNGs
+# ```@raw html
+# </details >
+# <br />
+# ```
+
+# ## Create sequence design
+# First we generate the minimal design of the experiment by specifying our conditions (a one-condition-two-levels design in our case)
+design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"]))
+generate_events(design)
+# Next we use the `SequenceDesign` and nest our initial design in it. "`SR_`" is code for an "`S`" (stimulus) event and an "`R`" (response) event - only single letter events are supported! The "`_`" is a signal for the onset generator to generate a bigger pause - no overlap between adjacent "`SR`" pairs.
+design = SequenceDesign(design, "SR_")
+generate_events(StableRNG(1), design)
+# The main thing that happened is that the design was repeated for every event (each 'letter') of the sequence, and an `event` column was added.
+
+# !!! hint
+#     More advanced sequences are possible as well, like "SR{1,3}", or "A[BC]". Infinite sequences are **not** possible like "AB*". 
+
+# Finally, let's repeat the current design 4 times
+design = RepeatDesign(design, 4)
+generate_events(StableRNG(1), design)
+
+# This results in 16 trials that nicely follow our sequence.
+
+# !!! hint
+#     There is a difference between `SequenceDesign(RepeatDesign)` and `RepeatDesign(SequenceDesign)` for variable sequences e.g. "A[BC]", where in the former case,  one sequence is drawn e.g. "AC" and applied to all repeated rows, in the latter, one sequence for each repeat is drawn.
+
+# ## Specify components for sequence events
+# Next we have to specify for both events `S` and `R` what the responses should look like.
+p1 = LinearModelComponent(;
+    basis = p100(),
+    formula = @formula(0 ~ 1 + condition),
+    β = [1, 0.5],
+)
+
+n1 = LinearModelComponent(;
+    basis = n170(),
+    formula = @formula(0 ~ 1 + condition),
+    β = [1, 0.5],
+)
+
+p3 = LinearModelComponent(;
+    basis = UnfoldSim.hanning(Int(0.5 * 100)), # sfreq = 100 for the other bases
+    formula = @formula(0 ~ 1 + condition),
+    β = [1, 0],
+)
+
+resp = LinearModelComponent(;
+    basis = UnfoldSim.hanning(Int(0.5 * 100)), # sfreq = 100 for the other bases
+    formula = @formula(0 ~ 1 + condition),
+    β = [1, 2],
+    offset = -10,
+)
+nothing ## hide
+
+# We combine them into a dictionary with a sequence-`Char` as key 
+components = Dict('S' => [p1, n1, p3], 'R' => [resp])
+
+# ## Simulate data
+# Given the design and the components. we specify onset and noise and simulate data
+
+
+data, evts = simulate(
+    StableRNG(1),
+    design,
+    components,
+    UniformOnset(offset = 40, width = 10),
+    NoNoise(),
+)
+nothing ## hide
+
+# Finally we can plot the results
+f, ax, h = lines(data)
+vlines!(ax, evts.latency[evts.event .== 'S'], color = (:darkblue, 0.5))
+vlines!(ax, evts.latency[evts.event .== 'R'], color = (:darkred, 0.5))
+ax.xlabel = "Time [samples]"
+ax.ylabel = "EEG [a.u]"
+xlims!(ax, 0, 500)
+f
+
+# As visible, the `R` response always follows the `S` response. Due to the "`_`" we have large breaks between the individual sequences.

docs/literate/reference/designtypes.jl

@@ -40,6 +40,7 @@ design_single_shuffled = SingleSubjectDesign(;
         :stimulus_type => ["natural", "artificial"],
         :contrast_level => range(0, 1, length = 3),
     ),
+
     event_order_function = shuffle,
 );
 # ```@raw html

docs/literate/reference/onsettypes.jl

@@ -1,6 +1,6 @@
 # # Overview: Onset types
 # The onset types determine the distances between event onsets in the continuous EEG signal. The distances are sampled from a certain probability distribution.
-# Currently, there are two types of onset distributions implemented: `UniformOnset` and `LogNormalOnset`.
+# Currently, there are two types of onset distributions implemented: `UniformOnset` and `LogNormalOnset`. Both are accompanied by their `UniformOnsetFormula` and `LogNormalOnsetFormula` conterparts, which allows to modify the overlap based on the design.
 
 # ### Setup
 # ```@raw html
@@ -282,3 +282,25 @@ end
 #        - if `offset` < `length(signal.basis)` -> there might be overlap, depending on the other parameters of the onset distribution
 
 # [^1]: Wikipedia contributors. (2023, December 5). Log-normal distribution. In Wikipedia, The Free Encyclopedia. Retrieved 12:27, December 7, 2023, from https://en.wikipedia.org/w/index.php?title=Log-normal_distribution&oldid=1188400077# 
+
+
+
+
+# ## Design-dependent `X-OnsetFormula`
+
+# For additional control we provide `UniformOnsetFormula` and `LogNormalOnsetFormula` types, that allow to control all parameters by specifying formulas
+o = UnfoldSim.UniformOnsetFormula(
+    width_formula = @formula(0 ~ 1 + cond),
+    width_β = [50, 20],
+)
+events = generate_events(design)
+onsets = UnfoldSim.simulate_interonset_distances(MersenneTwister(42), o, design)
+
+f = Figure()
+ax = f[1, 1] = Axis(f)
+hist!(ax, onsets[events.cond .== "A"], bins = range(0, 100, step = 1), label = "cond: A")
+hist!(ax, onsets[events.cond .== "B"], bins = range(0, 100, step = 1), label = "cond: B")
+axislegend(ax)
+f
+
+# Voila - the inter-onset intervals are `20` samples longer for condition `B`, exactly as specified.