using ClusterDepth
using Random
using CairoMakie
using UnfoldSim
using Unfold
using UnfoldMakie
using Statistics
Precompiling Unfold...
Info Given Unfold was explicitly requested, output will be shown live 
WARNING: Method definition check_data(Type{Unfold.ContinuousTimeTrait{UF}}, Type{UF}, AbstractArray{T, 1}) where {T, UF<:(Unfold.UnfoldModel{T} where T)} in module Unfold at /home/runner/.julia/packages/Unfold/i00Gm/src/fit.jl:225 overwritten at /home/runner/.julia/packages/Unfold/i00Gm/src/fit.jl:231.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
   4094.3 ms  ? Unfold
WARNING: Method definition check_data(Type{Unfold.ContinuousTimeTrait{UF}}, Type{UF}, AbstractArray{T, 1}) where {T, UF<:(Unfold.UnfoldModel{T} where T)} in module Unfold at /home/runner/.julia/packages/Unfold/i00Gm/src/fit.jl:225 overwritten at /home/runner/.julia/packages/Unfold/i00Gm/src/fit.jl:231.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
Precompiling UnfoldBSplineKitExt...
   4033.6 ms  ? Unfold
Info Given UnfoldBSplineKitExt was explicitly requested, output will be shown live 
┌ Warning: Module Unfold with build ID ffffffff-ffff-ffff-26ff-6647cd2eb0f9 is missing from the cache.
│ This may mean Unfold [181c99d8-e21b-4ff3-b70b-c233eddec679] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
    744.5 ms  ? Unfold → UnfoldBSplineKitExt
┌ Warning: Module Unfold with build ID ffffffff-ffff-ffff-26ff-6647cd2eb0f9 is missing from the cache.
│ This may mean Unfold [181c99d8-e21b-4ff3-b70b-c233eddec679] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
Precompiling UnfoldMakie...
   3936.6 ms  ? Unfold
    631.0 ms  ? Unfold → UnfoldBSplineKitExt
Info Given UnfoldMakie was explicitly requested, output will be shown live 
┌ Warning: Module Unfold with build ID ffffffff-ffff-ffff-26ff-6647cd2eb0f9 is missing from the cache.
│ This may mean Unfold [181c99d8-e21b-4ff3-b70b-c233eddec679] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
   7167.7 ms  ? UnfoldMakie
┌ Warning: Module Unfold with build ID ffffffff-ffff-ffff-26ff-6647cd2eb0f9 is missing from the cache.
│ This may mean Unfold [181c99d8-e21b-4ff3-b70b-c233eddec679] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
Precompiling UnfoldMakieUnfoldSimExt...
   4035.0 ms  ? Unfold
    639.0 ms  ? Unfold → UnfoldBSplineKitExt
   7079.4 ms  ? UnfoldMakie
Info Given UnfoldMakieUnfoldSimExt was explicitly requested, output will be shown live 
┌ Warning: Module UnfoldMakie with build ID ffffffff-ffff-ffff-11e4-09ee84bc74cc is missing from the cache.
│ This may mean UnfoldMakie [69a5ce3b-64fb-4f22-ae69-36dd4416af2a] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
   7130.7 ms  ? UnfoldMakie → UnfoldMakieUnfoldSimExt
┌ Warning: Module UnfoldMakie with build ID ffffffff-ffff-ffff-11e4-09ee84bc74cc is missing from the cache.
│ This may mean UnfoldMakie [69a5ce3b-64fb-4f22-ae69-36dd4416af2a] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541

How to use the ClusterDepth multiple comparison correction on multichannel data

This tutorial is adapted from the single-channel EEG example, and complements it with the HArtMuT NYhead model (https://github.com/harmening/HArtMuT) to simulate multiple channels.

First set up the EEG simulation as before, with one subject and 40 trials:

design =
    SingleSubjectDesign(conditions=Dict(:condition => ["car", "face"])) |>
    x -> RepeatDesign(x, 40);
p1 = LinearModelComponent(;
    basis=p100(; sfreq=250),
    formula=@formula(0 ~ 1),
    β=[1.0],
);

n170 = LinearModelComponent(;
    basis=UnfoldSim.n170(; sfreq=250),
    formula=@formula(0 ~ 1 + condition),
    β=[1.0, 0.5], # condition effect - faces are more negative than cars
);
p300 = LinearModelComponent(;
    basis=UnfoldSim.p300(; sfreq=250),
    formula=@formula(0 ~ 1 + condition),
    β=[1.0, 0], # no p300 condition effect
);

Now choose some source coordinates for each of the p100, n170, p300 that we want to simulate, and use the helper function closest_srcs to get the HArtMuT sources that are closest to these coordinates:

src_coords = [
    [20, -78, -10], #p100
    [-20, -78, -10], #p100
    [50, -40, -25], #n170
    [0, -50, 40], #p300
    [0, 5, 20], #p300
];

headmodel_HArtMuT = headmodel()
get_closest =
    coord ->
        UnfoldSim.closest_src(coord, headmodel_HArtMuT.cortical["pos"]) |>
        pi -> magnitude(headmodel_HArtMuT)[:, pi]

p1_l = p1 |> c -> MultichannelComponent(c, get_closest([-20, -78, -10]))
p1_r = p1 |> c -> MultichannelComponent(c, get_closest([20, -78, -10]))
n170_r = n170 |> c -> MultichannelComponent(c, get_closest([50, -40, -25]))
p300_do = p300 |> c -> MultichannelComponent(c, get_closest([0, -50, -40]))
p300_up = p300 |> c -> MultichannelComponent(c, get_closest([0, 5, 20]))

data, evts = simulate(
    MersenneTwister(1),
    design,
    [p1_l, p1_r, n170_r, p300_do, p300_up],
    UniformOnset(; offset=0.5 * 250, width=100),
    RedNoise(noiselevel=1);
    return_epoched=true,
);
Please cite: HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce

Plotting

This is what the data looks like, for one channel/trial respectively:

f = Figure()
Axis(f[1, 1], title="Single channel, all trials", xlabel="time", ylabel="y")
series!(data[1, :, :]', solid_color=(:black, 0.1))
lines!(mean(data[1, :, :], dims=2)[:, 1], color=:red)
hlines!([0], color=:gray)

Axis(f[2, 1], title="All channels, average over trials", xlabel="time", ylabel="y")
series!(mean(data, dims=3)[:, :, 1], solid_color=(:black, 0.1))
hlines!([0], color=:gray)
f
Example block output

And some topoplots:

positions = [
    Point2f(p[1] + 0.5, p[2] + 0.5) for
    p in to_positions(headmodel_HArtMuT.electrodes["pos"]')
]

df = UnfoldMakie.eeg_array_to_dataframe(
    mean(data, dims=3)[:, :, 1],
    string.(1:length(positions)),
);
Δbin = 20 # 20 samples / bin
plot_topoplotseries(
    df,
    bin_width=20,
    positions=positions,
    visual=(; enlarge=1, label_scatter=false),
)
Example block output

ClusterDepth

Now that the simulation is done, let's try out ClusterDepth and plot our results Note that this is a simple test of "activity" vs. 0

pvals = clusterdepth(data; τ=1.6, nperm=200);
fig, ax, hm = heatmap(transpose(pvals), colorscale=log10)
ax.title = "pvals";
ax.xlabel = "time";
ax.ylabel = "channel";
Colorbar(fig[:, end+1], hm);
fig
Example block output

This page was generated using Literate.jl.