AutoRegressiveNoise <: AbstractNoise

Not implemented.

ExponentialNoise <: AbstractNoise

Type for generating noise with exponential decay in AR spectrum.

Tip: To manually create noise samples use the simulate_noise function.

Fields

• noiselevel = 1 (optional): Factor that is used to scale the noise.
• ν = 1.5 (optional): Exponential factor of AR decay "nu".

Examples

julia> noise = ExponentialNoise()
ExponentialNoise
  noiselevel: Int64 1
  ν: Float64 1.5

julia> using StableRNGs

julia> simulate_noise(StableRNG(1), noise, 5)
5-element Vector{Float64}:
  -5.325200748641231
  -3.437402125380177
   2.7852625669058884
  -1.5381022393382109
 -14.818799857226612

See also PinkNoise, RedNoise, NoNoise, WhiteNoise.

Hartmut <: AbstractHeadmodel

Type for accessing the HArtMuT model (Harmening et al., 2022), a head model which includes not only brain sources but also muscular and ocular sources.

Note: Use Hartmut() to create an instance of this head model.

Fields

• artefacual::Any: Dict with artefactual sources (e.g. muscular and ocular) containing label, leadfield, orientation and position.
• cortical::Any: Dict with cortical sources containing label, leadfield, orientation and position.
• electrodes::Any: Dict with electrode labels and their positions in 3D space.

Examples

julia> h = Hartmut()
Please cite: HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce
HArtMuT-Headmodel
227 electrodes:  (AF3,AF3h...) - hartmut.electrodes
2004 source points: (Left Middle Temporal Gyrus, posterior division,...) - hartmut.cortical
4260 source points: (Muscle_DepressorLabii_left,...) - hartmut.artefactual

In addition to UnfoldSim.jl please cite:
HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce

# Access artefactual sources
julia> h.artefactual
Dict{String, Any} with 4 entries:
  "label"       => ["Muscle_DepressorLabii_left", "Muscle_DepressorLabii_left", "Muscle_DepressorLabii_left", "Muscle_DepressorLabii_left", "Muscle_DepressorLabii_left", "Muscle_DepressorLabii_left", "Muscle_DepressorLabii_left", "Mu…
  "leadfield"   => [-0.00148682 -0.00229078 … -0.307231 -0.321305; 0.0141537 0.0134523 … -0.141177 -0.140583; … ; 0.108747 0.1091 … 0.106614 0.117126; 0.0257134 0.0248186 … 0.00105064 -0.00433068;;; 0.0898798 0.0891799 … 2.02466 1.81…
  "orientation" => [0.54244 0.482395 -0.687789; 0.546949 0.507161 -0.666059; … ; 0.193693 -0.979684 0.0519768; 0.327641 -0.944196 -0.0338583]
  "pos"         => [-29.3728 58.6061 -120.829; -28.4183 59.0185 -121.448; … ; -21.8088 70.5968 -28.7679; -21.1545 70.1282 -31.4184]

LinearModelComponent <: AbstractComponent

A multiple regression component for one subject.

All fields can be named. Works best with SingleSubjectDesign.

Fields

• basis::Any: an object, if accessed, provides a 'basis function', e.g. hanning(40)::Vector, this defines the response at a single event. It will be weighted by the model prediction. Future versions will allow for functions, as of v0.3 this is restricted to array-like objects
• formula::Any: StatsModels formula object, e.g. @formula 0 ~ 1 + cond (left-hand side must be 0).
• β::Vector Vector of betas/coefficients, must fit the formula.
• contrasts::Dict (optional): Determines which coding scheme to use for which categorical variables. Default is empty which corresponds to dummy coding. For more information see https://juliastats.org/StatsModels.jl/stable/contrasts.

Examples

julia> LinearModelComponent(;
           basis = hanning(40),
           formula = @formula(0 ~ 1 + cond),
           β = [1., 2.],
           contrasts = Dict(:cond => EffectsCoding())
       )
LinearModelComponent
  basis: Array{Float64}((40,)) [0.0, 0.006474868681043577, 0.02573177902642726, 0.0572719871733951, 0.10027861829824952, 0.1536378232452003, 0.21596762663442215, 0.28565371929847283, 0.3608912680417737, 0.43973165987233853  …  0.43973165987233853, 0.3608912680417737, 0.28565371929847283, 0.21596762663442215, 0.1536378232452003, 0.10027861829824952, 0.0572719871733951, 0.02573177902642726, 0.006474868681043577, 0.0]
  formula: StatsModels.FormulaTerm{StatsModels.ConstantTerm{Int64}, Tuple{StatsModels.ConstantTerm{Int64}, StatsModels.Term}}
  β: Array{Float64}((2,)) [1.0, 2.0]
  contrasts: Dict{Symbol, EffectsCoding}

See also MixedModelComponent, MultichannelComponent.

LogNormalOnset <: AbstractOnset

Log-normal inter-event distances (in samples) using the Distributions.jl truncated LogNormal distribution (code and mathematical reference).

Be careful with large μ and σ values, as they are on logscale. σ>8 can quickly give you out-of-memory sized signals!
Tip: To manually generate inter-event distance samples use the simulate_interonset_distances function.

Fields

• μ: The mean of the log-transformed variable (the log-normal random variable's logarithm follows a normal distribution).
• σ: The standard deviation of the log-transformed variable.
• offset = 0 (optional): The minimal distance between events.
• truncate_upper = nothing (optional): Upper limit (in samples) at which the distribution is truncated.

Examples

julia> onset_distribution = LogNormalOnset(3, 0.25, 10, 25)
LogNormalOnset
  μ: Int64 3
  σ: Float64 0.25
  offset: Int64 10
  truncate_upper: Int64 25

See also UniformOnset, NoOnset.

MixedModelComponent <: AbstractComponent

A component that adds a hierarchical relation between parameters according to a Linear Mixed Model (LMM) defined via MixedModels.jl.

All fields can be named. Works best with MultiSubjectDesign.

Fields

• basis::Any: an object, if accessed, provides a 'basis function', e.g. hanning(40)::Vector, this defines the response at a single event. It will be weighted by the model prediction. Future versions will allow for functions, as of v0.3 this is restricted to array-like objects
• formula::Any: Formula-object in the style of MixedModels.jl e.g. @formula 0 ~ 1 + cond + (1|subject). The left-hand side is ignored.
• β::Vector Vector of betas (fixed effects), must fit the formula.
• σs::Dict Dict of random effect variances, e.g. Dict(:subject => [0.5, 0.4]) or to specify correlation matrix Dict(:subject=>[0.5,0.4,I(2,2)],...). Technically, this will be passed to the MixedModels.jl create_re function, which creates the θ matrices.
• contrasts::Dict (optional): Dict in the style of MixedModels.jl. Determines which coding scheme to use for which categorical variables. Default is empty which corresponds to dummy coding. For more information see https://juliastats.org/StatsModels.jl/stable/contrasts.

Examples

julia> MixedModelComponent(;
           basis = hanning(40),
           formula = @formula(0 ~ 1 + cond + (1 + cond|subject)),
           β = [1., 2.],
           σs= Dict(:subject => [0.5, 0.4]),
           contrasts=Dict(:cond => EffectsCoding())
       )
MixedModelComponent
  basis: Array{Float64}((40,)) [0.0, 0.006474868681043577, 0.02573177902642726, 0.0572719871733951, 0.10027861829824952, 0.1536378232452003, 0.21596762663442215, 0.28565371929847283, 0.3608912680417737, 0.43973165987233853  …  0.43973165987233853, 0.3608912680417737, 0.28565371929847283, 0.21596762663442215, 0.1536378232452003, 0.10027861829824952, 0.0572719871733951, 0.02573177902642726, 0.006474868681043577, 0.0]
  formula: StatsModels.FormulaTerm{StatsModels.ConstantTerm{Int64}, Tuple{StatsModels.ConstantTerm{Int64}, StatsModels.Term, StatsModels.FunctionTerm{typeof(|), Vector{StatsModels.AbstractTerm}}}}
  β: Array{Float64}((2,)) [1.0, 2.0]
  σs: Dict{Symbol, Vector{Float64}}
  contrasts: Dict{Symbol, EffectsCoding}

See also LinearModelComponent, MultichannelComponent.

MultiSubjectDesign <: AbstractDesign

A type for specifying the experimental design for multiple subjects (based on the given random-effects structure).

Tip: Check the resulting dataframe using the generate_events function.
Please note that the number of items n_items has to be a multiple of the number of between-item levels. The sample applies for n_subjects and the number of between-subject levels.

Fields

• n_subjects::Int: Number of subjects.
• n_items::Int: Number of items/stimuli (sometimes ≈ trials).
• subjects_between::Dict{Symbol,Vector}: Effects between subjects, e.g. young vs old.
• items_between::Dict{Symbol,Vector}: Effects between items, e.g. natural vs artificial images, (but shown to all subjects if not specified also in subjects_between).
• both_within::Dict{Symbol,Vector}: Effects completely crossed i.e. conditions/covariates that are both within-subject and within-item.
• event_order_function = (rng, x) -> x: Can be used to sort, or shuffle the events e.g. (rng, x) -> shuffle(rng, x) (or shorter just event_order_function = shuffle). The default is the identify function, i.e. not changing the order of the events.

Examples

# Declaring the same condition both between-subject and between-item results in a full between-subject/item design.
julia> design = MultiSubjectDesign(;
                       n_items = 10,
                       n_subjects = 30,
                       subjects_between = Dict(:cond => ["levelA", "levelB"]),
                       items_between = Dict(:cond => ["levelA", "levelB"]),
                       )
MultiSubjectDesign
  n_subjects: Int64 30
  n_items: Int64 10
  subjects_between: Dict{Symbol, Vector}
  items_between: Dict{Symbol, Vector}
  both_within: Dict{Symbol, Vector}
  event_order_function: #3 (function of type UnfoldSim.var"#3#7")

julia> generate_events(StableRNG(1), design)
150×3 DataFrame
 Row │ subject  cond    item   
     │ String   String  String 
─────┼─────────────────────────
   1 │ S01      levelA  I01
   2 │ S01      levelA  I03
   3 │ S01      levelA  I05
  ⋮  │    ⋮       ⋮       ⋮
 148 │ S30      levelB  I06
 149 │ S30      levelB  I08
 150 │ S30      levelB  I10
               144 rows omitted

See also SingleSubjectDesign, RepeatDesign

MultichannelComponent <: AbstractComponent

Projects a AbstractComponent to multiple "channels" via the projection vector.

Optionally, noise can be added to the source prior to projection. By default a MultichannelComponent can be constructed using one of the following options for projection:

• projection::AbstractVector: Directly pass a custom projection vector.
• projection::Pair{<:AbstractHeadmodel,String}: Generate a projection vector by specifying which headmodel to use and which sources should be active.

Fields

• component::AbstractComponent: The component that should be projected to the sensors.
• projection::AbstractVector or projection::Pair{<:AbstractHeadmodel,String}: Vector p that projects the (source) component c[t] (where t is time) to the sensors s. The length of p equals the number of sensors s. Typically, it is a slice of the leadfield matrix. out[s,t] = p[s]*c[t].
• noise::AbstractNoise (optional): Noise added in the source space. Default is NoNoise.

Examples

# Variant 1: Specify the projection vector manually
julia> c1 = LinearModelComponent(; basis = p100(), formula = @formula(0 ~ 1), β = [1]);

julia> mc1 = UnfoldSim.MultichannelComponent(c, [1, 2, -1, 3, 5, 2.3, 1])
MultichannelComponent
  component: LinearModelComponent
  projection: Array{Float64}((7,)) [1.0, 2.0, -1.0, 3.0, 5.0, 2.3, 1.0]
  noise: NoNoise NoNoise()

# Variant 2: Use a headmodel and specify a source
julia> c2 = LinearModelComponent(; basis = p300(), formula = @formula(0 ~ 1), β = [1]);

julia> hart = headmodel(type = "hartmut");
Please cite: HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce

julia> mc2 = UnfoldSim.MultichannelComponent(c2, hart => "Right Occipital Pole")
MultichannelComponent
  component: LinearModelComponent
  projection: Array{Float64}((227,)) [-0.03461859471337842, -0.04321094803502425, 0.0037088347968313525, -0.014722528968861278, -0.0234889834534478, 0.02731807504242923, 0.038863688452528036, 0.1190531258070562, -0.09956890221613562, -0.0867729334438599  …  0.37435404409695094, -0.020863789022627935, 0.25627478723535513, -0.05777985212119245, 0.37104376432271147, -0.19446620423767172, 0.2590764703721097, -0.12923837607416555, 0.1732886690359311, 0.4703016561960567]
  noise: NoNoise NoNoise()

See also LinearModelComponent, MixedModelComponent.

NoNoise <: AbstractNoise

A type for simulations without noise; return zeros instead of noise.

Tip: To manually create noise samples use the simulate_noise function.

Examples

julia> noise = NoNoise()
NoNoise()

julia> using StableRNGs

julia> simulate_noise(StableRNG(1), noise, 3)
3-element Vector{Float64}:
 0.0
 0.0
 0.0

See also PinkNoise, RedNoise, ExponentialNoise, WhiteNoise.

NoOnset <: AbstractOnset

For cases where the user wants to simulate epoched data without any overlap between consecutive events.

Examples

julia> onset_distribution = NoOnset()
NoOnset()

See also UniformOnset, LogNormalOnset.

PinkNoise <: AbstractNoise

A type for generating Pink Noise using the SignalAnalysis.jl implementation.

The noise values are sampled from a standard normal distribution 𝒩(μ=0, σ=1). That means that ~95% of the values are between ~-2 and ~2 (with noiselevel = 1).
Tip: To manually create noise samples use the simulate_noise function.

Fields

• noiselevel = 1 (optional): Factor that is used to scale the noise.
• func = SignalAnalysis.PinkGaussian (optional): Function that is used to create the noise samples. This field is for internal use and should not typically be modified directly by the user. Changes to this field may result in unexpected behavior.

Examples

julia> noise = PinkNoise(noiselevel = 3)
PinkNoise
  noiselevel: Int64 3
  func: UnionAll

julia> using StableRNGs

julia> simulate_noise(StableRNG(1), noise, 5)
5-element Vector{Float64}:
 2.578878369756878
 3.4972108606501786
 2.878568584946028
 2.2725654770788384
 3.5291669151888683

See also RedNoise, WhiteNoise, ExponentialNoise, NoNoise.

RealisticNoise <: AbstractNoise

Not implemented - planned to use Artifacts.jl to provide real EEG data to add.

RedNoise <: AbstractNoise

A type for generating Red Noise using the SignalAnalysis.jl implementation.

The noise values are sampled from a standard normal distribution 𝒩(μ=0, σ=1). That means that ~95% of the values are between ~-2 and ~2 (with noiselevel = 1).
Tip: To manually create noise samples use the simulate_noise function.

Fields

• noiselevel = 1 (optional): Factor that is used to scale the noise.
• func = SignalAnalysis.RedGaussian (optional): Function that is used to create the noise samples. This field is for internal use and should not typically be modified directly by the user. Changes to this field may result in unexpected behavior.

Examples

julia> noise = RedNoise(noiselevel = 2)
RedNoise
  noiselevel: Int64 2
  func: UnionAll

julia> using StableRNGs

julia> simulate_noise(StableRNG(2), noise, 3)
3-element Vector{Float64}:
 -0.34153942884005967
 -0.4651387715669636
 -0.4951538876376382

See also PinkNoise, WhiteNoise, ExponentialNoise, NoNoise.

RepeatDesign{T} <: AbstractDesign

Repeat a design (and the corresponding events DataFrame) multiple times to mimick repeatedly recorded trials.

Tip: Check the resulting dataframe using the generate_events function. Please note that when using an event_order_function(e.g. shuffle) in a RepeatDesign, the corresponding RNG is shared across repetitions and not deep-copied for each repetition. As a result, the order of events will differ for each repetition.

Fields

• design::T: The experimental design that should be repeated.
• repeat::Int = 1: The number of repetitions.

Examples

julia> using StableRNGs # For using the `generate_events` function in a reproducible way

julia> design_once =
           SingleSubjectDesign(;
               conditions = Dict(
                   :stimulus_type => ["natural", "artificial"],
                   :contrast_level => range(0, 1, length = 2),
               ),
               event_order_function = shuffle,
           );

julia> generate_events(StableRNG(1), design_once)
4×2 DataFrame
 Row │ contrast_level  stimulus_type 
     │ Float64         String        
─────┼───────────────────────────────
   1 │            1.0  natural
   2 │            1.0  artificial
   3 │            0.0  natural
   4 │            0.0  artificial

julia> design = RepeatDesign(design_once, 2)
RepeatDesign{SingleSubjectDesign}
  design: SingleSubjectDesign
  repeat: Int64 2

julia> generate_events(StableRNG(1), design)
8×2 DataFrame
 Row │ contrast_level  stimulus_type 
     │ Float64         String        
─────┼───────────────────────────────
   1 │            1.0  natural
   2 │            1.0  artificial
   3 │            0.0  natural
   4 │            0.0  artificial
   5 │            1.0  artificial
   6 │            0.0  natural
   7 │            1.0  natural
   8 │            0.0  artificial

See also SingleSubjectDesign, MultiSubjectDesign

Simulation

A type to store all "ingredients" for a simulation including their parameters.

Can either be created by the user or will be created automatically when calling the simulate function with the required "ingredients". Tip: Use the subtypes function to get an overview of the implemented "ingredients", e.g. subtypes(AbstractDesign).

Fields

• design::AbstractDesign: Experimental design.
• components::Vector{AbstractComponent}: Response function(s) for the events (e.g. the ERP shape in EEG).
• onset::AbstractOnset: Inter-onset distance distribution.
• noisetype::AbstractNoise: Noise type.

Examples

julia> design = SingleSubjectDesign(; conditions = Dict(:stimulus_type => ["natural", "artificial"]));

julia> component = LinearModelComponent(;
           basis = p100(),
           formula = @formula(0 ~ 1 + stimulus_type),
           β = [2, 3],
       );

julia> onset = UniformOnset();

julia> noise = PinkNoise();

julia> simulation = Simulation(design, component, onset, noise);

julia> using StableRNGs

julia> data, events = simulate(StableRNG(1), simulation);

julia> events
2×2 DataFrame
 Row │ stimulus_type  latency 
     │ String         Int64   
─────┼────────────────────────
   1 │ natural             20
   2 │ artificial          70

julia> data
85-element Vector{Float64}:
  0.8596261232522926
  1.1657369535500595
  0.9595228616486761
  ⋮
  0.9925202143746904
  0.2390652543395527
 -0.11672523788068771

SingleSubjectDesign <: AbstractDesign

A type for specifying the experimental design for a single subject (based on the given conditions).

Tip: Check the resulting dataframe using the generate_events function.
The number of trials/rows in the output of generate_events([rng, ]design) depends on the full factorial of your conditions.
To increase the number of repetitions, e.g. by 5, simply use RepeatDesign(SingleSubjectDesign(...), 5).
If conditions are omitted (or set to nothing), a single trial is simulated with a column :dummy and content :dummy - this is for convenience.

Fields

• conditions::Dict{Symbol,Vector}: Experimental conditions, e.g. Dict(:A => ["a_small","a_big"], :B => ["b_tiny","b_large"]).
• event_order_function = (rng, x) -> x: Can be used to sort by specifying sort, or shuffling by providing shuffle, or custom functions following the interface (rng, x) -> my_shuffle(rng,x). The default is the identify function, i.e. not changing the order of the events.

Examples

julia> using StableRNGs # For using the `generate_events` function in a reproducible way

julia> design =
           SingleSubjectDesign(;
               conditions = Dict(
                   :stimulus_type => ["natural", "artificial"],
                   :contrast_level => range(0, 1, length = 5),
               ),
           )
SingleSubjectDesign
  conditions: Dict{Symbol, Vector}
  event_order_function: #10 (function of type UnfoldSim.var"#10#14")

julia> generate_events(StableRNG(1), design)
10×2 DataFrame
 Row │ contrast_level  stimulus_type 
     │ Float64         String        
─────┼───────────────────────────────
   1 │           0.0   natural
   2 │           0.25  natural
   3 │           0.5   natural
  ⋮  │       ⋮               ⋮
   8 │           0.5   artificial
   9 │           0.75  artificial
  10 │           1.0   artificial
                       4 rows omitted

See also MultiSubjectDesign, RepeatDesign

UniformOnset <: AbstractOnset

Provide a Uniform Distribution for the inter-event distances (in samples).

Tip: To manually generate inter-event distance samples use the simulate_interonset_distances function.

Fields

• width = 50 (optional): Width of the uniform distribution (=> the "jitter"). Since the lower bound is 0, width is also the upper bound.
• offset = 0 (optional): The minimal distance between events. The maximal distance is offset + width.

Examples

julia> onset_distribution = UniformOnset(width = 25, offset = 5)
UniformOnset
  width: Int64 25
  offset: Int64 5

See also LogNormalOnset, NoOnset.

WhiteNoise <: AbstractNoise

A type for generating White Noise using randn - thus Gaussian noise.

The noise values are sampled from a standard normal distribution 𝒩(μ=0, σ=1). That means that ~95% of the values are between ~-2 and ~2 (with noiselevel = 1).
Tip: To manually create noise samples use the simulate_noise function.

Fields

• noiselevel = 1 (optional): Factor that is used to scale the noise.
• imfilter = nothing (optional): Use imfilter > 0 to smooth the noise using Image.imfilter with a Gaussian kernel with σ = imfilter.

Examples

julia> noise = WhiteNoise()
WhiteNoise
  noiselevel: Int64 1
  imfilter: Nothing nothing

julia> using StableRNGs

julia> simulate_noise(StableRNG(1), noise, 3)
3-element Vector{Float64}:
 -0.5325200748641231
  0.098465514284785
  0.7528865221245234

See also PinkNoise, RedNoise, ExponentialNoise, NoNoise.

Length is the product of all dimensions and equals the number of events in the corresponding events dataframe.

Return the dimensions of the experiment design.

hanning(duration, offset, sfreq)

Generate a (potentially shifted) hanning window with a certain duration.

Note: This function extends the DSP.hanning function using multiple dispatch.

Arguments

• duration: in s.
• offset: in s, defines hanning peak i.e. shift of the hanning window.
• sfreq: Sampling rate in Hz.

Returns

• Vector: Contains a shifted (i.e. zero-padded) hanning window.

Examples

julia> UnfoldSim.hanning(0.1, 0.3, 100)
25-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.9698463103929542
 0.75
 0.41317591116653485
 0.11697777844051105
 0.0

PuRF(; n = 10.1, tmax = 0.93, sfreq = 100)

Default generator for PuRF Pupil Response Function. The canonical PRF is a gamma function and implemented according to Denison 2020 equation (2) going back to Hoeks & Levelt, 1993.

The pupil response is evaluated at t = 0:1/sfreq:3*tmax. The response is normalized by the peak-maximum at tmax, thus typically a pupil-response of 1 is returned (disregarding numerical issues).

Keyword arguments:

• n = 10.1: shape parameter
• tmax = 0.93: peak maximum
• sfreq = 100: sampling frequency

Returns

• Vector: canonical pupil response with length(0:1/sfreq:3*tmax) entries.

Examples

julia> PuRF(; n = 5)
280-element Vector{Float64}:
 0.0
 2.0216617815131253e-8
 6.130689396024061e-7
 4.4118063684811444e-6
 1.761817666835793e-5
 ⋮
 0.012726253506722554
 0.012280989091455786
 0.011850525657416842
 0.011434405338911133
 0.011032182932283816

add_noise!(rng, noisetype::AbstractNoise, signal)

Generate and add noise to a signal.

Assumes that the signal can be linearized, that is, that the noise is stationary.

add_responses!(signal, responses::Vector, e, s, tvec, trial_idx)
add_responses!(signal, responses::Matrix, e, s, tvec, trial_idx)
add_responses!(signal, responses::AbstractArray, e, s, tvec, trial_idx)

Add (in-place) the given responses to the signal, for both 2D (1 channel) and 3D (X channel case). Helper function.

Arguments

• signal: Continuous EEG signal to be modified in place. Has the dimensions channels x continuous_time x subjects.
• responses::Union{Vector, Matrix, AbstractArray}: Responses to be added. In the multi-channel case, the dimensions are channels x maxlength(components) x length(simulation.design), else maxlength(components) x length(simulation.design). The data for all the subjects and their respective trials is concatenated.
• e: Index of the channel (in signal) for which to add the response.
• s: Index of the subject (in signal) for which to add the response.
• tvec: Time points (indices in signal) at which to add the response.
• trial_idx: Index of the particular trial (in responses) from where the response is to be added.

Returns

• Nothing: signal is modified in-place.

Examples

julia> signal, responses, tvec = zeros(5,15,2), ones(5,6), 1:5;

julia> UnfoldSim.add_responses!(signal, responses, 1, 2, tvec, 5);

julia> signal
5×15×2 Array{Float64, 3}:
[:, :, 1] =
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

[:, :, 2] =
 1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

apply_event_order_function(fun, rng, events)

Apply fun(rng, events), raise an error if function fun is wrongly defined. Convenience function to not repeat the error handling at multiple places.

closest_src(coords_list::AbstractVector{<:AbstractVector}, pos)

When applied to a vector of target coordinate vectors, return the index of the closest position in pos for each of the targets.

Returns

• Vector: Index of the closest position in pos for each vector in coords_list.

Examples

julia> positions = [0 0 0; 2 2 2; 3 3 3; 4 4 4];
julia> target_coordinate_vectors = [[0.5, 0, 0.5], [10, 1, 8]];

julia> UnfoldSim.closest_src(target_coordinate_vectors, positions)
2-element Vector{Int64}:
 1
 4

closest_src(head::Hartmut, label::String)

Return the index of the Harmut model source point that is closest (using Euclidean distance) to the mean position of the source points matching the given label.

Arguments

• head::Hartmut: Headmodel of type Hartmut.
• label::String: Label of the source point(s) of interest.

Returns

• Int: Index of the closest Hartmut source point.

Examples

julia> hartmut = Hartmut();
Please cite: HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce

julia> label = "Right Cingulate Gyrus, posterior division";

julia> UnfoldSim.closest_src(hartmut, label)
1875

closest_src(coords::Vector{<:Real}, pos)

Return the index of the position that is closest to the target coordinates (using the Euclidean distance).

Arguments

• coords::Vector{<:Real}: Target coordinate vector (typically with length 3).
• pos: Matrix that contains all possible positions. Its dimensions are the number of positions n times number of entries in the position coordinate vectors (usually 3) i.e. n x 3.

Returns

• Int: Row index of the position in pos that is closest to coords.

Examples

julia> positions = [0 0 0; 2 2 2; 3 3 3; 4 4 4]
4×3 Matrix{Int64}:
 0  0  0
 2  2  2
 3  3  3
 4  4  4

julia> target_coordinates = [0.5, 0, 0.5]
3-element Vector{Float64}:
 0.5
 0.0
 0.5

 julia> UnfoldSim.closest_src(target_coordinates, positions)
1

create_continuous_signal(rng, responses, simulation)

Simulate onset latencies and add together a continuous signal, based on the given responses and simulation parameters. Helper function.

Arguments

• rng: Random number generator, important to ensure reproducibility.
• responses: Responses to be combined with the given onsets.
• simulation: Simulation parameters, including design, components, onsets, and noisetype.

Returns

• (signal, latencies)::Tuple{Array, Array}:
  • signal contains the generated signal. Has the dimensions channels x continuous_time x subjects.
  • latencies contains the onset latencies.

Examples

julia> using StableRNGs # to get an RNG

julia> design = SingleSubjectDesign(; conditions = Dict(:cond => ["natural", "artificial"]));

julia> c1 = LinearModelComponent(; basis = p100(), formula = @formula(0 ~ 1 + cond), β = [1, 0.5]);

julia> c2 = LinearModelComponent(; basis = p300(), formula = @formula(0 ~ 1), β = [2]);

julia> simulation = Simulation(design, [c1, c2], UniformOnset(; width = 0, offset = 30), PinkNoise());

julia> responses = simulate_responses(StableRNG(1), [c1, c2], simulation)
45×2 Matrix{Float64}:
 0.0        0.0
 0.0        0.0
 ⋮          
 0.0233794  0.0233794
 0.0        0.0

julia> signal, latencies = UnfoldSim.create_continuous_signal(StableRNG(1), responses, simulation);

julia> signal
106-element Vector{Float64}:
 0.0
 0.0
 0.0
 ⋮
 0.023379444289913343
 0.0
 0.0

 julia> latencies
2-element Vector{Int64}:
 31
 61

epoch(data::AbstractVector, args...; kwargs...)

One channel case. The data is reshaped and then passed to the multi-channel epoch method.

epoch(
    data::AbstractArray{T,2},
    events,
    τ::Tuple{Number,Number},
    sfreq;
    eventtime::Symbol = :latency,
) where {T<:Union{Missing,Number}}

Helper function to segment continuous data into epochs based on the given events and the time window τ.

Adapted from Unfold.jl: https://github.com/unfoldtoolbox/Unfold.jl/blob/b3 a21c2bb7e93d2f45ec64b0197f4663a6d7939a/src/utilities.jl#L40

Arguments

• data::AbstractArray{T,2}: Continuous data with the dimensions channels x continuous_time.
• events: Events data frame (based on the design) with latencies.
• τ::Tuple{Number,Number}: Time window for epoching in s.
• sfreq: Sampling frequency in Hz.

Keyword arguments

• eventtime::Symbol = :latency: The name of the column in events that contains the latencies.

Returns

• Array: Epoched data with the dimensions channels x times x event (or times x event for the single channel case).

Examples

# One channel example
julia> data, events = UnfoldSim.predef_eeg();

julia> size(data), size(events)
((120199,), (2000, 3))

julia> UnfoldSim.epoch(data, events, (-0.2, 1), 100)
121×2000 Matrix{Float64}:
  0.114127   -0.105347     …  -0.125485   1.6383
  0.128198    0.0474874        0.0112935  1.28122
  0.0547917  -0.0832977       -0.126181   0.850062
  0.0992842  -0.230224        -0.0449072  0.496583
  0.024461   -0.175023        -0.0223837  0.170389
  0.165133   -0.000793527  …  -0.0278197  0.104454
  ⋮                        ⋱              
 -0.362249    2.91297          0.546699   0.0
 -0.199265    2.58394          0.171159   0.0
 -0.184075    2.34611         -0.0269841  0.0
 -0.13901     2.11971         -0.0552873  0.0
 -0.0674085   1.74561      …  -0.187959   0.0

generate_events(design::AbstractDesign)
generate_events(rng::AbstractRNG, design::AbstractDesign)

Generate a full-factorial events DataFrame based on the experimental conditions and covariates defined in the design.

Arguments

• design::AbstractDesign: Experimental design for which the events DataFrame should be created.
• rng::AbstractRNG (optional): Random number generator (RNG) to make the process reproducible. If none is given, MersenneTwister(1) will be used.

Returns

• DataFrame: Each row corresponds to one combination of condition/covariate levels which is often equivalent to one stimulus or trial.

generate_events(design::MultiSubjectDesign)
generate_events(rng::AbstractRNG, design::MultiSubjectDesign)

Generate the events Dataframe according to MixedModelsSim.jl's simdat_crossed function.

Afterwards apply design.event_order_function and finally sort by :subject.
Note: No condition can be named dv which is used internally in MixedModelsSim / MixedModels as a dummy left-side

Examples

julia> using Random # for shuffling
julia> using StableRNGs

julia> design = MultiSubjectDesign(;
                       n_items = 4,
                       n_subjects = 5,
                       both_within = Dict(:condition => ["red", "green"]),
                       event_order_function = shuffle,
                       );

julia> generate_events(StableRNG(1), design)
40×3 DataFrame
 Row │ subject  item    condition 
     │ String   String  String    
─────┼────────────────────────────
   1 │ S1       I2      red
   2 │ S1       I2      green
   3 │ S1       I3      green
  ⋮  │    ⋮       ⋮         ⋮
  38 │ S5       I3      red
  39 │ S5       I2      red
  40 │ S5       I4      green
                   34 rows omitted

UnfoldSim.generate_events([rng::AbstractRNG, ]design::RepeatDesign{T})

For a RepeatDesign, iteratively call generate_events for the underlying {T} design and concatenate the results.

In case of MultiSubjectDesign, sort by subject.
Please note that when using an event_order_function(e.g. shuffle) in a RepeatDesign, the corresponding RNG is shared across repetitions and not deep-copied for each repetition. As a result, the order of events will differ for each repetition.

generate_events(design::SingleSubjectDesign)
generate_events(rng::AbstractRNG, design::SingleSubjectDesign)

Generate the events DataFrame based on design.conditions and afterwards apply design.event_order_function.

If design.conditions is nothing, a single trial is simulated with a column :dummy and content :dummy - this is for convenience.

Examples

julia> using Random # for shuffling
julia> using StableRNGs

julia> design = SingleSubjectDesign(;
           conditions = Dict(:A => ["small", "large"], :B => range(1, 5, length = 3)),
           event_order_function = shuffle,
       );

# Variant 1: Without specifying an RNG, MersenneTwister(1) will be used for the shuffling specified as `event_order_function`.
julia> generate_events(design)
6×2 DataFrame
 Row │ A       B       
     │ String  Float64 
─────┼─────────────────
   1 │ large       5.0
   2 │ large       1.0
   3 │ large       3.0
   4 │ small       1.0
   5 │ small       3.0
   6 │ small       5.0

# Variant 2: Use a custom RNG.
julia> generate_events(StableRNG(1), design)
6×2 DataFrame
 Row │ A       B       
     │ String  Float64 
─────┼─────────────────
   1 │ large       5.0
   2 │ large       3.0
   3 │ small       1.0
   4 │ small       3.0
   5 │ large       1.0
   6 │ small       5.0

Return the citation string for HArtMuT.

hrf(;
TR = 1,
peak = 6.0,
post_undershoot = 16,
length = 32.0,
peak_width = 1.0,
post_undershoot_width = 1,
amplitude = 6,
shift = 0)

Generate a parameterized BOLD haemodynamic response function (HRF) kernel based on gamma-functions.

Implementation and default parameters were taken from the SPM-toolbox.

Note: TR = 1/sfreq

Keyword arguments

• TR = 1: repetition time, 1/sfreq.
• length = 32.0: total length of the kernel in seconds.
• amplitude = 6: maximal amplitude.
• peak = 6.0: peak timing.
• peak_width = 1.0: width of the peak.
• post_undershoot = 16: post-undershoot timing.
• post_undershoot_width = 1: post-undershoot width.
• shift = 0: shift the whole HRF.

Returns

• Vector: HRF BOLD response.

Examples

julia> hrf()
33-element Vector{Float64}:
  0.0
  0.0007715994433635659
  0.019784004131204957
  0.08202939459091822
  0.158157713522699
  ⋮
 -0.0006784790038792572
 -0.00042675060451877775
 -0.000263494738348278
 -0.00015990722628360688
 -9.548780093799345e-5

leadfield(hart::Hartmut; type = "cortical")

Return the leadfield for the (cortical or artefacual) sources of the HArtMuT model.

Keyword arguments

• type = "cortical": Defines whether the "cortical" or "artefactual" leadfield should be returned.

Returns

• Array{Float64, 3}: Leadfield values i.e. how much each source contributes to the potential measured at each electrode. The output dimensions are electrodes x sources x spatial dimension.

Examples

julia> h = Hartmut();
Please cite: HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce

julia> lf = leadfield(h);

julia> size(lf)
(227, 2004, 3)

# Access the leadfield for one spatial dimension
julia> lf[:,:,1]
227×2004 Matrix{Float64}:
  0.151429    0.0284341  …  -0.119927     -0.158114
  0.1732      0.0441432     -0.110515     -0.165316
  0.249592    0.10857       -0.000593027  -0.0206122
  0.245206    0.104854      -0.0200251    -0.055392
  0.126496    0.0118467     -0.128248     -0.146107
  0.253092    0.111688   …   0.02277       0.0184387
  ⋮                      ⋱                
  0.20306     0.138837       0.133486      0.18334
 -0.689154   -1.00904       -0.108276     -0.105398
  0.192729    0.148448       0.0924756     0.142322
 -1.26181    -1.59936       -0.140598     -0.133054
  0.213982    0.145953   …   0.115515      0.170698
  0.0731569   0.0794415      0.0485819     0.107929

magnitude(headmodel::AbstractHeadmodel)

Return the magnitude for the given headmodel based on the leadfield (and potentially the source orientations) specified in the headmodel.

If the headmodel includes source orientations these are used in the calculations, otherwise the leadfield is returned assuming that the source orientations are already included. Please note that (at least in the case of the HArtMuT model) the leadfield for the cortical sources is used.

Returns

• Matrix{Float64}: Contribution of each source to the potential measured at each electrode taking into account the orientation of the sources. The output dimensions are electrodes x sources.

Examples

# Using the HArtMut model as an example
julia> h = Hartmut();
Please cite: HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce

julia> magnitude(h)
227×2004 Matrix{Float64}:
  0.111706   …   0.301065
  0.0774359      0.332934
  ⋮          ⋱  
 -0.164539      -0.246555

magnitude(lf::AbstractArray{T,3}, orientation::AbstractArray{T,2}) where {T<:Real}

Return the magnitude of the leadfield lf along the given orientation.

Arguments

• lf::AbstractArray{T,3}: Leadfield with the dimensions electrodes x sources x spatial dimension.
• orientation::AbstractArray{T,2}: Source orientations with the dimensions sources x spatial dimensions.

Examples

# Specify the leadfield (often given by a headmodel)
julia> lf = cat([1 0; 0 1; 0 0], [1 1; 0 0; 0.5 0.5], dims=3);

# Specify the source orientations
julia> ori = [1.0 0; 0 1]

# Calculate the magnitude
julia> magnitude(lf, ori)
3×2 Matrix{Float64}:
 1.0  1.0
 0.0  0.0
 0.0  0.5

maxlength(c::Vector{<:AbstractComponent}) = maximum(length.(c))

Return the maximum of the individual components' lengths.

n170(; sfreq = 100)

Generate a Hanning window mimicking an N170 EEG component with a negative (!) peak at 170ms and a width of 150ms.

Keyword arguments

• sfreq = 100: Sampling frequency in Hz.

Returns

• Vector: Contains a shifted (i.e. zero-padded) hanning window.

Examples

julia> n170(; sfreq = 120)
28-element Vector{Float64}:
 -0.0
 -0.0
 -0.0
 -0.0
 -0.0
  ⋮
 -0.45386582026834904
 -0.2771308221117309
 -0.1304955413896705
 -0.03376388529782215
 -0.0

See also p100, p300, n400, hanning.

n400(; sfreq = 100)

Generate a Hanning window mimicking an N400 EEG component with a negative (!) peak at 400ms and a width of 400ms.

Keyword arguments

• sfreq = 100: Sampling frequency in Hz.

Returns

• Vector: Contains a shifted (i.e. zero-padded) hanning window.

Examples

julia> n400(; sfreq = 250)
150-element Vector{Float64}:
 -0.0
 -0.0
 -0.0
 -0.0
 -0.0
  ⋮
 -0.016025649301821876
 -0.009035651368646647
 -0.00402259358460233
 -0.0010066617640578368
 -0.0

See also p100, p300, n170, hanning.

n_channels(c::AbstractComponent)

Return the number of channels for the given component c. By default = 1.

n_channels(c::MultichannelComponent)

For MultichannelComponent return the length of the projection vector.

n_channels(c::Vector{<:AbstractComponent})

For a vector of MultichannelComponents, return the number of channels for the first component but assert all are of equal length.

orientation(hart::Hartmut; type = "cortical")

Return the orientations of the (cortical or artefacual) sources of the HArtMuT model.

The norm of the orientation vectors is 1 and the values are between -1 and 1.

Keyword arguments

• type = "cortical": Defines whether the "cortical" or "artefactual" orientations should be returned.

Returns

• Matrix{Float64}: Orientations in 3D space. The output dimensions are sources x spatial dimensions.

Examples

julia> h = Hartmut();
Please cite: HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce

julia> orientation(h)
2004×3 Matrix{Float64}:
 -0.921919   0.364292   0.131744
 -0.900757  -0.415843  -0.125345
 -0.954087   0.117479  -0.27553
 -0.814613  -0.55344    0.17352
 -0.790526   0.276849  -0.546281
  ⋮                    
 -0.962905   0.20498   -0.17549
 -0.828358   0.557468  -0.0552498
 -0.963785  -0.265607  -0.0239074
 -0.953909   0.203615   0.22045
 -0.75762    0.128027   0.640016

p100(; sfreq = 100)

Generate a Hanning window mimicking a P100 EEG component with a peak at 100ms and a width of 100ms.

Keyword arguments

• sfreq = 100: Sampling frequency in Hz.

Returns

• Vector: Contains a shifted (i.e. zero-padded) hanning window.

Examples

julia> p100(; sfreq = 200)
30-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.37725725642960045
 0.22652592093878665
 0.10542974530180327
 0.02709137914968268
 0.0

See also p300, n170, n400, hanning.

p300(; sfreq = 100)

Generate a Hanning window mimicking a P300 EEG component with a peak at 300ms and a width of 300ms.

Keyword arguments

• sfreq = 100: Sampling frequency in Hz.

Returns

• Vector: Contains a shifted (i.e. zero-padded) hanning window.

Examples

julia> p300(; sfreq = 150)
67-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.07937323358440934
 0.04518400232274078
 0.02025351319275137
 0.005089279059533658
 0.0

See also p100, n170, n400, hanning.

pad_array(arr::Vector, len::Int, val)    
pad_array(arr::Vector, len::Tuple, val)

Pads the input array arr with a specified value val either before or after the existing elements, based on the sign of len.

Arguments

• arr::Vector: The input array to be padded.
• len::Union{Int, Tuple}: The number of times that val should be added. If len is negative, the values are added before the existing array. Otherwise, they are added after the existing array. If len is a tuple, pad_array is called twice which enables padding the array before and after in one function call.
• val: The value to be used for padding.

Returns

• Vector: Padded vector with the new length length(arr) + sum(abs.(len)).

Examples

# Create an array that will be padded
julia> my_array = rand(5)
5-element Vector{Float64}:
 0.0420017254437951
 0.19179144973603235
 0.5388760239550549
 0.6973699906283798
 0.9966598131018376

 # Pad the array with zeros before the original array
julia> pad_array(my_array, -2, 0)
7-element Vector{Float64}:
 0.0
 0.0
 0.0420017254437951
 0.19179144973603235
 0.5388760239550549
 0.6973699906283798
 0.9966598131018376

# Pad the array with the value 5 before and after the original array
julia> pad_array(my_array, (-2, 1), 5)
8-element Vector{Float64}:
 5.0
 5.0
 0.0420017254437951
 0.19179144973603235
 0.5388760239550549
 0.6973699906283798
 0.9966598131018376
 5.0

predef_2x2(rng::AbstractRNG; <keyword arguments>)

Simulate data for a 2x2 design i.e. a design with two conditions with two levels each.

Note that this function is mainly used for demonstration and internal testing purposes or whenever a quick simulation is needed.

The most used keyword argument is: return_epoched = true which returns already epoched data. If you want epoched data without overlap, specify onset = NoOnset() and return_epoched = true.

Be careful if you modify n_items with n_subjects = 1, n_items has to be a multiple of 4 (or your equivalent conditions factorial, e.g. all combinations length).

In difference to predef_EEG, predef_2x2 is sample based (no sampling rate to be specified), and also has a 2x2 design, instead of a 2-categorical, 1-continuous design.

Keyword arguments

Design

• n_items = 100: Number of items.
• n_subjects = 1: Number of subjects.
• conditions = Dict(:A => ["a_small","a_big"], :B => ["b_tiny","b_large"]): Experimental conditions with their levels.
• event_order_function = shuffle: Random trial order.

Component / Signal

• signalsize = 100: Length of simulated hanning window.
• basis = hanning(signalsize): The actual "function". signalsize is only used here.
• β = [1, -0.5, .5, +1]: The parameters for the fixed effects.
• σs = Dict(:subject => [1, 0.5, 0.5, 0.5],:item => [1]): Only needed in n_subjects >= 2 cases; specifies the random effects.
• contrasts = Dict(:A => EffectsCoding(), :B => EffectsCoding()): Effect coding by default.
• formula = n_subjects == 1 ? @formula(0 ~ 1 + A*B) : @formula(dv ~ 1 + A*B + (A*B|subject) + (1|item)): Model formula with interaction.

Onset

• overlap = (0.5,0.2),
• onset = UniformOnset(; offset = signalsize * overlap[1], width = signalsize * overlap[2]), # Put offset to 1 for no overlap. put width to 0 for no jitter

Noise

• noiselevel = 0.2.
• noise = PinkNoise(; noiselevel = noiselevel).

Other parameters

• return_epoched = false: If true, already epoched data is returned. Otherwise, continuous data is returned.

Returns

• (data, events)::Tuple{Array, DataFrame}:

  • data contains the simulated EEG data. Its dimensionality depends on the status of return_epoched and whether nsubjects >=2: `continuoustime,continuoustime x nsubjects,times x size(design)[1]ortimes x size(design)[1] x n_subjects`.

  • events contains the event combinations based on the experimental design. Each row corresponds to one combination of condition levels which is often equivalent to one stimulus or trial.

Examples

julia> data, events = UnfoldSim.predef_2x2();

julia> events
100×3 DataFrame
 Row │ A        B        latency 
     │ String   String   Int64   
─────┼───────────────────────────
   1 │ a_small  b_large       62
   2 │ a_big    b_tiny       132
  ⋮  │    ⋮        ⋮        ⋮
  99 │ a_big    b_large     5883
 100 │ a_big    b_tiny      5935
                  96 rows omitted

julia> data
6035-element Vector{Float64}:
  0.32384561187165956
  0.4108799249488322
  0.4715514814540277
  0.5936253009785152
  ⋮
  0.047664408295942984
  0.051193074728432035
 -0.08951617593287822
 -0.14456000097460356

See also predef_eeg.

predef_eeg(; <keyword arguments>)
predef_eeg(rng; <keyword arguments>)
predef_eeg(rng, n_subjects; <keyword arguments>)

Simulate data for a P1/N1/P3 component complex.

Note that this function is mainly used for demonstration and internal testing purposes or whenever a quick simulation is needed.

In case n_subjects is defined - MixedModelComponents are generated (multi-subject simulation), else LinearModelComponents (single-subject simulation).

The most used keyword argument is: return_epoched = true which returns already epoched data. If you want epoched data without overlap, specify onset = NoOnset() and return_epoched = true.

Keyword arguments

Design

• n_repeats = 100: Number of times the experimental design is repeated. Only used in the single-subject case.
• n_items = 100: Number of items. Only used in the multi-subject case.
• event_order_function = shuffle: Random trial order. Use event_order_function = (rng, x) -> x to deactivate.
• conditions = Dict(:condition => ["car", "face"], :continuous => range(-5, 5, length = 10)): Conditions and covariates used in this predefined design.

Component / Signal

• sfreq = 100: Sampling frequency.
• p1 = (p100(; sfreq = sfreq), @formula(0 ~ 1), [5], Dict()): P1 with amplitude 5; no effects.
• n1 = (n170(; sfreq = sfreq), @formula(0 ~ 1 + condition), [5, 3], Dict()): N1 with amplitude 5, dummy-coded condition effect (levels "car", "face") of 3.
• p3 = (p300(; sfreq = sfreq), @formula(0 ~ 1 + continuous), [5, 1], Dict()): P3 with amplitude 5, continuous effect range [-5,5] with slope 1.

Onset

• overlap = (0.5,0.2), # convenient parameterization for the default onset::UniformOnset. (offset, width) in seconds. If you do not want any overlap, either use onset=NoOnset(), or put the offset to a value larger than the maximum used component length, e.g. overlap=(1,0.2). Put the width to 0 to have no jitter between events.
• onset = UniformOnset(; offset = sfreq * 0.5 * overlap[1], width = sfreq * 0.5 * overlap[2]),

Noise

• noiselevel = 0.2.
• noise = PinkNoise(; noiselevel = noiselevel).

Other parameters

• return_epoched = false: If true, already epoched data is returned. Otherwise, continuous data is returned.

Returns

• (data, events)::Tuple{Array, DataFrame}:

  • data contains the simulated EEG data. Its dimensionality depends on the status of return_epoched and whether it's a single- or multisubject simulation (1D, 2D or 3D array): continuous_time, continuous_time x n_subjects, times x size(design)[1] or times x size(design)[1] x n_subjects.

  • events contains the event combinations based on the experimental design. Each row corresponds to one combination of condition/covariate levels which is often equivalent to one stimulus or trial.

Examples

julia> data, events = UnfoldSim.predef_eeg();

julia> events
2000×3 DataFrame
  Row │ continuous  condition  latency 
      │ Float64     String     Int64   
──────┼────────────────────────────────
    1 │   2.77778   car             62
    2 │  -5.0       face           132
  ⋮   │     ⋮           ⋮         ⋮
 1999 │  -0.555556  face        120096
 2000 │  -2.77778   car         120154

julia> data
120199-element Vector{Float64}:
  0.31631798033146774
  0.40338935529989906
  0.46409775558165056
  0.5862082040156747
  ⋮
 -0.1879589005111152
 -0.3163314509311509
 -0.22230944464885682
 -0.01320095208877194

See also predef_2x2.

simulate(
rng::AbstractRNG,
design::AbstractDesign,
components,
onset::AbstractOnset,
noise::AbstractNoise = NoNoise();
return_epoched = false,
)

simulate(
design::AbstractDesign,
components,
onset::AbstractOnset,
noise::AbstractNoise = NoNoise();
return_epoched = false,
)

Simulate continuous or epoched signal, given design, [Array of] component, onset and optional noise and rng. Main simulation function.

Arguments

• rng::AbstractRNG (optional): Random number generator, important to ensure reproducibility.
• design::AbstractDesign: Desired experimental design.
• components: Component(s) for the desired signal.
• onset::AbstractOnset: Desired inter-onset distance distribution.
• noise::AbstractNoise = NoNoise() (optional): Desired noise.

Keyword arguments

• return_epoched::Bool = false: If set to true epoched data is returned, otherwise a continuous signal is returned (see also Notes below).

Returns

• (signal, events)::Tuple{Array, DataFrame}:
  • signal : Generated signal. Depending on the design, on the components and on return_epoched, the output can be a 1-D, 2-D, 3-D or 4-D Array. For example, a 4-D Array would have the dimensions channels x time x trials x subjects.
  • events: Generated events data frame with latencies.

Examples

Adapted from the quickstart tutorial in the UnfoldSim docs.

julia> using Random # to get an RNG

julia> design =
    SingleSubjectDesign(; conditions = Dict(:cond_A => ["level_A", "level_B"])) |>
    x -> RepeatDesign(x, 10);

julia> component = LinearModelComponent(; 
    basis = [0, 0, 0, 0.5, 1, 1, 0.5, 0, 0],
    formula = @formula(0 ~ 1 + cond_A),
    β = [1, 0.5],
);

julia> onset = UniformOnset(; width = 20, offset = 4);

julia> noise = PinkNoise(; noiselevel = 0.2);

# Variant 1: Use a custom RNG.
julia> data, events = simulate(MersenneTwister(2), design, component, onset, noise);

julia> data
293-element Vector{Float64}:
 -0.013583193323430123
  0.09159433856866195
  ⋮
 -0.25190584567097907
 -0.20179992275876316

julia> events
20×2 DataFrame
 Row │ cond_A   latency 
     │ String   Int64   
─────┼──────────────────
   1 │ level_A        9
   2 │ level_B       20
   3 │ level_A       27
   4 │ level_B       37
  ⋮  │    ⋮        ⋮
  18 │ level_B      257
  19 │ level_A      271
  20 │ level_B      284
         13 rows omitted

# Variant 2: Without specifying an RNG, MersenneTwister(1) will be used for the simulation.
julia> data1, events1 = simulate(design, component, onset, noise);
┌ Warning: No random generator defined, used the default (`Random.MersenneTwister(1)`) with a fixed seed. This will always return the same results and the user is strongly encouraged to provide their own random generator!

Notes

Some remarks on how the noise is added:

• If return_epoched = true and onset = NoOnset() the noise is added to the epoched data matrix.
• If onset is not NoOnset, a continuous signal is created and the noise is added to this i.e. this means that the noise won't be the same as in the onset = NoOnset() case even if return_epoched = true.
• The case return_epoched = false and onset = NoOnset() is not possible and therefore covered by an assert statement.

Additional remarks on the overlap of adjacent signals when return_epoched = true:

• If onset = NoOnset() there will not be any overlapping signals in the data because the onset calculation and conversion to a continuous signal is skipped.
• If an inter-onset distance distribution is given, a continuous signal(potentially with overlap) is constructed and partitioned into epochs afterwards.

simulate_and_add!(epoch_data::AbstractMatrix, c, simulation, rng)
simulate_and_add!(epoch_data::AbstractArray, c, simulation, rng)

Helper function to call simulate_component and add it to a provided Array epoch_data.

simulate_component(rng, c::AbstractComponent, simulation::Simulation)

By default call simulate_component with (rng, c::Abstractcomponent, design::AbstractDesign) instead of the whole simulation. This function exist solely to provide a "hook" if for a custom component something else than the design is necessary, e.g. a dependency on the onsets, noise or similar.

simulate_component(rng, c::LinearModelComponent, design::AbstractDesign)

Generate a linear model design matrix, weight it by the coefficients c.β and multiply the result with the given basis vector.

Returns

• Matrix{Float64}: Simulated component for each event in the events data frame. The output dimensions are length(c.basis) x length(design).

Examples

julia> design = SingleSubjectDesign(; conditions = Dict(:cond => ["natural", "artificial"]));

julia> c = UnfoldSim.LinearModelComponent([0, 1, 1, 0], @formula(0 ~ 1 + cond), [1, 2], Dict());

julia> using StableRNGs

julia> simulate_component(StableRNG(1), c, design)
4×2 Matrix{Float64}:
 0.0  0.0
 3.0  1.0
 3.0  1.0
 0.0  0.0

simulate_component(rng, c::MixedModelComponent, design::AbstractDesign, return_parameters = false)

Generate a MixedModel and simulate data according to the given parameters c.β and c.σs.

Keyword arguments

• return_parameters::Bool = false: Can be used to return the per-event parameters used to weight the basis function. Sometimes useful to inspect what is simulated.

Returns

• Matrix{Float64}: Simulated component for each event in the events data frame. The output dimensions are length(c.basis) x length(design).

Notes

1. MixedModels/Sim does not allow simulation of data without white noise of the residuals. Because we want our own noise, we use the following trick to remove the MixedModels-Noise:

Practically, we upscale the specified σs by factor 10000, and request a white-noise-level of σ = 0.0001. Internally in MixedModels/Sim, σs are relative to σ, and thus are normalized correctly, while keeping the noise 10000 times smaller than the random effects

We cannot exclude that this trick runs into strange numerical issues if the random effect σs are very large compared to the fixed effects.

1. Currently, it is not possible to use a different basis for fixed and random effects. If this is needed, some code-scaffold is available but commented out at the moment and requires a bit of implementation work.

Examples

julia> design = MultiSubjectDesign(; n_subjects = 2, n_items = 6, items_between = Dict(:cond => ["A", "B"]));

julia> c = UnfoldSim.MixedModelComponent([0, 1, 1, 0], @formula(0 ~ 1 + cond + (1|subject)), [1, 2], Dict(:subject => [2],), Dict());

julia> using StableRNGs

julia> simulate_component(StableRNG(1), c, design)
4×12 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
 -2.70645  -0.706388  -2.70632  -0.706482  -2.7066  -0.706424  0.325722  2.32569  0.325627  2.32564  0.325468  2.32565
 -2.70645  -0.706388  -2.70632  -0.706482  -2.7066  -0.706424  0.325722  2.32569  0.325627  2.32564  0.325468  2.32565
 -0.0      -0.0       -0.0      -0.0       -0.0     -0.0       0.0       0.0      0.0       0.0      0.0       0.0

julia> simulate_component(StableRNG(1), c, design, return_parameters = true)
1×12 Matrix{Float64}:
 -2.70645  -0.706388  -2.70632  -0.706482  -2.7066  -0.706424  0.325722  2.32569  0.325627  2.32564  0.325468  2.32565

simulate_component(rng, c::MultichannelComponent, design::AbstractDesign)

Return the projection of a MultichannelComponent c from "source" to "sensor" space.

Returns

• Array{Float64,3}: Projected simulated component for each event in the events data frame. The output dimensions are length(c.projection) x length(c.basis) x length(design).

Examples

julia> design = SingleSubjectDesign(; conditions = Dict(:cond => ["natural", "artificial"]));

julia> c = LinearModelComponent(; basis = p100(), formula = @formula(0 ~ 1 + cond), β = [1, 0.5]);
julia> mc = UnfoldSim.MultichannelComponent(c, [1, 2, -1, 3, 5, 2.3, 1], PinkNoise());

julia> using StableRNGs

julia> simulate_component(StableRNG(1), mc, design)
7×15×2 Array{Float64, 3}:
[:, :, 1] =
  0.859626   1.16574   0.959523   0.757522   1.17639   1.65156   1.3742    1.76706   2.76971   2.0306    1.17429   1.00922   1.09519   0.754659   2.25662
  1.71925    2.33147   1.91905    1.51504    2.35278   3.30312   2.7484    3.53412   5.53942   4.0612    2.34858   2.01845   2.19039   1.50932    4.51324
 -0.859626  -1.16574  -0.959523  -0.757522  -1.17639  -1.65156  -1.3742   -1.76706  -2.76971  -2.0306   -1.17429  -1.00922  -1.09519  -0.754659  -2.25662
  2.57888    3.49721   2.87857    2.27257    3.52917   4.95469   4.1226    5.30118   8.30913   6.09179   3.52287   3.02767   3.28558   2.26398    6.76985
  4.29813    5.82868   4.79761    3.78761    5.88194   8.25781   6.871     8.8353   13.8485   10.153     5.87145   5.04612   5.47597   3.7733    11.2831
  1.97714    2.68119   2.2069     1.7423     2.70569   3.79859   3.16066   4.06424   6.37033   4.67037   2.70087   2.32121   2.51894   1.73572    5.19022
  0.859626   1.16574   0.959523   0.757522   1.17639   1.65156   1.3742    1.76706   2.76971   2.0306    1.17429   1.00922   1.09519   0.754659   2.25662

[:, :, 2] =
  0.859626   1.16574   0.959523   0.757522   1.17639   1.65156   1.31571   1.56047   2.39471   1.54567   0.689367   0.634223   0.888605   0.69617   2.25662
  1.71925    2.33147   1.91905    1.51504    2.35278   3.30312   2.63142   3.12094   4.78942   3.09135   1.37873    1.26845    1.77721    1.39234   4.51324
 -0.859626  -1.16574  -0.959523  -0.757522  -1.17639  -1.65156  -1.31571  -1.56047  -2.39471  -1.54567  -0.689367  -0.634223  -0.888605  -0.69617  -2.25662
  2.57888    3.49721   2.87857    2.27257    3.52917   4.95469   3.94713   4.68142   7.18413   4.63702   2.0681     1.90267    2.66582    2.08851   6.76985
  4.29813    5.82868   4.79761    3.78761    5.88194   8.25781   6.57855   7.80236  11.9735    7.72837   3.44684    3.17112    4.44303    3.48085  11.2831
  1.97714    2.68119   2.2069     1.7423     2.70569   3.79859   3.02613   3.58909   5.50783   3.55505   1.58554    1.45871    2.04379    1.60119   5.19022
  0.859626   1.16574   0.959523   0.757522   1.17639   1.65156   1.31571   1.56047   2.39471   1.54567   0.689367   0.634223   0.888605   0.69617   2.25662

simulate_interonset_distances(rng, onset::AbstractOnset, design::AbstractDesign)

Generate the inter-onset distance vector by sampling from the respective distribution (in samples).

Arguments

• rng: Random number generator (RNG) to make the process reproducible.
• onset::AbstractOnset: Inter-onset distance distribution to sample from.
• design::AbstractDesign: Experimental design with conditions and covariates.

Returns

• Vector{Integer}: Inter-onset distances in samples. Note that these are distances between onsets and no latencies.

Examples

# Create an experimental design
julia> design_single = SingleSubjectDesign(;
           conditions = Dict(
               :stimulus_type => ["natural", "artificial"],
               :contrast_level => range(0, 1, length = 3),
           ),
       );

# Create an inter-onset distance distribution
julia> onset_distribution = LogNormalOnset(3, 0.5, 5, nothing);

julia> using StableRNGs

julia> simulate_interonset_distances(StableRNG(1), onset_distribution, design_single)
6-element Vector{Int64}:
 20
 26
 34
 18
 12
 23

See also simulate_onsets.

simulate_noise(rng, t::AbstractNoise, n::Int)

Generate noise samples of the given type t.

For details, see the documentation of the individual noise types. Use subtypes(AbstractNoise) for a list of the implemented noise types.

Arguments

• rng::AbstractRNG: Random number generator (RNG) to make the process reproducible.
• t::AbstractNoise: Instance of a noise type e.g. PinkNoise().
• n::Int: The number of noise samples that should be generated.

Returns

• Vector: Vector of length n containing the noise samples.

Examples

# Here we use White Noise as an example but it works in the same way for the other noise types.
julia> noise = WhiteNoise()
WhiteNoise
  noiselevel: Int64 1
  imfilter: Int64 0

julia> using StableRNGs

julia> simulate_noise(StableRNG(1), noise, 3)
3-element Vector{Float64}:
 -0.5325200748641231
  0.098465514284785
  0.7528865221245234

simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation)

Call simulate_interonset_distances to generate distances between events and then add them up to generate the actual latencies in samples.

Please note that this function is mainly for internal use in the context of simulate function calls.

Also note that the accumulation of onsets starts at 1 to avoid indexing problems in the case that the first sampled onset is 0.

Arguments

• rng: Random number generator (RNG) to make the process reproducible.
• onset::AbstractOnset: Inter-onset distance distribution which is passed to simulate_interonset_distances.
• simulation::Simulation: Simulation object which contains design, component(s), inter-onset distance distribution and noise.

Returns

Examples

# Create Simulation object
julia> design_single = SingleSubjectDesign(;
           conditions = Dict(
               :stimulus_type => ["natural", "artificial"],
               :contrast_level => range(0, 1, length = 3),
           ),
       );

julia> p1_component = LinearModelComponent(; basis = p100(), formula = @formula(0 ~ 1), β = [5]);

julia> simulation = Simulation(design_single, p1_component, UniformOnset(), NoNoise());

julia> using StableRNGs

# Simulate onsets for this simulation
julia> simulate_onsets(StableRNG(1), simulation.onset, simulation)
6-element Vector{Int64}:
  20
  70
  97
 110
 150
 182

See also simulate_interonset_distances.

simulate_responses(
    rng,
    components::Vector{<:AbstractComponent},
    simulation::Simulation)

Simulate multiple component responses and accumulate them on a per-event basis.

Returns

• epoch_data: Matrix (or Array in the multi-channel case) of combined simulated components. The output dimensions are maxlength(components) x length(simulation.design) for single-channel components and n_channels(components) x maxlength(components) x length(simulation.design) for multi-channel components.

Examples

julia> design = SingleSubjectDesign(; conditions = Dict(:cond => ["natural", "artificial"]));

julia> c1 = LinearModelComponent(; basis = p100(), formula = @formula(0 ~ 1 + cond), β = [1, 0.5]);
julia> c2 = LinearModelComponent(; basis = p300(), formula = @formula(0 ~ 1), β = [2]);

julia> simulation = Simulation(design, [c1, c2], UniformOnset(; width = 0, offset = 30), PinkNoise());

julia> using StableRNGs

julia> simulate_responses(StableRNG(1), [c1, c2], simulation)
45×2 Matrix{Float64}:
 0.0        0.0
 0.0        0.0
 0.0        0.0
 0.0        0.0
 0.0        0.0
 ⋮          
 0.352614   0.352614
 0.203907   0.203907
 0.0924246  0.0924246
 0.0233794  0.0233794
 0.0        0.0

weight_σs(σs::Dict, b_σs::Float64, σ_lmm::Float64)

Weight a σs Dict for MixedModels.jl by b_σs, a scaling factor typically from a basis.

Finally scales it again by σ_lmm, as a trick to simulate noise-free LMMs (see MixedModelsComponent)

In the future, we anticipate a function function weight_σs(σs::Dict,b_σs::Dict,σ_lmm::Float64) where each σs entry can be weighted individually by a matching b_σs, but it is not implemented.

Arguments

• σs::Dict = a Dict of random effects as output of MixedModels.create_re
• b_σs::Float64 = a scaling factor, typically one entry of a basis function from a component
• σ_lmm::Float64 = a scaling factor to simulate near-zero noise LMMs

Returns

`NamedTuple` of the weighted random effects