Visualize uncertainty in topoplots

Click to expand

using Base: channeled_tasks
using Unfold
using UnfoldMakie
using UnfoldSim
using DataFrames
using CairoMakie
using TopoPlots
using Statistics
using Random
using Animations

Representing uncertainty is one of the most difficult tasks in visualization. It is especially difficult for heatmaps and topoplots. Here we will present new ways to show uncertainty for topoplots.

Uncertainty in EEG data usually comes from subjects and trials:

Subjects can vary in phisological or behavioral characteristics;
Something can change between trials (electrode connection can get worse, etc.).

There are several measures of uncertainty. Here we will use standard deviation (x - mean(x) / N) and t-values (mean(x) / SE(x)).

Data input

Data for customized topoplots:

dat, positions = TopoPlots.example_data()
vec_estimate = dat[:, 340, 1];
vec_sd = dat[:, 340, 2]; # Here we use SD as uncertatinty measure

Data for animation:

df_toposeries, pos_toposeries =
    UnfoldMakie.example_data("bootstrap_toposeries"; noiselevel = 7);
df_toposeries1 = df_toposeries[df_toposeries.trial.<=15, :];
rng = MersenneTwister(1);

Please cite: HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce

This will generate DataFrame with 227 channels, 50 trials, 500 mseconds for bootstrapping. noiselevel is important for adding variability it your data.

Adjacent topoplot

In this case we already have two data vectors: vec_estimate with mean estimates and vec_sd with standard deviation.

Click to expand

begin
    f = Figure()
    ax = Axis(
        f[1, 1:2],
        title = "Time windows [340 ms]",
        titlesize = 24,
        titlealign = :center,
    )

    hidedecorations!(ax, label = false)
    hidespines!(ax)
    plot_topoplot!(
        f[1, 1],
        vec_estimate;
        positions = positions,
        visual = (; contours = false),
        axis = (; xlabel = ""),
        colorbar = (;
            label = "Voltage [µV]",
            labelsize = 24,
            ticklabelsize = 18,
            vertical = false,
            width = 180,
            position = :bottom,
        ),
    )
    plot_topoplot!(
        f[1, 2],
        vec_sd;
        positions = positions,
        visual = (; colormap = (:viridis), contours = false),
        axis = (; xlabel = "", xlabelsize = 24, ylabelsize = 24),
        colorbar = (;
            label = "Standard deviation",
            labelsize = 24,
            ticklabelsize = 18,
            vertical = false,
            width = 180,
            position = :bottom,
        ),
    )
end

GridPosition(GridLayout[1, 2] (3 children), GridLayoutBase.Span(1:1, 2:2), GridLayoutBase.Inner())

Uncertainty via marker size

We show uncertainty using donut-shaped electrode markers. The donut keeps the estimate color visible while the marker size reflects uncertainty — larger donuts mean higher uncertainty.

Click to expand

begin
    f = Figure()
    uncert_norm =
        (vec_sd .- minimum(vec_sd)) ./ (maximum(vec_sd) - minimum(vec_sd))
    uncert_scaled = uncert_norm * 30 .+ 10

    plot_topoplot!(
        f[1:4, 1],
        vec_estimate;
        positions,
        axis = (; xlabel = "Time point [340 ms]", xlabelsize = 24, ylabelsize = 24),
        topo_attributes = (;
            label_scatter = (; markersize = uncert_scaled, color = :transparent,
                strokecolor = :black,
                strokewidth = uncert_scaled .* 0.25)
        ),
        visual = (; colormap = :diverging_tritanopic_cwr_75_98_c20_n256, contours = false),
        colorbar = (; labelsize = 24, ticklabelsize = 18),
    )
    markersizes = round.(Int, range(extrema(uncert_scaled)...; length = 5))
    markerlables = round.(range(extrema(vec_sd)...; length = 5); digits = 2)

    group_size = [
        MarkerElement(
            marker = :circle,
            color = :transparent, strokecolor = :black, strokewidth = ms ÷ 5,
            markersize = ms) for ms in markersizes
    ]
    Legend(f[5, 1], group_size, ["$ms" for ms in markerlables], "Standard\ndeviation",
        patchsize = (maximum(markersizes) * 0.8, maximum(markersizes) * 0.8),
        framevisible = false,
        labelsize = 18, titlesize = 20,
        orientation = :horizontal, titleposition = :left, margin = (90, 0, 0, 0))
end

Legend()

Uncertainty via confidence intervals

Click to expand

begin
    f = Figure(; resolution = (900, 650))
    gf = f[1, 1] = GridLayout()
    pTopos = GridLayout()
    gf[1:2, 1:3] = pTopos

    pA = pTopos[1, 2]
    pB = pTopos[2, 1]
    pC = pTopos[2, 3]
    pcb = gf[:, 4]

    lims = begin
        p01, p99 = quantile(vec_estimate, [0.01, 0.99])
        m = max(abs(p01), abs(p99))
        Float32.((-m, m))
    end

    visual = (; limits = lims, colormap = cgrad(:RdBu, 10; categorical = true, rev = true))

    ticks5 = begin
        lo, hi = visual.limits
        pos = Float32[lo, lo/2, 0.0f0, hi/2, hi]
        lab = string.(round.(Float64.(pos); sigdigits = 2))
        (pos, lab)
    end

    plot_topoplot!(
        pA, vec_estimate;
        positions = positions,
        axis = (; xlabelsize = 24, xlabel = "Mean"),
        layout = (; use_colorbar = false),
        visual = visual,
    )

    plot_topoplot!(
        pB, vec_estimate .- vec_sd;
        positions = positions,
        axis = (; xlabelsize = 24, xlabel = "Mean - SE"),
        layout = (; use_colorbar = false),
        visual = visual,
    )

    plot_topoplot!(
        pC, vec_estimate .+ vec_sd;
        positions = positions,
        axis = (; xlabelsize = 24, xlabel = "Mean + SE"),
        layout = (; use_colorbar = false),
        visual = visual,
    )

    Colorbar(
        pcb;
        colormap = visual.colormap,
        limits = visual.limits,
        ticks = ticks5,
        label = "Voltage [µV]",
        labelsize = 24,
        ticklabelsize = 18,
        vertical = true,
        height = 300,
        flipaxis = true,
        labelrotation = -π/2,
    )

    rowgap!(pTopos, -90)
    colgap!(pTopos, -90)
    colgap!(gf, 10)
    colsize!(gf, 4, Auto(0.15))

end

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/kJl0u/src/scenes.jl:264

Uncertainty via animation

In this case, we need to boostrap the data, so we'll use raw data with single trials.

To show the uncertainty of the estimate, we will compute 10 different means of the boostrapped data. More specifically: 1) create N boostrapped data sets using random sampling with replacement across trials; 2) compute their means; 3) do a toposeries animation iterating over these means.

Click to expand for supportive functions

function for easing - smooth transition between frames in animation. update_ratio - transition ratio between time1 and time2. at - create animation object: 0 and 1 are time points, old and new are data vectors.

function ease_between(old, new, update_ratio; easing_function = sineio())
    anim = Animation(0, old, 1, new; defaulteasing = easing_function)
    return at(anim, update_ratio)
end

ease_between (generic function with 1 method)

Return (nchannels × nboot) matrix of bootstrap mean vectors, sampling independently per channel: μ + SE * randn().

function param_bootstrap_means(mean_vec::AbstractVector, se_vec::AbstractVector;
    n_boot::Int = 10, rng = MersenneTwister(1))
    T = float(promote_type(eltype(mean_vec), eltype(se_vec)))
    μ = convert(Vector{T}, mean_vec)
    se = convert(Vector{T}, se_vec)
    n_channels = length(μ)

    out = Matrix{T}(undef, n_channels, n_boot)
    for i_boot = 1:n_boot
        out[:, i_boot] = μ .+ se .* randn(rng, T, n_channels)
    end
    return out
end

param_bootstrap_means (generic function with 1 method)

n_boot = 20
vec_se = vec_sd ./ sqrt(15) # 15 subject according to paper
boot_means = param_bootstrap_means(vec_estimate, vec_se; n_boot = n_boot, rng = rng)
obs = Observable(boot_means[:, 1])
begin
    m = maximum(abs, boot_means)          # boot_means is nchan × n_boot
    vals = vec(boot_means)
    p01, p99 = quantile(vals, [0.01, 0.99])
    m = max(abs(p01), abs(p99))
    cr = Float32.((-m, m)) # symmetrical colorbars over zero
    f = Figure()
    plot_topoplot!(
        f[1, 1],
        obs;
        positions = positions,
        visual = (; contours = false, colormap = cgrad(:RdYlBu, 10; categorical = true), colorrange = cr),
        colorbar = (; labelsize = 24, ticklabelsize = 18, height = 350),
        axis = (; xlabel = "Time window 100", xlabelsize = 24, ylabelsize = 24),
    )
    f
end
record(f, "bootstrap_single_topo.mp4"; framerate = 12) do io
    recordframe!(io)  # first frame (original)
    for i_boot = 1:(n_boot-1)          # number of bootstrap targets
        new_v = boot_means[:, i_boot+1]
        old_v = copy(obs[])
        for u in range(0, 1, length = 10)   # easing steps
            obs[] = ease_between(old_v, new_v, u)
            recordframe!(io)
        end
    end
end

"bootstrap_single_topo.mp4"

This page was generated using Literate.jl.