effects(design::AbstractDict, model::UnfoldModel; typical = mean)

Calculates marginal effects for all term combinations in design.

Implementation based on Effects.jl package; likely could repackage in UnfoldEffects.jl; somebody wants to do it? This would make it easier to cross-maintain it to changes/bug fixes in the Effects.jl package. design is a dictionary containing those predictors (as keys) with levels (as values), that you want to evaluate. The typical refers to the value, which other predictors that are not specified in the dictionary, should take on.

For MixedModels, the returned effects are based on the "typical" subject, i.e. all random effects are put to 0.

Example

 julia> f = @formula 0 ~ categoricalA + continuousA + continuousB
 julia> uf = fit(UnfoldModel, (Any => (f, times)), data, events)
 julia> d = Dict(:categorical => ["levelA", "levelB"], :continuous => [-2, 0, 2])
 julia> effects(d, uf)

will result in 6 predicted values: A/-2, A/0, A/2, B/-2, B/0, B/2.

FileIO.load(file, ::Type{<:UnfoldModel}; generate_Xs=true)

Load UnfoldModel from a .jld2 file.

By default, the designmatrix is reconstructed. If it is not needed set generate_Xs=false which improves time-efficiency.

FileIO.save(file, uf::T; compress=false) where {T<:UnfoldModel}

Save UnfoldModel in a (by default uncompressed) .jld2 file.

For memory efficiency the designmatrix is set to missing. If needed, it can be reconstructed when loading the model.

coefnames(term)

coefnames of a TimeExpandedTerm concatenates the basis-function name with the kronecker product of the term name and the basis-function colnames. Separator is ' : ' Some examples for a firbasis: basis313 : (Intercept) : 0.1 basis313 : (Intercept) : 0.2 basis_313 : (Intercept) : 0.3 ...

fit(type::UnfoldModel,d::Vector{Pair},tbl::AbstractDataFrame,data::Array)
fit(type::UnfoldModel,f::FormulaTerm,tbl::AbstractDataFrame,data::Array{T,3},times)
fit(type::UnfoldModel,f::FormulaTerm,tbl::AbstractDataFrame,data::Array{T,2},basisfunction::BasisFunction)

Generates Designmatrix & fits model, either mass-univariate (one model per epoched-timepoint) or time-expanded (modeling linear overlap).

keyword arguments

• fit::Bool (default: true) - fit the model after constructing the designmatrix. Setting this to false is sometimes helpful if you only want to inspect the designmatrix.
• contrasts::Dict: (default: Dict()) contrast to be applied to formula. Example: Dict(:my_condition=>EffectsCoding()). More information here: https://juliastats.org/StatsModels.jl/stable/contrasts/
• eventcolumn::Union{Symbol,String} (default :event) - the column in tbl to differentiate the basisfunctions as defined in d::Vector{Pair}
• solver::function: (default: solver_default). The solver used for y=Xb, e.g. (X,y;kwargs...) -> solver_default(X,y;kwargs...). There are faster & alternative solvers available, see solver_predefined for a list of options, see solver benchmark in the online documentation. To use the GPU, you can provide the data as a CuArray after using CUDA. Please change the solver to e.g. solver_predef(X,y;solver=:qr) as lsmr+cuda => crash typically. It's worth though, speed increases >100x possible
• show_progress::Bool (default true) - show progress via ProgressMeter - passed to solver
• eventfields::Array: (optional, default[:latency]) Array of symbols, representing column names intbl`, which are passed to basisfunction event-wise. First field of array always defines eventonset in samples.

If a Vector[Pairs] is provided, it has to have one of the following structures: For deconvolution analyses (use Any=>(f,bf) to match all rows of tbl in one basis functions). Assumes data is a continuous EEG stream, either a Vector or a ch x time Matrix

f1 = @formula(0~1+my_condition)
[
 :A=>(f1,firbasis((-0.1,1),128), # sfreq = 128Hz
 :B=>(f2,firbasis((-3,2),128)
]

for mass-univariate analyses without deconvolution. Assumes data to be cut into epochs already (see Unfold.epoch). Follows eeglab standard ch x time x trials:

timesvector = range(-0.1,3,step=1/100)
[
 :A=>(f1,timesvector),
 :B=>(f2,timesvector)
]

Notes

• The type can be specified directly as well e.g. fit(type::UnfoldLinearModel) instead of relying on the automatic inference
• The data is reshaped if it is missing one dimension to have the first dimension then 1 "Channel".

Examples

Mass Univariate Linear

julia> data,evts = UnfoldSim.predef_eeg()
julia> data_e,times = Unfold.epoch(data=data,tbl=evts,τ=(-1.,1.9),sfreq=100) # cut the data into epochs. data_e is now ch x times x epoch

julia> f  = @formula 0~1+continuousA+continuousB
julia> model = fit(UnfoldModel,f,evts,data_e,times)
# or:
julia> model = fit(UnfoldModel,[Any=>(f,times)],evts,data_e)

Timexpanded Univariate Linear

julia> basisfunction = firbasis(τ=(-1,1),sfreq=10)
julia> model = fit(UnfoldModel,f,evts,data,basisfunction)
# or
julia> model = fit(UnfoldModel,[Any=>(f,basisfunction],evts,data)

StatsModels.modelmatrix(uf::UnfoldLinearModelContinuousTime, basisfunction = true)

Setting the optional second args to false, will return the modelmatrix without the timeexpansion / basisfunction applied.

modelmatrix(uf::UnfoldLinearModel)

returns the modelmatrix of the model. Concatenates them, except in the MassUnivariate cases, where a vector of modelmatrices is return

Compare with modelmatrices which returns a vector of modelmatrices, one per event

function predict(
    uf::UnfoldModel,
    f::Vector{<:FormulaTerm},
    evts::Vector{<:DataFrame};
    overlap::Bool = true,
    kwargs...
)

Returns a predicted ("y_hat = X*b") Array.

• uf is an <:UnfoldModel
• f is a (vector of) formulas, typically Unfold.formulas(uf), but formulas can be modified e.g. by effects.
• evts is a (vector of) events, can be Unfold.events(uf) to return the (possibly continuous-time) predictions of the model. Can be a custom even

kwargs:

if overlap = true (default), overlap based on the latency column of evts will be simulated, or in the case of !ContinuousTimeTrait just X*coef is returned.

if overlap = false, returns predictions without overlap (models with ContinuousTimeTrait (=> with basisfunction / deconvolution) only), via predict_no_overlap

if keep_basis or exclude_basis is defined, then predict_partial_overlap is called, which allows to selective introduce overlap based on specified (or excluded respective) events/basisfunctions

epoch_to and epoch_timewindow: calculate (partial) overlap controlled predictions, but returns them at the specified epoch_at event, with the times epoch_timewindow (default is taken from the basisfunction) in samples.

eventcolumn can be specified as well if different from the default event.

Hint: all kwargs can be Vector, or if e.g. string types are provided, will be put into a length==1 vector.

Output

• If overlap=false, returns a 3D-Array
• If overlap=true and epoch_to = nothing (default), returns a 2D-array
• If overlap=true and epoch_to != nothing, returns a 3D array

modelcols(term, tbl)

calculates the actual designmatrix for a timeexpandedterm. Multiple dispatch on StatsModels.modelcols

_modelcols(form::FormulaTerm, events)

_modelcols(forms::Vector,events::Vector)

A wrapper around StatsModels.modelcols that is only needed for easy multiple dispatch

apply_basisfunction(
    form,
    basisfunction,
    eventfields,
    eventname
)

timeexpand the rhs-term of the formula with the basisfunction

combine_yhat(list,single)

combines single into list, if either list or single contains missing, automatically casts the respective counter-part to allow missings as well

!!! Important: this is an ugly hack dating back to the time where UnfoldMixedModels was still an extension. We are overloading this function in UnfoldMixedModels.jl with a more specific type, to switch between MixedModels-Unfold types and not...

designmatrix(type, f, tbl; kwargs...)

call without basis function, continue with basisfunction = nothing

designmatrix(
    unfoldmodeltype,
    f,
    tbl,
    basisfunction;
    contrasts,
    eventname,
    kwargs...
)

designmatrix(type, f, tbl; kwargs...) Return a DesignMatrix used to fit the models.

Arguments

• type::UnfoldModel
• f::FormulaTerm: Formula to be used in this designmatrix
• tbl: Events (usually a data frame) to be modelled
• basisfunction::BasisFunction: basisfunction to be used in modeling (if specified)
• contrasts::Dict: (optional) contrast to be applied to formula
• eventfields::Array: (optional) Array of symbols which are passed to basisfunction event-wise.

First field of array always defines eventonset in samples. Default is [:latency]

Examples

julia>  designmatrix(UnfoldLinearModelContinuousTime,Dict(Any=>(f,basisfunction1),tbl)

designmatrix(
    T::Type{<:UnfoldModel},
    design_array::Vector{<:Pair},
    tbl;
    eventcolumn = :event,
    contrasts = Dict{Symbol,Any}(),
    kwargs...,

iteratively calls designmatrix for each event in the design_array, and returns a list of <:AbstractDesignMatrix

designmatrix(
    uf::UnfoldModel,
    tbl;
    eventcolumn = :event,
    contrasts = Dict{Symbol,Any}(),
    kwargs...,

Main function called from fit(UnfoldModel...), generates the designmatrix, returns a list of <:AbstractDesignMatrix

[X,y] = drop_missing_epochs(X, y::Array)

Helper function to remove epochs of y that contain missings. Drops them from both X and y. Often used in combination with Unfold.epoch

X can be anything that has two dimensions (Matrix, DataFrame etc)

empty_modelmatrix(d::AbstractDesignMatrix)

returns an empty modelmatrix of the type DesignMatrix type of d

epoch(data::Array{T,1},evts::DataFrame,τ::Tuple/Vector,sfreq;kwargs...,

Basic function to epoch data; all input also available as kwargs.

Additional kwarg: eventtime=:latency, which defines the column in evts that is used to cut the data (in samples). For uneven sample-times we use round()`

equalize_size(X, data)

Equates the length of data and designmatrix by cutting the shorter one

The reason we need this is because when generating the designmatrix, we do not know how long the data actually are. We only assume that event-latencies are synchronized with the data

firbasis(τ, sfreq; ...)
firbasis(τ, sfreq, name; interpolate, scale_duration)

Generate a sparse FIR basis around the τ timevector at sampling rate sfreq. This is useful if you cannot make any assumptions on the shape of the event responses. If unrounded events are supplied, they are split between samples. E.g. event-latency = 1.2 will result in a "0.8" and a "0.2" entry.

Advanced: second input can be duration in samples - careful: times(firbasis) always assumes duration = 1. Therefore, issues with LMM and predict will appear!

keyword arguments

interpolate (Bool, default false): if true, interpolates events between samples linearly. This results in predict functions to return a trailling 0scale_duration(Union{Bool,Interpolations-Interpolator}, default false): if true, scales the response by the fit-kwargseventfieldssecond entry. That is, the FIR becomes a stepfunction instead of a impulse response. if Interpolations.interpolator, e.g.Interpolations.Linear()- uses the fit-kwargseventfieldssecond entry to stretch the FIR kernel based onimresize`. This implements Hassall

Examples

Generate a FIR basis function from -0.1s to 0.3s at 100Hz

julia>  f = firbasis([-0.1,0.3],100)

Evaluate at an event occuring at sample 103.3

julia>  f(103.3)

firkernel(ev, times; interpolate, scale_duration)

Calculate a sparse firbasis

second input can be duration in samples - careful: times(firbasis) always assumes duration = 1. Therefore, issues with LMM and predict will appear!

Examples

julia>  f = firkernel(103.3,range(-0.1,step=0.01,stop=0.31))
julia>  f_dur = firkernel([103.3 4],range(-0.1,step=0.01,stop=0.31))

formulas(design::Vector{<:Pair})

returns vector of formulas, no schema has been applied (those formulas never saw the data). Also no timeexpansion has been applied (in the case of timecontinuous models)

get_basis_colnames(m)
get_basis_colnames(formulas)

returns list of colnames - e.g. times for firbasis.

get_basis_indices(uf, basisnames::Vector)

returns a boolean vector with length spanning all coefficients, which coefficient is defined by basisnames (vector of names)

get_basisnames(model::UnfoldModel)

Return the basisnames for all predictor terms as a vector.

The returned vector contains the name of the event type/basis, repeated by their actual coefficient number (after StatsModels.apply_schema / timeexpansion). If a model has more than one event type (e.g. stimulus and fixation), the vectors are concatenated.

hrfbasis(TR; parameters, name)

Generate a Hemodynamic-Response-Functio (HRF) basis with inverse-samplingrate "TR" (=1/FS)

Optional Parameters p: defaults {seconds} p(1) - delay of response (relative to onset) 6 p(2) - delay of undershoot (relative to onset) 16 p(3) - dispersion of response 1 p(4) - dispersion of undershoot 1 p(5) - ratio of response to undershoot 6 p(6) - onset {seconds} 0 p(7) - length of kernel {seconds} 32

Examples

Generate a HRF basis function object with Sampling rate 1/TR. And evaluate it at an event occuring at TR 103.3 with duration of 4.1 TRs

julia>  f = hrfbasis(2.3)
julia>  f(103.3,4.1)

hrfkernel(e, TR, p)

Calculate a HRF kernel. Input e can be [onset duration]

Examples

julia>  f = hrfkernel(103.3,2.3,[6. 16. 1. 1. 6. 0. 32.])

linearize(x)

Flatten a 1D array from of a 2D/3D array. Also drops the empty dimension

Returns a view of the Matrix y, according to the indices of the timeexpanded basisname

modelmatrices(X::AbstractDesignMatrix)
modelmatrices(X::Vector{<:AbstractDesignMatrix})
modelmatrices(modelmatrix::AbstractMatrix)

Returns the modelmatrices (also called designmatrices) separately for the events. This is similar to StatsModels.modelcols, but merely access the precomputed designmatrix. If the designmatrix needs to be computed, please use modelcols

Compare to modelmatrix which further concatenates the designmatrices (in the ContinuousTime case).

predict_no_overlap(, uf, coefs, f, evts)

in ContinuousTime case (typically the deconvolution model), we return idealized predictions without overlap between events.

in the Not-ContinuousTime case (typically the MassUnivariate model), we return predictions for each event independently. In that case, the function is unfortunately a missnomer, as overlap cannot be removed from mass-univariate models.

predict_partial_overlap(, uf, args; kwargs...)

Returns predicted time-continuous values, but only for a subset of events. This is achieved by excluding the part of the designmatrix that belongs to the basisfunctions/events you do not want to have in your model.

Typically called via predict, for configuration, keyword-arguments and usage see there.

One difference is, that we require the coefs(uf::UnfoldModel) already exctracted.

Due to the time-continuous nature, running it with a model not containing the ContinuousTimeTrait it will throw an error.

predicttable(model<:UnfoldModel,events=Unfold.events(model),args...;kwargs...)

Shortcut to call efficiently call (pseudocode) result_to_table(predict(...)).

Returns a tidy DataFrame with the predicted results. Loops all input to predict, but really only makes sense to use if you specify either:

overlap = false (the default) or epoch_to = "eventname".

prepare(X, data)

convert and permutedim input to follow the following output:

Ĥ, Y, X = prepare(X, data)

where Ĥ is used to save the beta, Y is the data in format ch x repeat x time (with size(time) = 1 if data is a Matrix), and X.

• if data is a CuArray, everything is transformed to CuArrays as well (via UnfoldCUDAExt.jl, CUDA needs to be loaded)
• same datatype between X and data is enforced

prepare_XTX(all)

instead of solving y = Xb, we solve X'Xb = X'y. This function calculates X'X and instantiates X'y to be used in the solver-step, to facilitate X'y calculations later, X' is also calculated.

prepare_pinv(all)

calculates pinv of the designmatrix for later use in the solver-step. This is helpful in case you have many chanels

result_to_table(model<:UnfoldModel, eff::AbstractArray, events::Vector{<:DataFrame})
result_to_table(
    eff::AbstractArray,
    events::Vector{<:DataFrame},
    times::Vector{<:Vector{<:Number}},
    eventnames::Vector)
result_to_table(
    eff::Vector{<:AbstractArray},
    events::Vector{<:DataFrame},
    times::Vector,
    eventnames::Vector,
)

Converts an array-result (prediction or coefficient) together with the events, to a tidy dataframe.

To support multi-event models, we expect everything to be put into Vectors - this should be refactored at some point to be compatible with broadcasting, but it is not right now.

args

eff: A vector that contains the array(s) to be converted to a tidy dataframe. Each event in your events dataframe your have it's own array (i.e. the array should have length(eff) == length(unique(events.event))). Each array should be 3D, with channel x time x predictor events: A vector of event-dataframes, each need to match size(eff,3). Each individual event (unique(events.event)) should have it's own dataframe. E.g. [[df_event1::DataFrame], [df_event1::DataFrame]]times: A vector of time-vectors withlength(eff), each time-vector needs to matchsize.(eff,2)eventnames`: A vector of eventnames, either symbols or strings, should be a single entry per event

solver_default(X, y; kwargs...)

default solvers.

• If data is continuous (2D), we solve Xb = y via lsmr
• If data is epoched (3D) we solve Xb = y via pinv

We highly recommend to check out solver_predefined for faster options by rather solving X'Xb = X'y via QR, cholesky, pinv or ``-solver. A benchmark is available in the online documentation.

Please see ?solver_main for keyword arguments of the solver (like stderror, multithreading, show_time, show_progress)

solver_predefined(X, y_in; solver, kwargs...)

helper function that returns solver with appropriate prepare-pipelines and fitting solver-functions. X is a (typically sparse) designmatrix, y is a 2D or 3D array.

solver : one of :cg, :pinv, :intern, :qr, :cholesky, :lsmr (default)

Only lsmr solves Xb = y via an iterative solver and should be more accurate in principle.

The other predefined-solvers solve X'Xb = X'y which is often computationally much cheaper, and because X'X can be precalculated, it should be cheaper to apply.

Testing this empirically is somewhat complicated, as depending on your sparsity structure (≈ your design) and the size of your data (sfreq & minutes) the best solver and the reached accuracy can change quite a bit.

GPU

All solvers except :lsmr support GPU calculations. For lsmr on the GPU try solver_krylov instead

Speed improved version of spdiagm, takes a single float value instead of a vector, like a version of spdiagm that takes in a UniformScaling

e.g.

sz = 5 ix = [1,3,10] spdiagm_diag(sz,(.-ix.=>1)...)

Helper function to decide whether all bases have the same number of columns per event

timeexpand_cols(basisfunction, bases, ncolsBasis, ncolsX)

calculates in which rows the individual event-basisfunctions should go in Xdc

see also timeexpandrows timeexpandvals

timeexpand_rows(onsets, bases, shift, ncolsX)

calculates in which rows the individual event-basisfunctions should go in Xdc

timeexpandrows timeexpandvals

times(model<:UnfoldModel)

returns arrays of time-vectors, one for each basisfunction / parallel-fitted-model (MassUnivarite case)

wrapper to make apply_schema mixed models as extension possible

Note: type is not necessary here, but for LMM it is for multiple dispatch reasons!