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

Important

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)

6×5 DataFrame

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:08
  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