Skip to content

Commit

Permalink
merge free
Browse files Browse the repository at this point in the history
  • Loading branch information
behinger committed Apr 16, 2024
2 parents da88104 + 7d1a06c commit 2066c52
Show file tree
Hide file tree
Showing 6 changed files with 138 additions and 148 deletions.
44 changes: 21 additions & 23 deletions docs/src/HowTo/custom_solvers.md
Original file line number Diff line number Diff line change
@@ -1,48 +1,46 @@
# [Custom Solvers](@id custom_solvers)
# [Alternative Solvers](@id custom_solvers)

### Setup some data

```@Example main
using Unfold
using UnfoldMakie, CairoMakie
using UnfoldSim,
dat,evts = UnfoldSim.predef_eeg(;noiselevel=10,return_epoched=true)
f = @formula 0~1+condition+continuous
designDict = Dict(Any=>(f,range(0,1,length=size(dat,1))))
using UnfoldSim
dat, evts = UnfoldSim.predef_eeg(; noiselevel = 10, return_epoched = true)
f = @formula 0 ~ 1 + condition + continuous
designDict = Dict(Any => (f, range(0, 1, length = size(dat, 1))))
```

### GPU Solvers
GPU solvers can speed up your modelfit drastically! up to factor of 30 has been observed already
GPU solvers can significantly speed up your model fitting, with observed improvements of up to a factor of 30!

```julia
using Krylov,CUDA # necessary to load the right package extension
gpu_solver =(x,y)->Unfold.solver_krylov(x,y;GPU=true)
m = Unfold.fit(UnfoldModel,designDict,evts,dat,solver=gpu_solver)
using Krylov, CUDA # necessary to load the right package extension
gpu_solver =(x, y) -> Unfold.solver_krylov(x, y; GPU = true)
m = Unfold.fit(UnfoldModel, designDict, evts, dat, solver = gpu_solver)
```
We can't run it on the docs though, so try it yourself! If you need something else than CUDA, write an issue, we cant test it with something else right now...

To test it, you will need to run it yourself as we cannot run it on the docs. If you require a different graphicscard vendor than NVIDA/CUDA, please create an issue. Currently, we are unable to test it due to lack of hardware.

### Robust Solvers
Robust solvers automatically account for outlying trials. They come at a severe computational cost though!
Robust solvers automatically account for outlier trials, but they come at a significant computational cost.
```@Example main
using RobustModels # necessary to load the Unfold st
package extension
se_solver =(x,y)->Unfold.solver_robust(x,y)
m = Unfold.fit(UnfoldModel,designDict,evts,dat,solver=se_solver)
results =coeftable(m)
plot_erp(results;stderror=true)
using RobustModels # necessary to load the Unfold package extension
se_solver = (x, y) -> Unfold.solver_robust(x, y)
m = Unfold.fit(UnfoldModel, designDict, evts, dat, solver = se_solver)
results = coeftable(m)
plot_erp(results; stderror = true)
```

### Back2Back regression
```@Example main
b2b_solver = (x, y) -> Unfold.solver_b2b(x, y;ross_val_reps = 5)
dat_3d = permutedims(repeat(dat,1,1,20),[3 1 2])
m = Unfold.fit(UnfoldModel, designDict, evts,dat_3d; solver=b2b_solver)
b2b_solver = (x, y) -> Unfold.solver_b2b(x, y; ross_val_reps = 5)
dat_3d = permutedims(repeat(dat, 1, 1, 20), [3 1 2])
m = Unfold.fit(UnfoldModel, designDict, evts, dat_3d; solver = b2b_solver)
results = coeftable(m)
plot_erp(results)
```
These are the decoding-results for `conditionA` while taking into account `conditionB` - and vice versa.
These are the decoding results for `conditionA` while considering `conditionB`, and vice versa.


36 changes: 22 additions & 14 deletions docs/src/HowTo/multiple_events.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
# How to model multiple events

In case of overlapping data, we often need to plot multiple events.

When dealing with overlapping data, it is often necessary to plot multiple events.

### Load Example Data
```@example main
Expand All @@ -11,34 +10,43 @@ using DataFrames
using StatsModels
using MixedModels
include(joinpath(dirname(pathof(Unfold)), "../test/test_utilities.jl") ) # to load data
include(joinpath(dirname(pathof(Unfold)), "../test/test_utilities.jl")) # to load data
dat, evts = loadtestdata("test_case_4b");
evts[1:5,:]
```
As you can see, there are two events here. EventA and EventB. Both are in the column `type` in which the toolbox looks for different events by default.
The `type` column of table `evts` contains two conditions: EventA and EventB. By default, the toolbox will search for these conditions.

### Specify formulas and basisfunctions

```@example main
bf1 = firbasis(τ=(-0.4,.8),sfreq=50)
bf2 = firbasis(τ=(-0.2,1.2),sfreq=50)
bf1 = firbasis(τ = (-0.4, 0.8), sfreq = 50)
bf2 = firbasis(τ = (-0.2, 1.2), sfreq = 50)>>>>>>> 7d1a06c0f1c6aa93ce83f7d828a317c45dfaa674
```
For each event, we have to specify a basisfunction and a formula. We could use the same basis and the same formulas though
For each event, a basis function and formula must be specified. The same basis and formulas may be used.
```@example main
f = @formula 0~1
f = @formula 0 ~ 1
```

Next, we have to specify for each event, what is the formula and what is the basisfunction we want to use
For each event, we must specify the formula and basis function to be used.
```@example main
bfDict = ["eventA"=>(f,bf1),
"eventB"=>(f,bf2)]
bfDict = [ "eventA" => (f, bf1),
"eventB" => (f, bf2) ]
```

Finally, fitting & plotting works the same way as always
```@example main

m = Unfold.fit(UnfoldModel,bfDict,evts,dat;solver=(x,y) -> Unfold.solver_default(x,y;stderror=true),eventcolumn="type")
```@example main
m = Unfold.fit(
UnfoldModel,
bfDict,
evts,
dat,
solver = (x, y) -> Unfold.solver_default(x, y; stderror = true),
eventcolumn = "type",
)
results = coeftable(m)
plot_erp(results;stderror=true,mapping=(;col=:group))
plot_erp(results; stderror = true, mapping = (; col = :group))
```
21 changes: 10 additions & 11 deletions docs/src/HowTo/standarderrors.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,22 +5,21 @@
```@Example main
using Unfold
using UnfoldMakie, CairoMakie
using UnfoldSim,
dat,evts = UnfoldSim.predef_eeg(;noiselevel=10,return_epoched=true)
f = @formula 0~1+condition+continuous
designDict = Dict(Any=>(f,range(0,1,length=size(dat,1))))
using UnfoldSim
dat, evts = UnfoldSim.predef_eeg(; noiselevel = 10, return_epoched = true)
f = @formula 0 ~ 1 + condition + continuous
designDict = Dict(Any => (f, range(0, 1, length = size(dat, 1))))
```

Similar to @Ref(custom_solver), it is possible to specify a solver than also calculates the standard errors of the (single-subject) estimates.
It is possible to specify a solver that calculates the standard errors of the estimates for a single subject as it possible for [custom solvers](@ref custom_solvers).

```@Example main
se_solver =(x,y)->Unfold.solver_default(x,y,stderror=true)
m = Unfold.fit(UnfoldModel,designDict,evts,dat,solver=se_solver)
results =coeftable(m)
plot_erp(results;stderror=true)
se_solver = (x, y) -> Unfold.solver_default(x, y, stderror = true)
m = Unfold.fit(UnfoldModel, designDict, evts, dat, solver = se_solver)
results = coeftable(m)
plot_erp(results; stderror = true)
```
!!! warning
**In case of overlap-correction:** Use single-subject SE on your own risk. Because EEG data are autocorrelated your standarderrors will be typically too small.
**In case of overlap-correction:** Use single-subject standard errors on your own risk. EEG data is autocorrelated, which means that standard errors are typically too small.

73 changes: 38 additions & 35 deletions docs/src/tutorials/lm_mu.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# [Mass Univariate Linear Models (no overlap correction)](@id lm_massunivariate)

In this notebook we will fit regression models to (simulated) EEG data. We will see that we need some type of overlap correction, as the events are close in time to each other, so that the respective brain responses overlap.
In this notebook we will fit regression models to simulated EEG data. We will see that we need some type of overlap correction, as the events are close in time to each other, so that the respective brain responses overlap.
If you want more detailed introduction to this topic check out [our paper](https://peerj.com/articles/7838/).


Expand All @@ -9,7 +9,7 @@ If you want more detailed introduction to this topic check out [our paper](https
```@example Main
using DataFrames
using Unfold
using UnfoldMakie,CairoMakie # for plotting
using UnfoldMakie, CairoMakie # for plotting
using UnfoldSim
Expand All @@ -23,40 +23,43 @@ data, evts = UnfoldSim.predef_eeg()
nothing # hide
```
## Inspection
The data has little noise and the underlying signal is a pos-neg-pos spike pattern
The data has only little noise. The underlying signal pattern is a positive-negative-positive spike.
```@example Main
times = range(1/50,length=200,step=1/50)
plot(times,data[1:200])
vlines!(current_axis(),evts[evts.latency.<=200,:latency]./50) # show events, latency in samples!
times = range(1/50, length=200, step=1/50)
f = Figure()
plot(f[1, 1], data[1:200], times)
vlines!(evts[evts.latency .<= 200, :latency] ./ 50) # show events, latency in samples!
f
```

To inspect the event dataframe we use
```@example Main
show(first(evts,6,),allcols=true)
show(first(evts, 6), allcols = true)
```
Every row is an event. Note that `:latency` is commonly the timestamp in samples, whereas `:onset` would typically refer to seconds.
Every row is an experimental event. Note that `:latency` refers to time in samples, whereas `:onset` would typically refer to seconds.


## Traditional Mass Univariate Analysis
For a mass univariate analysis, will
1. epoch the data
2. specify a formula
3. fit a linear model to each time point & channel
4. visualize the results.
To perform a mass univariate analysis, you must complete the following steps:

1. Split data into epochs
2. Specify a formula
3. Fit a linear model to each time point & channel
4. Visualize the results.

#### 1. Epoch the data
We have to cut the data into small regular trials.

#### 1. Split data into epochs
Initially, you have data with a duration that represents the whole experimental trial. You need to cut the data into small regular epochs related to the some event, e.g. start of fixation.

```@example Main
# we have multi channel support
data_r = reshape(data,(1,:))
data_r = reshape(data, (1,:))
# cut the data into epochs
data_epochs,times = Unfold.epoch(data=data_r,tbl=evts,τ=(-0.4,0.8),sfreq=50);
data_epochs, times = Unfold.epoch(data = data_r, tbl = evts, τ = (-0.4, 0.8), sfreq = 50);
size(data_epochs)
```

`τ` specifies the epoch size. To convert `τ` to samples, we need to specify the sampling rate.
- `τ` specifies the epoch size.
- `sfreq` - sampling rate, converts `τ` to samples.


```@example Main
Expand All @@ -66,48 +69,48 @@ typeof(data_epochs)
In julia, `missing` is supported throughout the ecosystem. Thus, we can have partial trials and they will be incorporated / ignored at the respective functions. Helpful functions are the julia-base `disallowmissing` and the internal `Unfold.drop_missing_epochs` functions



#### 2. Specify a formula
We define a formula that we want to apply to each point in time
Define a formula to be applied to each time point.

```@example Main
f = @formula 0~1+condition+continuous # 0 as a dummy, we will combine wit data later
f = @formula 0 ~ 1 + condition + continuous # 0 as a dummy, we will combine with data later
nothing # hide
```

#### 3. fit a linear model to each time-point & channel
#### 3. Fit a linear model to each time point & channel

We fit the `UnfoldLinearModel`.
Fit the `UnfoldLinearModel`.
```@example Main
m = fit(UnfoldModel,f,evts,data_epochs,times);
m = fit(UnfoldModel, f, evts, data_epochs, times);
nothing #hide
```

Alternatively, we could also call it using:
Alternative way to call this model is below. This syntax allows you to fit multiple events at once. For example, replacing `Any` with `:fixation =>...` will fit this model specifically to the fixation event type.
```@example Main
m = fit(UnfoldModel,Dict(Any=>(f,times)),evts,data_epochs);
m = fit(UnfoldModel, Dict(Any=>(f, times)), evts, data_epochs);
nothing #hide
```
A syntax which allows for fitting multiple events at the same time (by replacing Any with e.g. `:fixation =>...` it will fit that model specifically to the fixation event type)

We can inspect the object
Inspect the fitted model:
```@example Main
m
```
And see that there are several helpful functions to recover `design`, `designmatrix`, raw-`modelfit` and most importantly, `coeftable`.
Note these functions to discover the model: `design`, `designmatrix`, `modelfit` and most importantly, `coeftable`.

!!! info
There are of course further methods, e.g. `coef`, `ranef`, `Unfold.formula`, `modelmatrix` which might be helpful at some point, but not important now
There are of course further methods, e.g. `coef`, `ranef`, `Unfold.formula`, `modelmatrix` which might be helpful at some point, but not important now.

Using `coeftable`, we can get a *tidy*-dataframe, very useful for your further analysis.
Using `coeftable`, we can get a *tidy* DataFrames, very useful for your further analysis.

```@example Main
first(coeftable(m),6)
first(coeftable(m), 6)
```

#### 4. Visualize the results
Tidy-Dataframes make them easy to visualize using e.g. AlgebraOfGraphics.jl. We simpliefied this even more, by implementing a wrapper in `UnfoldMakie`
Tidy DataFrames are easy to visualize using e.g. AlgebraOfGraphics.jl. Function `plot_erp` from `UnfoldMakie`makes it even easier.

```@example Main
results = coeftable(m)
plot_erp(results)
```
As you can see here, a lot is going on, even in the baseline-period! This is because the signal was simulated with overlapping events. Head over to the next tutorial to find out how to remedy this.
As you can see, there is a lot going on, even in the baseline period! This is because the signal was simulated with overlapping events. In the next tutorial you will learn how to fix this.
Loading

0 comments on commit 2066c52

Please sign in to comment.