better fit docs (#221)

* better fit docs

* docs: add maxvanmigem as a contributor for bug (#222)

* docs: update README.md

* docs: update .all-contributorsrc

---------

Co-authored-by: allcontributors[bot] <46447321+allcontributors[bot]@users.noreply.github.com>

* added contrast tutorial; better docs

* forgot to undo a change sum/effects

---------

Co-authored-by: allcontributors[bot] <46447321+allcontributors[bot]@users.noreply.github.com>

.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"
+      ]
     }
   ]
 }

README.md

@@ -143,6 +143,7 @@ You are very welcome to raise issues and start pull requests!
       <td align="center" valign="top" width="14.28%"><a href="https://github.com/suddha-bpn"><img src="https://avatars.githubusercontent.com/u/7974144?v=4?s=100" width="100px;" alt="suddha-bpn"/><br /><sub><b>suddha-bpn</b></sub></a><br /><a href="#bug-suddha-bpn" title="Bug reports">🐛</a></td>
       <td align="center" valign="top" width="14.28%"><a href="https://github.com/vladdez"><img src="https://avatars.githubusercontent.com/u/33777074?v=4?s=100" width="100px;" alt="Vladimir Mikheev"/><br /><sub><b>Vladimir Mikheev</b></sub></a><br /><a href="#bug-vladdez" title="Bug reports">🐛</a> <a href="#doc-vladdez" title="Documentation">📖</a></td>
       <td align="center" valign="top" width="14.28%"><a href="https://github.com/carmenamme"><img src="https://avatars.githubusercontent.com/u/100191854?v=4?s=100" width="100px;" alt="carmenamme"/><br /><sub><b>carmenamme</b></sub></a><br /><a href="#doc-carmenamme" title="Documentation">📖</a></td>
+      <td align="center" valign="top" width="14.28%"><a href="https://github.com/maxvanmigem"><img src="https://avatars.githubusercontent.com/u/115144441?v=4?s=100" width="100px;" alt="Maximilien Van Migem"/><br /><sub><b>Maximilien Van Migem</b></sub></a><br /><a href="#bug-maxvanmigem" title="Bug reports">🐛</a></td>
     </tr>
   </tbody>
 </table>

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).

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",

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))

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)
 ```
 
 """