Merge pull request #11 from s-ccs/docfixes

Docfixes

Author: behinger

Project.toml

@@ -4,13 +4,11 @@ authors = ["Benedikt V. Ehinger", "Maanik Marathe"]
 version = "0.2.1"
 
 [deps]
-Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
-Images = "0.25"
 Random = "1"
 SparseArrays = "1"
 StatsBase = "0.33, 0.34"

docs/make.jl

@@ -2,7 +2,7 @@ using ClusterDepth
 using Documenter
 using Glob
 using Literate
-DocMeta.setdocmeta!(ClusterDepth, :DocTestSetup, :(using ClusterDepth); recursive=true)
+DocMeta.setdocmeta!(ClusterDepth, :DocTestSetup, :(using ClusterDepth); recursive = true)
 
 
 GENERATED = joinpath(@__DIR__, "src")
@@ -13,17 +13,17 @@ for subfolder ∈ ["explanations", "HowTo", "tutorials", "reference"]
 
 end
 makedocs(;
-    modules=[ClusterDepth],
-    authors="Benedikt V. Ehinger, Maanik Marathe",
-    repo="https://github.com/s-ccs/ClusterDepth.jl/blob/{commit}{path}#{line}",
-    sitename="ClusterDepth.jl",
-    format=Documenter.HTML(;
-        prettyurls=get(ENV, "CI", "false") == "true",
-        canonical="https://s-ccs.github.io/ClusterDepth.jl",
-        edit_link="main",
-        assets=String[],
+    modules = [ClusterDepth],
+    authors = "Benedikt V. Ehinger, Maanik Marathe",
+    repo = "https://github.com/s-ccs/ClusterDepth.jl/blob/{commit}{path}#{line}",
+    sitename = "ClusterDepth.jl",
+    format = Documenter.HTML(;
+        prettyurls = get(ENV, "CI", "false") == "true",
+        canonical = "https://s-ccs.github.io/ClusterDepth.jl",
+        edit_link = "main",
+        assets = String[],
     ),
-    pages=[
+    pages = [
         "Home" => "index.md",
         "Tutorials" => [
             "An EEG Example" => "tutorials/eeg.md",
@@ -36,7 +36,4 @@ makedocs(;
     ],
 )
 
-deploydocs(;
-    repo="github.com/s-ccs/ClusterDepth.jl",
-    devbranch="main",
-)
+deploydocs(; repo = "github.com/s-ccs/ClusterDepth.jl", devbranch = "main")

docs/src/reference/type1.jl

@@ -15,7 +15,11 @@ using DataFrames
 # ## Setup Simulation
 # Let's setup a simulation using UnfoldSim.jl. We simulate a simple 1x2 design with 20 subjects
 n_subjects = 20
-design = MultiSubjectDesign(n_subjects=n_subjects, n_items=2, items_between=Dict(:condition => ["small", "large"]))
+design = MultiSubjectDesign(
+    n_subjects=n_subjects,
+    n_items=2,
+    items_between=Dict(:condition => ["small", "large"]),
+)
 first(generate_events(design), 5)
 
 # 
@@ -33,19 +37,32 @@ signal = MixedModelComponent(;
 # Let's move the actual simulation into a function, so we can call it many times.
 # Note that we use (`RedNoise`)[https://unfoldtoolbox.github.io/UnfoldSim.jl/dev/literate/reference/noisetypes/] which has lot's of Autocorrelation between timepoints. nice!
 function run_fun(r)
-    data, events = simulate(MersenneTwister(r), design, signal, UniformOnset(; offset=5, width=4), RedNoise(noiselevel=1); return_epoched=true)
+    data, events = simulate(
+        MersenneTwister(r),
+        design,
+        signal,
+        UniformOnset(; offset=5, width=4),
+        RedNoise(noiselevel=1);
+        return_epoched=true,
+    )
     data = reshape(data, size(data, 1), :)
     data = data[:, events.condition.=="small"] .- data[:, events.condition.=="large"]
 
-    return data, clusterdepth(data'; τ=quantile(TDist(n_subjects - 1), 0.95), nperm=1000)
+    return data,
+    clusterdepth(data'; τ=quantile(TDist(n_subjects - 1), 0.95), nperm=1000)
 end;
 
 # ## Understanding the simulation
 # let's have a look at the actual data by running it once, plotting condition wise trials, the ERP and histograms of uncorrected and corrected p-values
 data, pval = run_fun(5)
 conditionSmall = data[:, 1:2:end]
 conditionLarge = data[:, 2:2:end]
-pval_uncorrected = 1 .- cdf.(TDist(n_subjects - 1), abs.(ClusterDepth.studentt(conditionSmall .- conditionLarge)))
+pval_uncorrected =
+    1 .-
+    cdf.(
+        TDist(n_subjects - 1),
+        abs.(ClusterDepth.studentt(conditionSmall .- conditionLarge)),
+    )
 sig = pval_uncorrected .<= 0.025;
 
 # For the uncorrected p-values based on the t-distribution, we get a type1 error over "time":
@@ -67,8 +84,8 @@ sig = allowmissing(sig)
 sig[sig.==0] .= missing
 @show sum(skipmissing(sig))
 lines!(sig, color=:gray, linewidth=4)
-lines!(ax, mean(conditionSmall, dims=2)[:, 1], solid_color=:red)
-lines!(ax, mean(conditionLarge, dims=2)[:, 1], solid_color=:blue)
+lines!(ax, mean(conditionSmall, dims=2)[:, 1], color=:red)
+lines!(ax, mean(conditionLarge, dims=2)[:, 1], color=:blue)
 
 hist!(Axis(f[3, 1], title="uncorrected pvalues"), pval_uncorrected, bins=0:0.01:1.1)
 hist!(Axis(f[3, 2], title="clusterdepth corrected pvalues"), pval, bins=0:0.01:1.1)
@@ -81,12 +98,13 @@ res = fill(NaN, reps, 2)
 Threads.@threads for r = 1:reps
     data, pvals = run_fun(r)
     res[r, 1] = mean(pvals .<= 0.05)
-    res[r, 2] = mean(abs.(ClusterDepth.studentt(data)) .>= quantile(TDist(n_subjects - 1), 0.975))
+    res[r, 2] =
+        mean(abs.(ClusterDepth.studentt(data)) .>= quantile(TDist(n_subjects - 1), 0.975))
 end;
 # Finally, let's calculate the percentage of simulations where we find a significant effect somewhere
 mean(res .> 0, dims=1) |> x -> (:clusterdepth => x[1], :uncorrected => x[2])
 
-# Nice, correction seems to work in principle :) Clusterdepth is not exactly 5%, but with more repetitions we should get there (e.g. with 5000 repetitions, we get 0.051%).
+# Nice, correction seems to work in principle :) Clusterdepth is not necessarily exactly 5%, but with more repetitions we should get there (e.g. with 5000 repetitions, we got 0.051%).
 
 # !!! info
-#       if you look closely, the `:uncorrected` value (around 60%) is not as bad as the 99% promised in the introduction. This is due to the correlation between the tests introduced by the noise. Indeed, a good exercise is to repeat everything, but put `RedNoise` to `WhiteNoise`
+#       if you look closely, the `:uncorrected` value (can be around 60%) is not as bad as the 99% promised in the introduction. This is due to the correlation between the tests introduced by the noise. Indeed, a good exercise is to repeat everything, but put `RedNoise` to `WhiteNoise`

docs/src/reference/type1_troendle.jl

@@ -50,23 +50,22 @@ end;
 res = any(pvals_all[:, :, :] .<= 0.05, dims=3)[:, :, 1]
 mean(res .> 0, dims=1) |> x -> (:troendle => x[1], :uncorrected => x[2])
 
-# Nice. Troendle fits perfectly and the uncorrected is pretty close to what we calculated above!
+# Nice. Troendle corrects the data and we are below 0.05. The  uncorrected simulations are close to what we calculated above!
 
 # Finally we end this with a short figure to get a better idea of how this data looks like and a histogram of the p-values
 
 f = Figure()
-ax = f[1, 1] = Axis(f)
+ax = f[1, 1] = Axis(f, title="t-values")
 
 
 
 lines!(ax, abs.(ClusterDepth.studentt(data)))
-heatmap!(Axis(f[2, 1]), data)
-series!(Axis(f[2, 2]), data[:, 1:7]')
-h1 = scatter!(Axis(f[1, 2]; yscale=log10), pvals, label="troendle")
+heatmap!(Axis(f[1, 2], title="heatmap of data"), data)
+series!(Axis(f[2, 2], title="data: subset of left plot"), data[:, 1:7]')
+h1 = scatter!(Axis(f[2, 1]; yscale=log10, title="troendle pvals"), pvals, label="troendle")
 
 hlines!([0.05, 0.01])
 
-hist!(Axis(f[3, 1]), pvals_all[:, 1, :][:], bins=0:0.01:1.1)
-hist!(Axis(f[3, 2]), pvals_all[:, 2, :][:], bins=0:0.01:1.1)
+hist!(Axis(f[3, 1], title="troendle corrected"), pvals_all[:, 1, :][:], bins=0:0.01:1.1)
+hist!(Axis(f[3, 2], title="uncorrected"), pvals_all[:, 2, :][:], bins=0:0.01:1.1)
 f
-

docs/src/tutorials/demo.jl

@@ -3,34 +3,34 @@ using Random
 using CairoMakie
 
 
-n_t =40 # timepoints
+n_t = 40 # timepoints
 n_sub = 50
 n_perm = 5000
 
 snr = 0.5 # signal to nois
 
 ## add a signal to the middle
-signal = vcat(zeros(n_t÷4), sin.(range(0,π,length=n_t÷2)), zeros(n_t÷4))
+signal = vcat(zeros(n_t ÷ 4), sin.(range(0, π, length = n_t ÷ 2)), zeros(n_t ÷ 4))
 
 ## same signal for all subs
-signal = repeat(signal,1,n_sub)
+signal = repeat(signal, 1, n_sub)
 
 
 ## add noise
-data = randn(MersenneTwister(123),n_t,n_sub).+ snr .* signal
+data = randn(MersenneTwister(123), n_t, n_sub) .+ snr .* signal
 
 ## by default assumes τ=2.3 (~alpha=0.05), and one-sample ttest
 @time pvals = clusterdepth(data);
 
 f = Figure()
-ax = f[1,1] = Axis(f)
+ax = f[1, 1] = Axis(f)
 
 
 lines!(abs.(ClusterDepth.studentt(data)))
-h1 = scatter(f[1,2],pvals;axis=(;yscale=log10),label="troendle")
+h1 = scatter(f[1, 2], pvals; axis = (; yscale = log10), label = "troendle")
 
-pvals2 = clusterdepth(data;pval_type=:naive)
-h2 = scatter!(1.2:40.2,pvals2,color="red",label="naive")
+pvals2 = clusterdepth(data; pval_type = :naive)
+h2 = scatter!(1.2:40.2, pvals2, color = "red", label = "naive")
 #hlines!(([0.05]))
 axislegend()
-f
+f

docs/src/tutorials/eeg-multichannel.jl

@@ -6,16 +6,18 @@ using Unfold
 using UnfoldMakie
 using Statistics
 
-# ## How to use clusterDepth multiple comparison correction on multichannel data
+# # How to use the ClusterDepth multiple comparison correction on multichannel data
 
-# This tutorial is adapted from the first EEG example and uses the HArtMuT NYhead model (https://github.com/harmening/HArtMuT) to simulate multiple channels. 
+# 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);
+design =
+    SingleSubjectDesign(conditions=Dict(:condition => ["car", "face"])) |>
+    x -> RepeatDesign(x, 40);
 p1 = LinearModelComponent(;
     basis=p100(; sfreq=250),
     formula=@formula(0 ~ 1),
-    β=[1.0]
+    β=[1.0],
 );
 
 n170 = LinearModelComponent(;
@@ -39,46 +41,66 @@ src_coords = [
 ];
 
 headmodel_HArtMuT = headmodel()
-get_closest = coord -> UnfoldSim.closest_src(coord, headmodel_HArtMuT.cortical["pos"]) |> pi -> magnitude(headmodel_HArtMuT; type="perpendicular")[:, pi]
+get_closest =
+    coord ->
+        UnfoldSim.closest_src(coord, headmodel_HArtMuT.cortical["pos"]) |>
+        pi -> magnitude(headmodel_HArtMuT; type="perpendicular")[:, 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, events = simulate(MersenneTwister(1), design, [p1_l, p1_r, n170_r, p300_do, p300_up],
+data, events = 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);
+    RedNoise(noiselevel=1);
+    return_epoched=true,
+);
 
 
 # ## 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)
+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)
+series!(mean(data, dims=3)[:, :, 1], solid_color=(:black, 0.1))
 hlines!([0], color=:gray)
 f
 
 # And some topoplots: 
-positions = [Point2f(p[1] + 0.5, p[2] + 0.5) for p in to_positions(headmodel_HArtMuT.electrodes["pos"]')]
+positions = [
+    Point2f(p[1] + 0.5, p[2] + 0.5) for
+    p in to_positions(headmodel_HArtMuT.electrodes["pos"]')
+]
 
-df = UnfoldMakie.eeg_matrix_to_dataframe(mean(data, dims=3)[:, :, 1], string.(1:length(positions)));
+df = UnfoldMakie.eeg_array_to_dataframe(
+    mean(data, dims=3)[:, :, 1],
+    string.(1:length(positions)),
+);
 Δbin = 20 # 20 samples / bin
-plot_topoplotseries(df, Δbin; positions=positions, visual=(; enlarge=1, label_scatter=false))
+plot_topoplotseries(
+    df,
+    bin_width=20,
+    positions=positions,
+    visual=(; enlarge=1, label_scatter=false),
+)
 
 # ## 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=100);
-fig, ax, hm = heatmap(transpose(pvals))
+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
+loopvectorization