Mass Univariate Linear Mixed Models
using Unfold
using UnfoldSim
using MixedModels # important to load to activate the UnfoldMixedModelsExtension
using UnfoldMakie, CairoMakie # plotting
using DataFrames
using CategoricalArrays
You have to run using MixedModels
before or after loading Unfold to activate the MixedModels abilities!
This notebook is similar to the Mass Univariate Linear Models (no overlap correction) tutorial, but fits mass-univariate mixed models - that is, one model over all subjects, instead of one model per subject. This allows to include item effects, for example.
Mass Univariate Mixed Models
Again we have 4 steps:
- Split data into epochs
- Specify a formula
- Fit a linear model to each time point & channel
- Visualize the results.
1. Epoching
data, evts = UnfoldSim.predef_eeg(10; return_epoched = true) # simulate 10 subjects
data = reshape(data, 1, size(data, 1), :) # concatenate the data into a long EEG dataset
times = range(0, length = size(data, 2), step = 1 / 100)
transform!(evts, :subject => categorical => :subject); # :subject must be categorical, otherwise MixedModels.jl complains
The events
dataFrame has an additional column (besides being much taller): subject
first(evts, 6)
Row | subject | item | continuous | condition | latency |
---|---|---|---|---|---|
Cat⦠| String | Float64 | String | Int64 | |
1 | S01 | I032 | -3.88889 | face | 64 |
2 | S01 | I003 | -2.77778 | car | 131 |
3 | S01 | I064 | -1.66667 | car | 183 |
4 | S01 | I024 | -1.66667 | car | 233 |
5 | S01 | I002 | -3.88889 | car | 284 |
6 | S01 | I013 | -2.77778 | face | 354 |
2. Formula specification
We define the formula. Importantly, we need to specify a random effect. We use zerocorr
to speed up the calculation.
f = @formula 0 ~ 1 + condition * continuous + zerocorr(1 + condition * continuous | subject);
3. Model fitting
We can now run the LinearMixedModel at each time point.
m = fit(UnfoldModel, f, evts, data, times)
Progress: 4%|ββ | ETA: 0:00:09
channel: 1
time: 2
Progress: 100%|βββββββββββββββββββββββββββββββββββββββββ| Time: 0:00:00
channel: 1
time: 45
4. Visualization of results
Let's start with the fixed effects. We see the condition effects and some residual overlap activity in the fixed effects.
results = coeftable(m)
res_fixef = results[isnothing.(results.group), :]
plot_erp(res_fixef)
And now comes the random effect:
res_ranef = results[results.group .== :subject, :]
plot_erp(res_ranef)
Statistics
Check out the LMM p-value tutorial