diff --git a/.all-contributorsrc b/.all-contributorsrc index f452a97d..b197f4c5 100644 --- a/.all-contributorsrc +++ b/.all-contributorsrc @@ -116,6 +116,15 @@ "contributions": [ "doc" ] + }, + { + "login": "maxvanmigem", + "name": "Maximilien Van Migem", + "avatar_url": "https://avatars.githubusercontent.com/u/115144441?v=4", + "profile": "https://github.com/maxvanmigem", + "contributions": [ + "bug" + ] } ] } diff --git a/README.md b/README.md index acbec651..254895b5 100644 --- a/README.md +++ b/README.md @@ -143,6 +143,7 @@ You are very welcome to raise issues and start pull requests! suddha-bpn
suddha-bpn
🐛 Vladimir Mikheev
Vladimir Mikheev
🐛 📖 carmenamme
carmenamme
📖 + Maximilien Van Migem
Maximilien Van Migem
🐛 diff --git a/docs/literate/HowTo/contrasts.jl b/docs/literate/HowTo/contrasts.jl new file mode 100644 index 00000000..da8e1a17 --- /dev/null +++ b/docs/literate/HowTo/contrasts.jl @@ -0,0 +1,39 @@ +using CairoMakie +using Unfold +using UnfoldMakie +using UnfoldSim + + +## Contrast coding +# Unfold.jl uses the `StatsModels` package for the formula interface. This allows for a wide range of contrast coding schemes. For a full tutorial, please see [the StatsModels docs](https://juliastats.org/StatsModels.jl/stable/contrasts/). +# Please read their tutorial, as a motivation of why one would change the contrast coding scheme is outside of the realms of this package and more a basic linear regression question. + + +# !!! hint +# Given we have a nice `effects` implementation (mimicking `emmeans` and similar packages), coding schema is typically less important. + +# Here we will show a simple example of how to change the contrast coding scheme. We will use the `condition` variable, which has two levels, `A` and `B`. We will change the contrast coding from `Dummy` aka `Reference` aka 0/1 coding to `Sum` coding, which is the default in R. +eeg, evts = UnfoldSim.predef_eeg(noiselevel = 0) +f = @formula 0 ~ 1 + condition +basis = firbasis((-0.1, 0.6), 100) +m_dummy = fit(UnfoldModel, f, evts, eeg, basis) +m_effec = + fit(UnfoldModel, f, evts, eeg, basis; contrasts = Dict(:condition => EffectsCoding())) + + +# we could directly inspect the designmatrix +modelmatrix(m_dummy, false)[1][1:5, :] + +# and the effects coding +modelmatrix(m_effec, false)[1][1:5, :] + +# To confirm the difference in the actual fit, let's visualize them +c_d = coeftable(m_dummy) +c_e = coeftable(m_effec) +c_d.group .= "Dummy Coding" +c_e.group .= "Effects Coding" +c = vcat(c_d, c_e) + +plot_erp(c; mapping = (; color = :coefname, col = :group)) + +# As expected, the effects-coding slope of `condition: face` is half the size of the dummy-coding one (because -1/1 coding was used). diff --git a/docs/make.jl b/docs/make.jl index 5cb66a4a..a3002a4e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -40,6 +40,7 @@ makedocs( "Marginal effects (focus on non-linear predictors)" => "generated/HowTo/effects.md", #"Time domain basis functions"=>"generated/HowTo/timesplines.md", "Save and load Unfold models" => "generated/HowTo/unfold_io.md", + "Change contrasts / coding schema" => "generated/HowTo/contrasts.md", ], "Explanations" => [ "About basisfunctions" => "./explanations/basisfunctions.md", diff --git a/src/designmatrix.jl b/src/designmatrix.jl index d4fd2077..1a0a70ed 100644 --- a/src/designmatrix.jl +++ b/src/designmatrix.jl @@ -310,6 +310,12 @@ function extend_to_larger(modelmatrix1::SparseMatrixCSC, modelmatrix2::SparseMat end return hcat(modelmatrix1, modelmatrix2) end + + +""" + StatsModels.modelmatrix(uf::UnfoldLinearModelContinuousTime, basisfunction = true) +Setting the optional second args to false, will return the modelmatrix without the timeexpansion / basisfunction applied. +""" function StatsModels.modelmatrix(uf::UnfoldLinearModelContinuousTime, basisfunction = true) if basisfunction return modelmatrix(designmatrix(uf)) diff --git a/src/fit.jl b/src/fit.jl index 13df77eb..2767b91f 100644 --- a/src/fit.jl +++ b/src/fit.jl @@ -6,32 +6,52 @@ Generates Designmatrix & fits model, either mass-univariate (one model per epoched-timepoint) or time-expanded (modeling linear overlap). -- `eventcolumn` (Symbol/String, default :event) - the column in `tbl::AbstractDataFrame` to differentiate the basisfunctions as defined in `d::Vector{Pair}` -- `show_progress` (Bool, default true) - show Progress via ProgressMeter +## keyword arguments +- `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_defaut`). The solver used for `y=Xb`, e.g. `(X,y;kwargs...) -> solver_default(X,y;kwargs...)` +- `show_progress::Bool` (default `true`) - show progress via ProgressMeter - passed to `solver` +- `eventfields::Array: (optional, default `[:latency]`) Array of symbols, representing column names in `tbl`, 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: -`[:A=>(f,basisfunction), :B=>(f2,bf2)]` - for deconvolutin analyses (use `Any=>(f,bf)` to match all rows of `tbl` in one basis functins) -`[:A=>(f,timesvector), :B=>(f2,timesvector)]` - for mass univariate analyses. If multiple rERPs are calculated at the same time, the timesvectors must be the same - +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` +```julia +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`: +```julia +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 inferred -- The data is reshaped if it is missing one dimension to have the first dimesion then `1` "Channel". +- 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-repl -julia> data,evts = loadtestdata("testCase1") -julia> data_r = reshape(data,(1,:)) -julia> data_e,times = Unfold.epoch(data=data_r,tbl=evts,τ=(-1.,1.9),sfreq=10) # cut the data into epochs. data_e is now ch x times x epoch +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 # 1 +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-repl -julia> basisfunction = firbasis(τ=(-1,1),sfreq=10,name="A") -julia> model = fit(UnfoldModel,[Any=>(f,basisfunction],evts,data_r) +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) ``` """